Skip to main content
Download PDF

Expression profile of long noncoding RNAs and comprehensive analysis of lncRNA-cisTF-DGE regulation in condyloma acuminatum

BMC Medical Genomics volume 17, Article number: 167 (2024) Cite this article

Abstract

Objective

To identify differentially expressed long noncoding RNAs (lncRNAs) in condyloma acuminatum (CA) and to explore their probable regulatory mechanisms by establishing coexpression networks.

Methods

High-throughput RNA sequencing was performed to assess genome-wide lncRNA expression in CA and paired adjacent mucosal tissue. The expression of candidate lncRNAs and their target genes in larger CA specimens was validated using real-time quantitative reverse transcriptase polymerase chain reaction (RT‒qPCR). Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used for the functional enrichment analysis of these candidate lncRNAs and differential mRNAs. The coexpressed mRNAs of the candidate lncRNAs, calculated by Pearson’s correlation coefficient, were also analysed using GO and KEGG analysis. In addition, the interactions among differentially expressed lncRNAs (DElncRNAs)-cis-regulatory transcription factors (cisTFs)-differentially expressed genes (DEGs) were analysed and their network was constructed.

Results

A total of 546 lncRNAs and 2553 mRNAs were found to be differentially expressed in CA compared to the paired control. Functional enrichment analysis revealed that the DEGs coexpressed with DElncRNAs were enriched in the terms of cell adhesion and keratinocyte differentiation, and the pathways of ECM-receptor interaction, local adhesion, PI3K/AKT and TGF-ß signaling. We further constructed the network among DElncRNAs-cisTFs-DEGs and found that these 95 DEGs were mainly enriched in GO terms of epithelial development, regulation of transcription or gene expression. Furthermore, the expression of 3 pairs of DElncRNAs and cisTFs, EVX1-AS and HOXA13, HOXA11-AS and EVX1, and DLX6-AS and DLX5, was validated with a larger number of specimens using RT‒qPCR.

Conclusion

CA has a specific lncRNA profile, and the differentially expressed lncRNAs play regulatory roles in mRNA expression through cis-acting TFs, which provides insight into their regulatory networks. It will be useful to understand the pathogenesis of CA to provide new directions for the prevention, clinical treatment and efficacy evaluation of CA.

Peer Review reports

Introduction

Condyloma acuminatum (CA), mostly caused by low-risk HPV6 or HPV11, is characterized by repeated recurrence [1] and rapid proliferation [2]. The mechanisms by which HPV promotes cell proliferation and inhibits apoptosis of host cells include the following: on the one hand, HPV E6/E7 protein reduces apoptosis and autophagy by inhibiting p53/Rb; on the other hand, it activates the PI3K/Akt/mTOR and p38 MAPK/ERK pathways to promote cellular proliferation [3]. By recognizing EGFR, the HPV E5 protein activates the PI3K/Akt/MEK/Erk1/2 pathway [4] and enhances the transgenic activity of E7 [5]. E6/E7 of low-risk HPV also exerts similar activity to that of high-risk E6/E7 and affects common downstream pathways [6]. However, the upstream regulatory mechanisms of these pathogenic pathways for CA remain unknown.

Long noncoding RNAs (lncRNAs), more than 200 nucleotides, are described as RNAs with little or no protein-coding capability. These lncRNAs are involved in almost all important regulatory processes, including chromosome silencing or remodelling, transcriptional or posttranscriptional processing and intranuclear transportation [7, 8]. However, the significance of the expression and function of lncRNAs and their probable mechanism in CA has not been illuminated.

We speculated that lncRNAs expressed differentially in CA might regulate downstream target genes and participate in the occurrence and development of CA. To verify this hypothesis, high-throughput RNA sequencing was used to detect the genome-wide expression levels of lncRNAs and mRNAs between CA and their paired adjacent normal mucosa. We further validated the candidate lncRNAs in a larger number of CA specimens using real-time quantitative reverse transcriptase polymerase chain reaction (qRT‒PCR). Moreover, bioinformatics for differentially expressed genes (DEGs) using Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation were performed to explore the probable regulatory mechanisms of candidate lncRNAs. Furthermore, we analysed the interactions among differentially expressed lncRNAs (DElncRNAs)-cis-regulatory transcription factors (cisTFs)-DEGs, focusing on cis-regulatory mechanisms. Studying lncRNAs in CA may be helpful to provide a novel understanding of the pathogenesis of CA and to explore potential intervention targets for CA.

Materials and methods

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 (2018-086 and 2022 -795). Written informed consent was obtained from all subjects. Female CA patients aged 20–40 years old with positive HPV6 or HPV11 genotyping were included. A total of 5 pairs of HPV-positive CA specimens (CA group) and their adjacent mucosal tissue (CON group) were used for high-throughput transcriptional sequencing. Another 40 CA specimens (CA group) and 15 HPV-negative vaginal mucosal tissues (CON group) were collected for further verification using qRT‒PCR. An HPV GenoArray Test Kit (Hybribo, China) was used to detect HPV genotyping. CA patients with tumors and other immunocompromised diseases were excluded.

RNA extraction and transcriptome sequencing

All samples were stored in liquid nitrogen immediately after sampling. Then, the TRIzol reagent (Invitrogen, CA, USA) was used for extraction of total RNA. The RNA amount and purity of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) and the RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0. The Epicenter RiboZero Gold Kit (Illumina, San Diego, USA) was used for deletion of ribosomal RNA. Then, the final cDNA library was established by reverse transcription of the poly(A)- or poly(A) + RNA following the protocol for the mRNA-seq specimen preparation kit (Illumina, San Diego, USA). The average insert size for the final cDNA library was 300 ± 50 bp. Then, the 2 × 150 bp paired-end sequencing (PE150) was performed using the Illumina HiSeq 2000 platform at LC-Bio Biotech, Ltd. (Hangzhou, China). The sequencing data were submitted to the GEO database with the accession number GSE172140.

Read alignment and DEG analysis

Several strategies were applied to obtain clean reads as follows: The adapter of mRNA was removed and reads > 17 bp were retained; The first three bases were cut off; The bases < Q20 were removed and reads > 16 bp were retained; The bases containing 30% bases < Q20 were removed; NN ( unsure bases in sequencing ), poly G and poly A were discarded. Clean reads were aligned to the human GRch38 genome by HISAT2 [9], allowing 4 mismatches. Based on the uniquely mapped reads, the StringTie (https://ccb.jhu.edu/software/stringtie) was used to perform expression level for all mRNAs from input libraries by calculating the fragments per kilobase of transcript per million fragments mapped (FPKM) for each gene. DEGs were identified using edgeR software [10] and screened with a cut-off of fold change (FC) \(\ge\) 2 or \(\le\) 0.5 and false discovery rate (FDR) < 0.05 as the thresholds to determine whether genes were differentially expressed.

LncRNA prediction

Four software programs of Coding Potential Calculators (CPC2) [11], ORF Length and GC Content (LGC) [12], Coding-Non-Coding Index (CNCI) [13] and Coding Potential Assessment Tool (CPAT) [14], were used to predict or identify credible lncRNAs and for subsequent analysis and processing. In brief, the workflow of lncRNA prediction and further analysis of lncRNA trans-/cis-targets is shown in Fig. 1. Meanwhile, four types of novel transcripts were also removed, as shown in Fig. 1.

Fig. 1
figure 1

Workflow of lncRNA prediction and analysis in this research

Full size image

Functional enrichment analysis

Functional enrichment of GO terms and KEGG pathways for all DEGs was performed using the KOBAS 2.0 server [15]. 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).

LncRNA cis -regulatory target

The threshold of colocation was set as 100 kb upstream and downstream of lncRNAs in the cis-regulatory relationship pair [16]. For screening the lncRNA target pairs, the Pearson correlation coefficient between lncRNA and mRNA colocalization was calculated, and the absolute value of the correlation coefficient was set as greater than 0.6 and P value < 0.01. Then, the intersection of the two datasets of colocalization and coexpression was used to obtain the cis target of lncRNA. Then, lncRNA cis-regulatory transcription factors (cisTFs) were filtered out according to a catalog of 1796 transcription factors (TFs) retrieved from HumanTFDB (http://bioinfo.life.hust.edu.cn/HumanTFDB#!/).

Coexpression analysis between TFs and DEGs

Next, we also calculated the Pearson correlation coefficient between cisTFs and DEGs, and the absolute value of the correlation coefficient between cisTF-DEG relationship pairs was also set as correlation coefficient > 0.95 and P value < 0.01. Annotated cisTF-target pairs were downloaded from three databases: Transcription Factor Target Gene Database (http://tfbsdb.systemsbiology.net/download), TRRUST v2 Database (https://www.grnpedia.org/trrust/) and encode_chip TF targets Database. Furthermore, we generated a regulatory network between lncRNAs, cisTF, and DEmRNAs using Cytoscape software [17]. In brief, the excel of DElncRNA_TF_Target was imported into the Cytoscape, and the Source node, Target node, Source Node Attribute and Target Node Attribute according to the four column name were selected to layout the network of DElncRNA-cisTF-DEGs.

RT‒PCR for quantification of candidate lncRNAs

Primers for candidate lncRNAs and GAPDH were synthesized by RiboBio (Guangzhou, China). A 10 µL reaction volume containing 5 µL of 2 × SYBR® Premix Ex Taq™ and 0.2 µL of 50 × ROX reference dye (TaKaRa) was used for the amplification reactions. The ∆Ct data were collected automatically, and the data of -∆∆Ct were calculated by -∆∆Ct = average ∆Ct of the negative control group-∆Ct of the treated group. The relative expression of a target gene was calculated using 2-∆∆Ct. Each treatment was performed in triplicate.

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 RPM (Reads per million mapped reads) for each gene: RPM = ExonMappedReads × 106 ÷ TotalMappedReads. Hierarchical clustering based on Euclidean distance was performed using the pheatmap package. Student’s t test was used to compare the expression of target genes between CA group and CON group. Other statistical results were obtained by R software, where a P value < 0.05 was considered statistically significant.

Results

Differential expression profiles of lncRNAs and mRNAs in CA specimens

As shown in Fig. 2a, b and a total of 10,343 known lncRNAs and 2154 novel predicted lncRNAs were detected (expressed in at least 1 sample with FPKM > 0). The Venn diagram in Fig. S1a shows the overlapping lncRNAs by the four methods of CPC2, LGC, CNCI and CPAT, and the pie chart in Figure S1b based on the lncRNA type distribution shows that 3609 were intronic lncRNAs (84%), 471 were antisense lncRNAs (11%), and 216 were intergenic lncRNAs (5%). The distribution of exon length and density of the length distribution of noncoding RNAs and protein-coding RNAs are displayed in Fig. S1c-d.

Principal component analysis (PCA) for the expression pattern in these ten samples showed that all expressed genes and DEGs in CA and CON samples were clearly separated (Fig. 2c and f), indicating that the gene expression profile was different between CA and CON samples. Next, the differences in lncRNAs and mRNAs with a cut-off of FC ≥ 2 or ≤ 0.5 between the CA and CON groups were displayed using volcano plot filtering (Fig. 2g and h). A total of 2553 mRNAs were differentially expressed, with 1074 upregulated and 1479 downregulated (Fig. 2g). Meanwhile, a total of 546 lncRNAs, 171 upregulated and 375 downregulated, were found to be differentially expressed between CA and CON (Fig. 2h). The top 15 significantly upregulated and downregulated lncRNAs are displayed using a heatmap in Fig. 2i.

Fig. 2
figure 2

Identification of differentially expressed lncRNAs in condyloma acuminatum. (a-b) Venn diagram showing the number of detected lncRNAs. The known and new predicted lncRNAs were detected (expressed in at least 1 sample with FPKM > 0). (c-f) PCA based on the FPKM value of all mRNAs, all lncRNAs, differentially expressed mRNAs and differentially expressed lncRNAs. The ellipse for each group is the confidence ellipse. (g-h) Volcano plot showing differentially expressed mRNAs and lncRNAs between the CA group and the CON group, with cut-offs of DEseq2. FDR < 0.05 and FC ≥ 2 or ≤ 0.5. (i) Heatmap showing the expression profile of the top 15 significantly upregulated and downregulated lncRNAs

Full size image

Functional analysis of DEGs and DEGs coexpressed with DElncRNAs using GO or KEGG annotations

The functions of most lncRNAs have not been investigated completely until now, and analysing differential mRNAs with GO functional analysis and KEGG pathway analysis might be helpful to recognize the role of lncRNAs. In the present study, GO terms of BP for upregulated DEGs in CA were involved in the cell cycle, cell division and mitotic cell cycle (Fig. 3a), and those for downregulated DEGs in CA were involved in cell adhesion and angiogenesis (Fig. 3b). Meanwhile, the top 10 pathways for upregulated and downregulated DEGs using KEGG analysis are displayed in Fig. S2c-2d, including cell cycle and ECM-receptor interaction, PI3K/AKT signaling and TGF-ß signaling pathways.

Furthermore, the Pearson correlation coefficient between DEGs and DElncRNAs was calculated based on the transcriptional sequencing data, and the DEGs coexpressed with the candidate lncRNAs with correlation coefficients > 0.90 (P < 0.01) were annotated by GO terms or KEGG. Figure 3c displays the DElncRNA enrichment by CA compared with the CON samples and the number of coexpressed DEmRNAs. Interestingly, GO and KEGG annotations showed that DEGs coexpressed with DElncRNAs were also enriched in cell adhesion, keratinocyte differentiation (Fig. 3d), and the pathways involved in ECM-receptor interaction, PI3K/AKT signaling, TGF-ß signaling and local adhesion (Fig. 3e). This result indicated that these DElncRNAs in CA tissues might play potential regulatory roles in these pathways.

Fig. 3
figure 3

Analysis of differential expression of lncRNAs and mRNAs in condyloma acuminatum. (a-b) Bar plot showing the most enriched GO terms (biological process) for the upregulated and downregulated DEGs. (c) Scatter plot showing DElncRNA enrichment by CA compared with the CON samples and the number of coexpressed DEmRNAs. Red points represent upregulated lncRNAs involved in coexpression pairs, and blue points represent downregulated lncRNAs. Cut-offs of P value < 0.01 and Pearson coefficient > 0.9 were applied to identify the coexpression pairs. (d-e) Bubble diagram exhibiting the most enriched GO terms (biological process) and KEGG pathways of the DEmRNAs coexpressed with DElncRNAs

Full size image

Construction of the regulatory network among DElncRNAs-cisTF-DEGs in CA

LncRNAs play critical roles in gene regulation through mechanisms in cis or in trans [18]. One of our objectives in the present study was to identify potential cis-regulation of DElncRNAs. Figure 4a shows a positive relationship between DElncRNAs and cis-regulatory DEGs, whose levels are displayed using a heatmap in Fig. 4b. Next, we analysed the GO term enrichment and KEGG of these cisDEGs by lncRNAs and found that they were mainly enriched in transcription-related biological functions (Fig. 4c) and that most cisTFs belonged to the HOX family (Fig. 4d). Thus, we speculated that a large number of TFs, especially those in the HOX family, may be regulated by lncRNAs in CA. Through the intersection of these cisTF potential target genes from the database and cisTF coexpressed DEGs in our data, 95 DEGs were identified as possible targets of these cisTFs (Fig. 4e). The GO term annotation analysis for these 95 DEGs revealed that these genes were mainly enriched in GO terms of epithelial development, regulation of transcription and regulation of gene expression (Fig. 4f). In addition, a regulatory network among DElncRNAs-cisTF-DEGs was constructed, as shown in Fig. 4g.

In Fig. 5a-b, the expression of lncRNAs and cisTFs from the DElncRNA-cisTF-DEG network, as shown in Fig. 4g, is displayed using a heatmap. Based on the number of lncRNA-cisTF-DEG pairs and the probable functions of target lncRNAs or their cisTFs, we selected several lncRNAs of EVX1-AS, HOXA11-AS and DLX6-AS and showed the read distribution for these target lncRNAs and their cisTFs (Fig. 5c and e and S3a). Furthermore, based on the FPKM value, the expression of EVX1-AS, its cisTF HOXA13, and its target DEGs EMP1 and LEXM all decreased in the CA group (Fig. 5d). The expression of HOXA11-AS, EVX1 and PDGFRA was also attenuated in the CA group (Fig. 5f), while that of DLX6-AS and DXL5 was increased in the CA group (Fig. S3b). Then, the expression levels of these target lncRNAs and cis-regulatory DEGs were validated in 40 cases of CA compared with 15 cases of normal vaginal mucosal tissues using real-time qRT‒PCR. The qRT‒PCR results revealed a significant decrease in EVX1-AS, HOXA13, EMP1, LEXM, HOXA11-AS, EVX1 and PDGFRA and an obvious increase in DLX6-AS and DXL5 in CA tissues compared to CON tissues (Fig. 6a-c), which was in agreement with the RNA-seq analysis. Both our sequencing data and the qRT‒PCR verification results indicated the probable cis-regulatory mechanisms of lncRNAs-cisTFs-DEGs in CA pathogenesis.

Fig. 4
figure 4

Cis-regulatory TFs of DElncRNAs associated with condyloma acuminatum. (a) Scatter plot showing the log2 FC of DElncRNAs in CA compared with CON samples and their cis-regulatory genes. (b) Heatmap showing lncRNAs and cis-regulatory genes. (c) Bar plot showing the most enriched GO terms (biological process) of the differential lncRNA cis-target DEGs. (d) Bar plot showing the number of lncRNA-cisTF pairs. (e) Venn diagram showing the number of overlapping pairs of cisTF-target pairs from the database and cisTF-DEG coexpression pairs in the present study. (f) Bar plot showing the most enriched GO terms (biological process) of the cisTF target DEGs. (g) Network diagram showing the cisTF-targeted DEGs regulated by cisTF and DElncRNAs

Full size image
Fig. 5
figure 5

Expression of DElncRNAs, cisTFs and targeted DEGs in CA and CON samples. (a-b) The heatmap showing the expression profile of the DElncRNAs and cisTFs in Fig. 4g. (c-d) The read distribution showed lncRNA EVX1-AS and its cis-regulated TF HOXA13. The upward and green reads in the upper part represent for the red transcripts in the bottom part, and the downward and yellow reads represent for the blue transcripts in the bottom part. Boxplot showing the expression of lncRNAs, cisTFs and target DEGs. *** P < 0.001. (e-f) Read distribution showing lncRNA HOXA11-AS and its cis-regulated TF EVX1. The upward and green reads in the upper part represent for the red transcripts in the bottom part, and the downward and yellow reads represent for the blue transcripts in the bottom part. Boxplot showing the expression of lncRNAs, their regulated cisTF and target DEGs. **P < 0.01, *** P < 0.001

Full size image
Fig. 6
figure 6

Validation of target lncRNAs and their cis-regulatory TFs and DEGs in condyloma acuminata. (a-c) Bar plots showing the relative expression levels of the lncRNAs, cisTFs and TF-target DEGs using qRT‒PCR. Error bars represent the mean ± SEM. * P < 0.05, ** P < 0.01, *** P < 0.001

Full size image

Discussion

LncRNAs are a novel class of regulatory noncoding RNAs that participate in almost all biological or pathological processes [18]. Studies on HPV-related lncRNAs have mainly focused on the pathogenesis, invasion and metastasis of cervical cancers (CCs) and head and neck squamous cell carcinoma (HNSCC) caused by high-risk HPV. As early as 2016, Chinese scholars reported that in HPV18 + SiHa or HPV16 + HeLa cell lines, certain lncRNAs led to differential expression of their target mRNAs involved in several crucial processes of DNA repair, cellular cycle, proliferation and apoptosis [19]. In 2017, another Chinese group revealed the noncoding RNA profiles in HPV16-positive CC, compared them to adjacent tissue and constructed a coexpression and ceRNA network among differential lncRNAs. circRNAs, miRNAs and mRNAs [20]. In HNSCC, 132 different lncRNAs were identified under different HPV-infected states of HNSCC, among which HOTAIR, PROM1, CCAT1 and MUC19 were negatively correlated to myeloid-derived suppressor cell recruitment in HPV-associated HNSCC [21]. In recent years, the functions of several specific HPV-related lncRNAs have been identified, most of which are mediated by oncogenic HPV E6/E7 and promote cervical cancer progression [22,23,24]. Recently, the increase in lnc-FANCI-2 mediated primarily by E7 was shown to be due to YY1’s interaction with an E7 CR3 core motif and activation of the promoter of lnc-FANCI-2 [25]. In addition, HPV18 delayed replicative senescence of human keratinocytes through downregulation of antisense noncoding mitochondrial RNA-2 [26].

There have been very few studies on the transcriptome in CA specimens mainly caused by low-risk HPV. Dr. Zhang and coworkers compared the gene expression between the stable HPV-6bE7 HaCaT cell line and parental cells [27], and Gao et al. reported proteomic changes in CA tissue caused by mild hyperthermia treatment [28]. However, neither of them detected lncRNAs that were differentially expressed in CA. Although Tarkhan AH and coworkers reported lncRNA profiles in LR-HPV infection in 2022, their specimens were limited to common warts without a larger cohort for validation [29]. In the present study, we not only revealed transcriptional profiles, including lncRNAs and mRNAs but also analysed the probable regulatory pathway and cis-regulatory mechanism by lncRNAs in CA for the first time, further confirming the sequencing data in larger samples using qRT‒PCR. Herein, we used CPC2, LGC, CNCI [13] and CPAT [14] to predict or identify credible lncRNAs. CPC2 is a tool to evaluate the protein coding potential of transcript by comparing the sequence of the transcript with the known protein database [11]. LGC can characterize and identify lncRNAs based on the relationship between open read frame length and GC content, and can accurately distinguish lncRNAs from protein-coding RNA in a cross-species manner [12]. CNCI is a method to distinguish coding transcripts with non-coding ones by adjacent nucleotide triad features. The tool does not rely on known annotation files and can effectively predict incomplete transcripts and antisense transcripts. When the transcript score is < 0, it is noncoding [13]. CPAT is an analytical method to determine encoding and non-coding ability of transcripts by calculating Ficketter score and Hexmaer score based on ORF length and coverage [14]. In our analysis, a total of 10,343 known lncRNAs were detected and 2154 novel predicted lncRNAs were identified. Amongst them, a total of 546 lncRNAs were found to be dysregulated in CA specimens compared to adjacent mucosal tissues. In addition, PCA for all expressed mRNAs or lncRNAs and DEGs in CA separated from those in CON, indicating that HPV infection changed the expression pattern of host genes. Notably, the functional GO term and pathway analysis of DEGs coexpressed with DElncRNAs in CA were enriched in cell adhesion, keratinocyte differentiation, and the pathways involved in ECM-receptor interaction, PI3K/AKT signaling pathway, TGF-ß signaling pathway and local adhesion, which indicated that they played crucial roles in CA pathogenesis, consistent with our previous study [30].

For the lncRNA field, discerning functional noncoding RNAs from a vast transcriptome is a principal priority and challenge. It is a potential strategy to classify and elucidate lncRNA function based on the cis- or trans-regulatory pattern [18]. Cis-acting lncRNAs can be positioned and oriented relative to their nearby genes, such as lincRNAs around transcription factor start sites [31]. Furthermore, DEGs cis-regulated by these differential TFs (especially belonging to the HOX family) were enriched in GO terms of epithelial development and regulation of transcription or gene expression, which indicated the probable function of lncRNAs differentially expressed in CA.

We further confirmed the expression of 3 pairs of DElncRNAs and cisTFs, EVX1-AS and HOXA13, HOXA11-AS and EVX1, and DLX6-AS and DLX5. In our analysis, we chose the Pearson correlation coefficient > 0.6 and P value < 0.01 as set a cutoff to screen the lncRNA target pairs. Their PCR validation results were consistent with the sequencing data. For lncRNA EVX1-AS, there are very few data about its function. It was only reported that EVX1-AS might be correlated with colon cancer in the LncRNADisease database [32]. Although EVX1-AS and EVX1 have been reported to localize to the same genomic locus and are in fact transcribed from alternative promoters in close proximity [33, 34]. the Pearson correlation coefficient between EVX1 and EVX1-AS in our study was 0.59, and the P value was 0.071. Therefore, the pair of EVX1 and EVX1-AS was not included in our study. It is worth thinking about and confirm whether there is a better way to screen the pairs of lncRNA-cisTF and how to set the screen cutoff. Here, in our study, EVX1-AS and its cis-regulatory TF-HOXA13 both decreased in CA. HOXA13, a HOX transcription factor, usually functions as an oncogene [35]. In terms of HOXA13 in viral infection or inflammation, HOXA13 participates in inflammation activation. HBV DNA polymerase suppresses HBV replication via activation of the CREB1-HOTTIP-HOXA13 axis, followed by attenuation of HBV replication to result in persistent infection [36]. During influenza infection in the absence of the IFN-α/ß receptor, inflammatory and apoptotic responses may be initiated via induction of Ing1, Nr4a1, Polr2a, or Hoxa13 [37]. LncRNA-HOTTIP regulates the UV-mediated cellular response to UV through the coordinated activation of its neighboring gene Hoxa13 [38]. Thus, we speculated that attenuated EVX1-AS in CA might reduce the inflammatory response via cis-downregulated HOXA13 to allow HPV to evade host recognition. For lncRNA HOXA11-AS, a meta analysis showed up-regulated HOX11-AS in various cancers plays a tumor promoter [39], and indicates a risk factor for poor clinical outcome [40]. In the regulation of inflammation, HOXA11-AS mediated CaOx crystal-induced renal inflammation via the miR-124-3p/MCP-1 axis [41], while inhibition of HOXA11-AS repressed neuroinflammation and neuronal apoptosis in Parkinson’s disease model through through miR-124-3p-FSTL1-NF-κB axis [42] and miR-98-5p/EphA4 [43]. In addition, induction of HOXA11-AS by hypoxia/inflammation upregulated of HIF-1α and C/EBPβ further to promote epithelial-mesenchymal transition ability of nephroblastoma [44]. All these evidences indicate that HOXA-11 acts as a proinflammatory factor. In our results, whether the down-regulated HOX11-AS in CA also attenuates inflammation in CA to promote immune evasion, it needs further functional validation.

The expression of the pair of DLX6-AS and DLX5 increased in CA in the RNA-seq data and qRT‒PCR results. Both DLX6 and DLX-5 act as oncogenes to promote cellular proliferation and inhibit apoptosis [45,46,47,48]. Interestingly, DLX6-AS1 and DLX5 play their roles by targeting the PI3K/AKT/mTOR signaling pathway [49]. Whether increased DLX6-AS and DLX5 in CA participate in regulating the proliferation and apoptosis of HPV-infected cells needs to be confirmed by carrying out functional gain- and loss-of function experiments to explore the function of candidate lncRNAs in CA.

Conclusions

CA has a specific lncRNA and mRNA profile, and dysregulated lncRNAs in CA may be potential novel candidates. Functional enrichment analysis revealed that the DEGs coexpressed with DElncRNAs were involved in cell adhesion, keratinocyte differentiation, and the pathways involved in ECM-receptor interaction, PI3K/AKT signaling pathway, TGF-ß signaling pathway and local adhesion. We further constructed a network among DElncRNAs-cisTFs-DEGs and found that some differential lncRNAs play cis-regulatory roles on TFs, especially those belonging to the HOX family. In the future, we will explore the relationships between the candidate lncRNAs of EVX1-AS, HOXA11-AS, DLX6-AS and CA and their probable regulatory mechanisms. It will be useful to understand the pathogenesis of CA to provide new directions for the prevention, clinical treatment and efficacy evaluation of CA.

Data availability

Sequence data that support the findings of this study have been deposited in GEO database with the accession number GSE172140.

Abbreviations

CA:

Condyloma acuminatum

CNCI:

Coding-Non-Coding Index

CPC2:

Coding Potential Calculators

CPAT:

Coding Potential Assessment Tool

DEG:

Differentially expressed gene

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

LGCl:

ORF Length and GC Content

lncRNA:

Long noncoding RNAs

PCA:

Principal component analysis

RPM:

Reads per million mapped reads

RT-qPCR:

Real-time quantitative reverse transcriptase polymerase chain reaction

TF:

Transcription factor

References

  1. Le Poole C, Denman CJ, Arbiser JL. Immunosuppression may be present within condyloma acuminata. J Am Acad Dermatol. 2008;59:967–74.

    Article  PubMed  Google Scholar 

  2. Yang Z, Xiong H, Wei S, Liu Q, Gao Y, Liu L, et al. Yes-Associated protein promotes the Development of Condyloma Acuminatum through EGFR Pathway activation. Dermatology. 2020;236:454–66.

    Article  CAS  PubMed  Google Scholar 

  3. Zhang L, Wu J, Ling MT, Zhao L, Zhao KN. The role of the PI3K/Akt/mTOR signalling pathway in human cancers induced by infection with human papillomaviruses. Mol Cancer. 2015;14:87.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Kim SH, Juhnn YS, Kang S, Park SW, Sung MW, Bang YJ, et al. Human papillomavirus 16 E5 up-regulates the expression of vascular endothelial growth factor through the activation of epidermal growth factor receptor, MEK/ ERK1,2 and PI3K/Akt. Cell Mol Life Sci. 2006;63:930–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Ren S, Gaykalova DA, Guo T, Favorov AV, Fertig EJ, Tamayo P, et al. HPV E2, E4, E5 drive alternative carcinogenic pathways in HPV positive cancers. Oncogene. 2020;39:6327–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Oh ST, Longworth MS, Laimins LA. Roles of the E6 and E7 proteins in the life cycle of low-risk human papillomavirus type 11. J Virol. 2004;78:2620–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Statello L, Guo CJ, Chen LL, Huarte M. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol. 2021;22:96–118.

    Article  CAS  PubMed  Google Scholar 

  8. Anastasiadou E, Jacob LS, Slack FJ. Non-coding RNA networks in cancer. Nat Rev Cancer. 2018;18:5–18.

    Article  CAS  PubMed  Google Scholar 

  9. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  11. Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, et al. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45:W12–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Wang G, Yin H, Li B, Yu C, Wang F, Xu X, et al. Characterization and identification of long non-coding RNAs based on feature relationship. Bioinformatics. 2019;35:2949–56.

    Article  CAS  PubMed  Google Scholar 

  13. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41:e166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Li Y, Wang L. RNA coding potential prediction using alignment-free logistic regression model. Methods Mol Biol. 2021;2254:27–39.

    Article  CAS  PubMed  Google Scholar 

  15. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39:W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Yang P, Xu ZP, Chen T, He ZY. Long noncoding RNA expression profile analysis of colorectal cancer and metastatic lymph node based on microarray data. Onco Targets Ther. 2016;9:2465–78.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. 2018;172:393–407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Yang L, Yi K, Wang H, Zhao Y, Xi M. Comprehensive analysis of lncRNAs microarray profile and mRNA-lncRNA co-expression in oncogenic HPV-positive cervical cancer cell lines. Oncotarget. 2016;7:49917–29.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Wang H, Zhao Y, Chen M, Cui J. Identification of Novel Long non-coding and circular RNAs in human papillomavirus-mediated cervical Cancer. Front Microbiol. 2017;8:1720.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Ma X, Sheng S, Wu J, Jiang Y, Gao X, Cen X, et al. LncRNAs as an intermediate in HPV16 promoting myeloid-derived suppressor cell recruitment of head and neck squamous cell carcinoma. Oncotarget. 2017;8:42061–75.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Lai SY, Guan HM, Liu J, Huang LJ, Hu XL, Chen YH, et al. Long noncoding RNA SNHG12 modulated by human papillomavirus 16 E6/E7 promotes cervical cancer progression via ERK/Slug pathway. J Cell Physiol. 2020;235:7911–22.

    Article  CAS  PubMed  Google Scholar 

  23. He H, Liu X, Liu Y, Zhang M, Lai Y, Hao Y, et al. Human papillomavirus E6/E7 and long noncoding RNA TMPOP2 mutually upregulated Gene expression in Cervical Cancer cells. J Virol. 2019;93:e01808–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Barr JA, Hayes KE, Brownmiller T, Harold AD, Jagannathan R, Lockman PR, et al. Long non-coding RNA FAM83H-AS1 is regulated by human papillomavirus 16 E6 independently of p53 in cervical cancer cells. Sci Rep. 2019;9:3662.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Liu H, Xu J, Yang Y, Wang X, Wu E, Majerciak V, et al. Oncogenic HPV promotes the expression of the long noncoding RNA lnc-FANCI-2 through E7 and YY1. Proc Natl Acad Sci U S A. 2021;118:e2014195118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Villota C, Varas-Godoy M, Jeldes E, Campos A, Villegas J, Borgna V, et al. HPV-18 E2 protein downregulates antisense noncoding mitochondrial RNA-2, delaying replicative senescence of human keratinocytes. Aging. 2018;11:33–47.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Zhang B, Chen X, Zhou Q, Song Y, Sun S, Cheng H. Human gene expression microarray analysis of the HPV 6bE7-HaCaT stable cell line. Gene. 2018;657:60–8.

    Article  CAS  PubMed  Google Scholar 

  28. Sun YZ, Li JF, Wei ZD, Jiang HH, Hong YX, Zheng S, et al. Proteomic and bioinformatic analysis of condyloma acuminata: mild hyperthermia treatment reveals compromised HPV infectivity of keratinocytes via regulation of metabolism, differentiation and anti-viral responses. Int J Hyperth. 2019;36:383–93.

    Article  Google Scholar 

  29. Tarkhan AH, Al-Eitan LN, Alkhatib RQ, Alghamdi MA. Network analysis of long non-coding RNA expression profiles in common warts. Heliyon. 2022;8:e11790.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Liu X, Zhang Y, Wang S, Liu G, Ruan L. Loss of miR-143 and miR-145 in condyloma acuminatum promotes cellular proliferation and inhibits apoptosis by targeting NRAS. R Soc Open Sci. 2018;5:172376.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Gil N, Ulitsky I. Regulation of gene expression by cis-acting long non-coding RNAs. Nat Rev Genet. 2020;21:102–17.

    Article  CAS  PubMed  Google Scholar 

  32. Gao M, Guo Y, Xiao Y, Shang X. Comprehensive analyses of correlation and survival reveal informative lncRNA prognostic signatures in colon cancer. World J Surg Oncol. 2021;19:104.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Luo S, Lu JY, Liu L, Yin Y, Chen C, Han X, et al. Divergent lncRNAs regulate Gene expression and lineage differentiation in pluripotent cells. Cell Stem Cell. 2016;18:637–52.

    Article  CAS  PubMed  Google Scholar 

  34. Bell CC, Amaral PP, Kalsbeek A, Magor GW, Gillinder KR, Tangermann P, et al. The Evx1/Evx1as gene locus regulates anterior-posterior patterning during gastrulation. Sci Rep. 2016;6:26657.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Wen Y, Shu F, Chen Y, Chen Y, Lan Y, Duan X, et al. The prognostic value of HOXA13 in solid tumors: a meta-analysis. Clin Chim Acta. 2018;483:64–8.

    Article  CAS  PubMed  Google Scholar 

  36. Zhao X, Fan H, Chen X, Zhao X, Wang X, Feng Y, et al. Hepatitis B Virus DNA polymerase restrains viral replication through the CREB1/HOXA distal transcript antisense RNA homeobox A13 Axis. Hepatology. 2021;73:503–19.

    Article  CAS  PubMed  Google Scholar 

  37. Goodman AG, Zeng H, Proll SC, Peng X, Cillóniz C, Carter VS, et al. The alpha/beta interferon receptor provides protection against influenza virus replication but is dispensable for inflammatory response signaling. J Virol. 2010;84:2027–37.

    Article  CAS  PubMed  Google Scholar 

  38. Liang M, Hu K. Involvement of lncRNA-HOTTIP in the repair of Ultraviolet Light-Induced DNA damage in Spermatogenic cells. Mol Cells. 2019;42:794–803.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Li N, Yang M, Shi K, Li W. Long non-coding RNA HOXA11-AS in human cancer: a meta-analysis. Clin Chim Acta. 2017;474:165–70.

    Article  CAS  PubMed  Google Scholar 

  40. Lu CW, Zhou DD, Xie T, Hao JL, Pant OP, Lu CB, et al. HOXA11 antisense long noncoding RNA (HOXA11-AS): a promising lncRNA in human cancers. Cancer Med. 2018;7(8):3792–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Li Y, Yan G, Zhang J, Chen W, Ding T, Yin Y, et al. LncRNA HOXA11-AS regulates calcium oxalate crystal-induced renal inflammation via miR-124-3p/MCP-1. J Cell Mol Med. 2020;24(1):238–49.

    Article  CAS  PubMed  Google Scholar 

  42. Cao H, Han X, Jia Y, Zhang B. Inhibition of long non-coding RNA HOXA11-AS against neuroinflammation in Parkinson’s disease model via targeting mir-124-3p mediated FSTL1/NF-κB axis. Aging. 2021;13(8):11455–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Zhao L, Wang Z, Chen H, Du Y, Ma W, Tao Q, et al. Effects of lncRNA HOXA11-AS on Sevoflurane-Induced neuronal apoptosis and inflammatory responses by regulating miR-98-5p/EphA4. Mediators Inflamm. 2023;2023:7750134.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Zhu S, Zhou R, Tang X, Fu W, Jia W. Hypoxia/inflammation-induced upregulation of HIF-1α and C/EBPβ promotes nephroblastoma cell EMT by improving HOXA11-AS transcription. Heliyon. 2024;10(6):e27654.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Zhao H, Xu Q. Long non-coding RNA DLX6-AS1 mediates proliferation, invasion and apoptosis of endometrial cancer cells by recruiting p300/E2F1 in DLX6 promoter region. J Cell Mol Med. 2020;24:12572–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Huang Y, Ni R, Wang J, Liu Y. Knockdown of lncRNA DLX6-AS1 inhibits cell proliferation, migration and invasion while promotes apoptosis by downregulating PRR11 expression and upregulating miR-144 in non-small cell lung cancer. Biomed Pharmacother. 2019;109:1851–9.

    Article  PubMed  Google Scholar 

  47. Zhang TJ, Xu ZJ, Gu Y, Wen XM, Ma JC, Zhang W, et al. Identification and validation of prognosis-related DLX5 methylation as an epigenetic driver in myeloid neoplasms. Clin Transl Med. 2020;10:e29.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Zhang X, Bian H, Wei W, Wang Q, Chen J, Hei R, et al. DLX5 promotes osteosarcoma progression via activation of the NOTCH signaling pathway. Am J Cancer Res. 2021;11:3354–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. Zhang JJ, Xu WR, Chen B, Wang YY, Yang N, Wang LJ, et al. The up-regulated lncRNA DLX6-AS1 in colorectal cancer promotes cell proliferation, invasion and migration via modulating PI3K/AKT/mTOR pathway. Eur Rev Med Pharmacol Sci. 2019;23:8321–31.

    PubMed  Google Scholar 

Download references

Acknowledgements

We thank Wuhan Ruixing Biotechnology Co. Ltd for their support for data analysis.

Funding

This work was supported by the Natural Science Foundation of Zhejiang Province (LY22H110001) and the Medical Science and Technology Project of Zhejiang Province (2020KY913).

Author information

Authors and Affiliations

  1. Department of Urology, the First Affiliated Hospital, Zhejiang University School of Medicine, 79# Qingchun Road, Hangzhou, Zhejiang Province, 310003, China

    Bo Xie

  2. Department of Dermatology, the First Affiliated Hospital, Zhejiang University School of Medicine, 79# Qingchun Road, Hangzhou, Zhejiang Province, 310003, China

    Yinhua Wu, Su Wang, Liming Ruan & Xiaoyan Liu

  3. Department of Dermatology, Beilun People’s Hospital of Ningbo City, 1288# Lushan East Road, Ningbo, Zhejiang Province, 310058, China

    Liming Ruan

Authors
  1. Bo Xie
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yinhua Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Su Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Liming Ruan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiaoyan Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

X.L. and L.R. designed the study. B.X. and S.W. conducted the experiment. X.L. and B.X. wrote the manuscript. X.L., X.B., Y.W., S.W. and L.R. collected the clinical specimens and the data. All authors reviewed the manuscript.

Corresponding authors

Correspondence to Liming Ruan or Xiaoyan Liu.

Ethics declarations

Ethicsapproval and consent to participate

The clinical materials were collected after patients’ consent, and the study was approved by the Institutional Review Board (IRB) of the First Affiliated Hospital, Zhejiang University School of Medicine. Written informed consent was obtained from all subjects.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Xie, B., Wu, Y., Wang, S. et al. Expression profile of long noncoding RNAs and comprehensive analysis of lncRNA-cisTF-DGE regulation in condyloma acuminatum. BMC Med Genomics 17, 167 (2024). https://doi.org/10.1186/s12920-024-01938-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-024-01938-z

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us