Alteration of RNA m6A methylation mediates aberrant RNA binding protein expression and alternative splicing in condyloma acuminatum

Patients and specimens

Approval for this study was obtained from the Institutional Review Board of the First Affiliated Hospital, School of Medicine of Zhejiang University (2022-795). Written informed consent was obtained from all subjects. The flow of the procedural steps in our study is shown in Fig. 1. CA female patients aged 20-40 years old with positive HPV genotyping were included, and those complicated with tumors, other immunocompromised diseases and other sexually transmitted diseases were excluded. A total of five pairs of HPV6 or 11-positive CA specimens (CA group) and their adjacent mucosal tissues (CON group) were collected for m6A immunoprecipitation and m6A sequencing. Furthermore, an additional 40 HPV-positive CA specimens (CA group) and 15 HPV-negative vaginal mucosal tissues (CON group) were collected for further verification using quantitative real-time PCR (RT-qPCR). The clinical characteristics of CA patients were shown in Table S1. All samples were removed using surgical scissors and were stored in liquid nitrogen immediately and kept in liquid nitrogen for not more than 1 month for further detection. Swab samples (Sarier et al., 2023) were obtained from the warts before excision to get the exfoliated cells, and DNA were delivered for HPV genotyping using a 23 HPV Genotyping Real-time PCR Kit (Hybribo, China), including 17 high-risk HPV types of HPV 16, 18, 31, 33, 35, 39, 45, 51, 52, 53, 56, 58, 59, 66, 68, 73, 82, and six low-risk HPV types of HPV 6, 11, 42, 43, 44, and 81.

M6A immunoprecipitation and m6A sequencing

Mature polyadenylated mRNAs were isolated from total RNA according to the GenElute™ mRNA Miniprep Kit (Sigma, Burlington, MA, USA) instructions and fragmented to a length of 100 nt using RNA fragmentation buffer. Then, m6A enrichment was performed by incubating Pierce™ ChIP-grade Protein A/G Magnetic Beads (26162; Thermo Fisher Scientific, Waltham, MA, USA) washed with IP buffer with 12.5 µg of anti-m6A antibody (202003; Synaptic Systems, Germany). Then the m6A antibody-binding RNA was eluted from the magnetic beads using m6A elution buffer, and the supernatant was collected for purification. Finally, the IP and Input samples were reversely transcribed to prepare cDNA libraries according to the protocol of the RNA-Seq Library Preparation Kit, and the quality-checked libraries were then subjected to PE150 sequencing using the Illumina HiSeq 2000 platform at LC-Bio Biotech, Ltd. (Hangzhou, China). The raw sequence data were deposited in the Gene Expression Omnibus (GEO) database under the accession number GSE172140.

Analysis of m6A-seq data

Cutadapt was used to remove adaptors and low-quality bases from raw reads, and reads less than 16 nt were discarded to obtain high-quality clean reads. The number of reads covered by the IP libraries was significantly higher than that of the Input, which could result in “peaks” and be used to indicate the methylation loci in the genome. The final differential peaks were called and screened through the ABLIRC pipeline (Xia et al., 2017). For each peak, the unique number of reads, the lengths of reads, and the observed maximum peak height were counted. Meanwhile, the computer simulations were used to generate reads with the same number and length for each gene. These outputting reads were mapped onto the same genes, to generate random maximum peak heights from overlapping reads. Finally, all observed peaks whose height is higher than the random maximum peak with a P <0.05 were selected. In addition, Homer software is used to show the shape of m6A peak.

Differentially expressed gene analysis

Clean reads were aligned to the human GRch38 genome by Tophat2 (Kim et al., 2013), allowing four mismatches. Based on the uniquely mapped reads, the fragments per kilobase of transcript per million fragments mapped (FPKM) for each gene were calculated to evaluate the gene expression levels. Differentially expressed genes (DEGs) were identified using edgeR software (Robinson, McCarthy & Smyth, 2010) and screened with a cutoff of FC (≥ 2 or ≤ 0.5) and FDR <0.05 as the thresholds to determine whether genes were differentially expressed.

Alternative splicing analysis

Alternative splicing analyses (ASEs) and RASEs between samples were defined and quantified by using the ABLas pipeline (Jin et al., 2017; Xia et al., 2017). Briefly, ABLas detected 10 major splicing types based on splice junction reads, including exon skipping (ES), alternative 5′ splicing site (A5SS), alternative 3′ splicing site (A3SS), intron retention (IR), mutually exclusive exons (MXE), mutually exclusive 5′ UTRs (5pMXE), mutually exclusive 3′ UTRs (3pMXE), box exons, A3SS&ES and A5SS&ES. To clarify the RBP-regulated ASEs, Student’s t-test was used to calculate the statistical P value after ASE detection for each sample, and the RASE ratio was defined based on the change in alternative reads and model reads of samples, whereas those with a significant change in splicing ratio and P <0.05 were considered to be significant ASEs regulated by RBP.

Coexpression analysis

The coexpression regulation patterns between RBP expression and RASEs were determined by calculating their Pearson correlation coefficients (PCCs) and corrected P values, and could be further classified as positive, negative and noncorrelation based on the PCC values.

Functional enrichment analysis

Functional enrichment of GO terms and KEGG pathways for all DEGs was performed using the KOBAS 2.0 server (Xie et al., 2011). The hypergeometric test and Benjamini −Hochberg FDR controlling procedure were used to define the enrichment of each term. Reactome pathway analysis was also performed online using the Reactome database (http://reactome.org).

Expression validation of DEGs and RASEs by RT-qPCR

RT-qPCR was conducted on an additional 40 HPV-positive CA specimens (CA group) and 15 HPV-negative vaginal mucosal tissues (CON group) to elucidate the validity of the RNA-seq data. Primer sequences of target DEGs are shown in Table S2 and were synthesized by RiboBio (Guangzhou, China). GAPDH was used as a reference gene. The PCR conditions were as follows: denaturing for 10 min at 95 °C, then 40 cycles including denaturation for 15 s at 95 °C, with annealing and extension for 1 min at 60 °C. Three technique replicates were performed for PCR amplification of each sample.

In addition, RT-qPCR assays were applied to analyze ASEs. A boundary-spanning primer and an opposing primer in a constitutive exon were used for the sequence encompassing the junction of constitutive or alternative exons, and detecting alternative isoforms, respectively. The boundary-spanning primer of alternative exons was designed according to “model exon” to detect model splicing or according to “altered exon” to detect altered splicing. The primer sequences for detecting ASE of RASAL2 and YWHAZ are also included in Table S2.

Statistical analysis

The R package factoextra was used for PCA, and the samples were clustered using the first two principal components. Next-generation sequencing data and genome annotations were visualized using in-house scripts (sogen) after normalization of the reads by FPKM for each gene. Hierarchical clustering based on Euclidean distance was performed using the pheatmap package. Student’s t test was used to compare AS ratio changes for target gene splicing events. Other statistical results were obtained by R software, where a P value <0.05 was considered statistically significant.

Distribution characteristics and functional analysis of m6A modification in CA and paired adjacent normal tissues

A total of 20 libraries consisting of five replicates for both input and IP samples in the CA and CON groups were sequenced (Table S3). The principal component analysis (PCA) with five pairs of IP and input samples after normalizing all gene expression levels is shown in Fig. 2A. The expression level of all genes clearly distinguished the samples in the four groups were clustered. Meanwhile, the heatmap clustering analysis of sample correlation was performed based on the normalized mapped reads on each gene, showing immunoprecipitated samples clustered (Fig. S1A). Most of input and IP reads of each sample were mapped to introns, CDS, 3′UTR, and intergenic regions by analysis of the distribution of unique mapped reads (Fig. S1B). When compared with the enrichment in input controls, m6A-IP reads were enriched around both the translation start and stop codons in both CA and CON samples (Fig. 2B). Notably, the increase in CON IP samples was much more pronounced than that in CA IP samples (Fig. 2B). This results indicated that HPV infection decreased the m6A level around the translation start and stop codons of mRNA in CA.

Then, we identified m6A peaks from m6A-seq data by using the ABLIRC pipeline with input reads as backgrounds (Li et al., 2019). Peaks detected in at least three CA or CON samples were retained, which resulted in 33,885 and 12,778 m6A peaks for CA and CON, respectively (Fig. S1C). A plot of the distribution of m6A peaks showed a preferable location in the intron region (54.7%) (Fig. S1D). A de novo motif search was performed by using the HOMER algorithm for the merged m6A peaks of CON or CA samples. Motif analysis showed that both the CON and CA samples had the canonical GGACU motif of human cells in the top third motif (Figs. S1E–S1F).

Concerning the function of CA- or CON-specific m6A peak-associated genes, GO analysis revealed that CON-specific m6A peak-associated genes were enriched in terms of DNA-dependent regulation of transcription and intracellular signal transduction (Fig. 2C), while CA-specific m6A peak-associated genes were enriched in terms of mitotic cell cycle, mitosis, stem cell maintenance and Notch signaling pathway (Fig. 2D). These results indicated that alteration of the m6A modification caused by HPV infection could participate in the regulation on cell cycle-associated genes to influence the cell proliferation or apoptosis of HPV-infected cells.

Differentially expressed genes with m6A modification between CA and paired adjacent normal tissues

The DEGs between CA and CON samples were further analyzed to explore the potential functional consequence of m6A alteration in the RNA. PCA for these ten samples showed that CA and CON samples were clearly separated (Fig. 3A), indicating that the gene expression profile was different between CA and CON samples.

To identify the DEGs between the CA and CON groups, edgeR was used with a cutoff of fold change ≥ 2 or ≤ 0.5 and a 5% false discovery rate (FDR). In total, compared to CON, there were 1,519 upregulated or 1,791 downregulated DEGs in CA tissues (Fig. 3B). Hierarchical clustering of normalized FPKM values of DEGs showed a clear separation between the CA and CON groups and a high consistency for the five replicate datasets (Fig. 3C), indicating that HPV infection significantly altered the transcription expression level in CA tissue.

Gene Ontology (GO) functional analysis was performed to annotate all 3,257 DEGs to explore the potential roles of DEGs. The upregulated genes in CA tissues were significantly enriched in terms of keratinization, mitotic cell cycle, inflammatory response, chemotaxis and cell adhesion, et al. (Fig. 3D). The downregulated genes were significantly enriched in terms of extracellular matrix organization, cell adhesion and innate immune response (Fig. 3E). It was suggested that HPV infection affects the expression of genes associated with the cell cycle and immune response in CA tissues.

Furthermore, CON- or CA-specific m6A peak-associated genes overlapped with the DEGs to reveal whether m6A alteration changes the expression of genes. It was shown that 626 CA- and 1030 CON-specific m6A peak-associated genes were significantly differentially expressed in CA tissues from CON tissues (Figs. 3F and 3H). Then, further GO analysis for DEGs with CON- or CA-specific m6A peaks uncovered that DEGs with CA-specific m6A peaks were also enriched in terms inducing the cell cycle, response to cytokines and cell adhesion (Fig. 3G), and the DEGs with CON-specific m6A peaks were enriched in terms inducing cell adhesion, angiogenesis and regulation of apoptotic processes (Fig. 3I). These results indicated that HPV infection would alter the expression of m6A-modified genes functioning in the cell cycle, cell adhesion and cellular apoptosis.

Alteration of the m6A level of RNA binding protein genes affects their expression in CA

RBPs play key roles in the regulation of HPV gene expression in HR-HPV-induced intracervical neoplasia and cervical cancer (Kajitani & Schwartz, 2022). Thus, we focused on CON- or CA-specific m6A peak-associated RBP genes that were differentially expressed between CA and CON samples. First, a list of 1,550 known RBPs was created from two previous studies (Baltz et al., 2012; Castello et al., 2012). Overlap analysis showed that 49 DEGs in our study were RBPs (Fig. 4A). The RBP genes expressed at an expression level of FPKM >0.5 in at least one sample were obtained. A heatmap plot of these RBP genes showed that ten RBPs (TLR7, MATR3, RNASEK, RBP7, TST, RBM20, RNASE4, PDCD4, ZCCHC24, BICC1, NSUN7 and TLR3) had lower expression in CA samples than in CON samples (Fig. 4B). However, compared to CON samples, the expression of EMG1, RBM38, ENOX1, BRCA1, RNF17, RBP1, ZC3H12D, OASL, PCBP3, DQX1, DZIP1L, AZGP1, ADARB1, DNMT3B, ZC3H12C, DNMT1 and RNASEH2A increased in CA samples (Fig. 4B).

Then, an overlap analysis was performed to identify the CON- or CA-specific m6A peak-associated differentially expressed RBP genes. The results showed that one differentially expressed RBP gene (DZIP1L) had a CA-specific m6A peak, and six differentially expressed RBP genes (ENOX1, ZC3H12D, DNMT1, CSDC2, DAZL, NR0B1) had a CON-specific m6A peak (Fig. 4C). DNMT1 and ZC3H12D had a m6A peak in CON samples (Figs. 5A–5B), while DZIP1L had a m6A peak in CA samples (Fig. 5C). Interestingly, all of these three RBP genes of DNMT1, ZC3H12D and DZIP1L showed higher expression in CA than CON samples (Figs. 5D–5F). To verify the expression of these three RBP genes with CON-specific or CA-specific m6A peaks, RT-qPCR was conducted to quantify their changes in additional CON and CA tissue samples in larger numbers. The RT-qRCR results revealed a significant increase in DNMT1, ZC3H12D and DZIP1L (Figs. 5G–5I) and a decrease in NR0B1, DAZL, ENOX1 and CSDS2 (Figs. S2A–S2D) in CA tissues compared to CON tissue, which was consistent with the RNA-seq data (Figs. 5D–5F, and Figs. S2E–S2H). These results suggested that m6A modification may be an important regulatory contributor to RBPs in CA tissues further influencing their protein expression.

Differentially expressed RBPs with altered m6A methylation potentially regulate AS in condyloma acuminatum

RBPs also play vital roles in the regulation of RNA processing, including AS of both host cellular RNA and HPV RNA (Fredericks et al., 2015; Nilsson, Wu & Schwartz, 2018). We wondered whether aberrant expression of RBPs caused changes in the AS pattern of some genes in CA tissue. In our data, a total of 1,703 RASEs were identified with 96 intron-retention AS events (IR) and 1,607 other types of AS events (NIR) (Fig. 6A). Most of the NIR events were A5SS (498), A3SS (372), cassette Exon (257) and Exon Skipping (254) (Fig. 6A).

Then, the percentage spliced in (PSI) value of these RASEs was calculated for each sample. A heatmap of the unsupervised hierarchical clustering of the PSI values of these RASEs showed that five CA samples and five CON samples could be clustered together (Fig. S3A). Genome location analysis showed that these RASEs occurred on a total of 1,224 genes (referring to RASGs). Overlap analysis of DEGs and RASGs showed that only 135 genes had changes in both expression and AS patterns between CA and CON samples (Fig. 6B). GO term analysis was conducted to explore the potential function of RASGs and showed that RASGs were enriched in the terms of cell–cell adhesion and extracellular matrix organization (Fig. 6C). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that RASGs were enriched in pathways including adherens junction, focal adhesion and NF-kappa B signaling pathways (Fig. S3B). These results indicated that HPV infection could change the AS pattern of host cellular RNA in CA tissue.

To explore whether these seven differentially expressed RBPs with CON-specific and CA-specific m6A peaks affect the AS in CON and CA samples, coexpression analysis was performed based on the expression values of RBPs and PSI values of RASEs in all the ten input samples. The RASGs with RASEs coexpressed with these seven RBPs (DZIP1L, ENOX1, ZC3H12D, DNMT1, CSDC2, DAZL, NR0B1) were extracted (correlation coefficient >0.7 and P value <0.01). Then, GO and KEGG analyses were performed to explore the potential functions of these RASGs. The results showed that these RASGs were significantly enriched in GO terms including cellular component disassembly involved in apoptotic process and regulation of apoptotic process (Fig. 6D). KEGG analysis showed that these RASGs were enriched in pathways including regulation of autophagy, ECM-receptor interaction and adherens junction (Fig. S3C). Furthermore, a coexpression network among the seven RBPs and RASEs, RASGs, and the top enriched GO terms of RASGs was constructed, as shown in Fig. 6E.

Finally, we focused on the three RBPs of DNMT1, ZC3H12D and DZIP1L that had been confirmed using RT-qPCR and their probable regulatory roles on in RASEs and RASGs. Interestingly, the top enriched GO terms of RASGs that might be regulated by DMNT1 were also involved in the regulation of apoptotic processes and cellular component disassembly involved in apoptotic processes (Fig. 7A). Then, some genes in the coexpression network with changes in AS patterns were presented as differences in PSI values among CA and CON samples. Our results showed that PSI values of A5SS events on RASAL2 and YWHAZ were higher in the CON group than in the CA group (Figs. 7B–7C), which was confirmed in additional CON and CA tissue samples using RT-qPCR (Figs. 7D–7E). Notably, m6A peaks on the same site of these three genes were detected for both the CA and CON samples (Figs. 7B–7C), which indicated that m6A methylation did not affect AS. In addition, we analyzed the expression level of DMNT1 in different groups with various clinical characteristics, and found that the more DMNT1, the more times or duration of treatment (Figs. S4A–S4B). It was suggested that highly expressed DMNT1 indicated a greater likelihood of relapse in the course of CA treatment. Taken together, these results indicated that differentially expressed DMNT1 with CON-specific m6A peaks could result in changes in the AS pattern of some genes involved in terms of the regulation of apoptotic processes, thus participating in the pathological process of CA.