Identification of Differentially Expressed Genes in Pelvic Organ Prolapse by RNA-Seq.

BACKGROUND Pelvic organ prolapse (POP) brings major health issues for women, affecting 40% of postmenopausal women, and directly affects bladder and bowel function, as well as quality of life. In light of the projected growth in demand for care for pelvic floor disorders, determining the etiology and progression of POP has important public health implications. MATERIAL AND METHODS Uterosacral ligaments (USLs) samples of POP patients and normal controls were enrolled for RNA-Seq, and functional annotation analysis and Protein-Protein interaction (PPI) networks construction were performed for differentially expressed genes (DEGs). RESULTS A total of 81 DEGs were identified between POP and normal control, and distinctly classify all samples into normal and POP group by hierarchical clustering. Sixty-six DEGs demonstrated the same expression pattern among the POP samples with different stages. For those DEGs, canonical Wnt receptor signaling pathway was the most significantly enriched GO term (P value=3.33E-07), and neuroactive ligand-receptor interaction was the most significantly enriched pathway (P value=1.24E-03). In The PPI networks of 81 dysregulated genes, significant hub proteins contained TOP2A (Degree=54), KCNA5 (Degree=22) and PLA2G2A (Degree=19), suggesting their important role in the development of POP. CONCLUSIONS This RNA-seq analysis identified a POP signature of 81 genes, and some ECM-related genes, including COMP, NDP, and SNAI2 might participate in the pathology of POP and be applied as potential therapeutic targets.

Background

Pelvic organ prolapse (POP) develops as the loss of the support provided by pelvic fascia and muscles resulting in the bulging or herniation of the uterus, bladder, and/or rectum [1]. POP are major health issues for women, affecting 40% of postmenopausal women, and directly affects bladder and bowel function, as well as quality of life [2]. In light of the projected growth in demand for care of pelvic floor disorders over the coming decades [3], determining the etiology and progression of POP has important public health implications.

Given the recent studies into the etiology of POP, the pathophysiology of POP is multifactorial, including genetic predisposition, vaginal parity, connective tissue disorders, and factors associated with elevated intra-abdominal pressure (e.g., obesity [4] and chronic constipation) [5-7]. Recently published meta-analyses revealed the association between single nucleotide polymorphisms of COL3A1 and POP, providing overwhelming evidence of the genetic influences on POP [8,9]. However, these factors do not fully explain the etiology and progression of POP, due to some inconsistent evidence. For example, advanced POP has been observed in nulliparous women, and the end-stage disease often occurs several decades after the end of childbearing. Therefore, the underlying mechanism by which this failure occurs remains unclear.

The uterosacral ligaments (USLs) are condensations of endopelvic fascia, and contribute to primary uterine support [10], the disruption of structural components of these ligaments may lead to a loss of support and eventually cause POP. A decrease of the connective tissue extracellular matrix (ECM) proteins, such as collagen and elastic in USLs of POP in the form of quality and quantity, has been discovered [11,12].

High-throughput detection of global genomic expressions has allowed a systematic understanding of the complex biological processes of diseases. RNA-Seq, a next-generation sequencing technology, overcomes the limitation of conditional microarrays to measure global genomic expressions with high resolution and low cost [13,14]. The gene expression profiles between women with and without POP have been investigated by previous microarray studies, while RNA-seq approach has never been used to characterize changes of the POP transcriptome. Prolapse was most commonly defined as ≥ Stage II based on POP-Q assessment [15]. In the present study, we hypothesized that the global gene expression profile of POP with varying stages was somewhat different. To that end, we used the RNA-seq approach to compare the gene expression profiles of early- and advanced-stage POP (Stage II–IV) patients with those of normal controls to uncover the potential molecular mechanisms contributing to the pathogenesis of POP.

Material and Methods

Patients and tissues

The tissue specimens were obtained from the Affiliated Union Hospital of Fujian Medical University. Patients were defined as women with stage II to stage IV prolapse based on the POP quantification (POP-Q) examination, and controls were uterine leiomyomata patients without prolapse gynecological disorders or history of prolapse surgery. The mean age was 55 years for the control group and 61 years for the POP group. USLs specimens were excised during pelvic reconstructive surgery in patients or during gynecologic surgery for nonprolapse indications in controls, and immediately snap-frozen individually in liquid nitrogen and stored at −80°C. This study was approved by the Institutional Review Board at the Affiliated Union Hospital of Fujian Medical University on 5 June 2013, and the registration number was 2013KY018. All women provided informed consent.

RNA isolation and sequencing

Total RNA was extracted from USLs specimens using Trizol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions, and purified using an RNeasy mini kit (Qiagen, Manchester, UK). The quality and quantity of the purified RNA were determined using a Nanodrop ND-2000 spectrophotometer (Thermo Scientific, Wilmington, DE). RNA integrity was further verified electrophoretically through a 1.5% (w/v) agarose gel. Messenger RNA (mRNA) was purified from total RNA using oligo-d(T) probes for polyA selection with a TruSeq RNA library preparation kit (Illumina, San Diego, CA). The purified mRNA was subsequently fragmented into sizes of around 200 bp, and, after being subjected to complementary DNA (cDNA), the fragmented RNA was further converted into double-stranded (ds) cDNA. The short cDNA then underwent end repair, adapter ligation, and gel purification (2% TAE) to isolate 300nt fragments. The QIAquick PCR was then performed to assess the relative concentration of the library (Bioanalyzer 2100; Agilent Technologies). A HiSeqTM 2500 platform (Illumina) was used to perform sequencing.

Differential expression analysis

The raw sequencing data was filtered using SeqPrep () and Sickle () software to remove low-quality sequences, including ambiguous nucleotides and adaptor sequences. TopHat v1.3.1 was used to perform the alignment of cleaned sequencing reads to the UCSC human reference genome (build hg19).

The original alignment file produced was processed using Cufflinks v1.0.3 [16] to calculate the abundance of the transcripts, and Fragments per Kilobase of exon per Million fragments mapped (FPKM) was used to determine the transcription abundance of each gene. The reference gene transfer format (GTF) annotation file that was used for Cufflinks was retrieved from the Ensembl database (Homo_sapiens.GRCh37.63.gtf). Statistical analysis was carried out using paired t-tests. After applying Benjamini-Hochberg correction for multiple test, the false discovery rate (FDR) <0.05 was selected as the criteria for significant differences. Hierarchical clustering of DEGs was performed using the “pheatmap” function of the R/Bioconductor package [17].

Functional enrichment analysis of differentially expressed genes

Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were both used to investigate the biological function of those differentially expressed genes (DEGs), and these analyses was achieved by Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7, which is a web-based functional annotation tool4, and the cut-off of FDR was set at 0.05.

PPI Network Construction

Protein-protein interactions (PPIs) research can reveal the protein functions of DEGs at the molecular level. The identification of genome-wide protein interaction would help to interpret the cellular regulatory mechanisms [18]. We constructed a PPI network by using Biological General Repository for Interaction Datasets (BioGRID) (), and visualized the distribution characteristics of DEGs in the PPI network with Cytoscape [19]. In the PPI network, nodes stand for proteins, and edges stand for interactions between 2 proteins.

Results

Analysis of transcriptome sequencing

USLs specimens from different POP stages (stage II, III, and IV) were subjected to RNA sequencing. In total, 2.19×107, 2.10×107, 2.25×107 sequencing reads from POP specimens, and 1.89×107 reads from normal specimen were generated, and 1.90×107(86.89%), 1.82×107(86.58%), 1.94×107 (86.21%) reads from POP specimens, and 1.64×107 (86.65%) reads from normal specimen were uniquely aligned to the UCSC human reference genome (hg.19) (Table 1).

Table 1

RNA-Seq results.

Sample	Total reads	Mapped reads	Uniq reads	Mapped rate	Uniq rate
Normal	1.89E+07	1.73E+07	1.64E+07	91.58%	86.65%
Stage II	2.19E+07	2.01E+07	1.90E+07	91.78%	86.89%
Stage III	2.10E+07	1.93E+07	1.82E+07	91.57%	86.58%
Stage IV	2.25E+07	2.04E+07	1.94E+07	90.94%	86.21%

Differentially expressed genes

After applying a filter of FDR <0.05 cut-off value, genes that were significantly up- or down-regulated in the stage II, III, and IV in comparison with normal USLs were identified. As illustrated in Figure 1, the signature could distinctly classify all samples into normal and POP groups by hierarchical clustering. Within the stage II POP, 33 genes were significantly up-regulated and 16 genes were down-regulated. Similarly, 22 and 12 genes were up- and down-regulated, respectively, in the stage III POP. In the stage IV, 26 genes were significantly up-regulated and 9 genes were down-regulated. Twenty-six genes were specifically differentially expressed in the early-stage POP. Ten genes – CHRDL2, CPXM1, SALL1, KIAA1644, STAC2, CHODL, ZCCHC12, RPRM, PCP4, and LINC00890 – were consistently significantly differentially expressed during different POP stages). A total of 81 genes were differentially expressed, with 58 up-regulated and 23 down-regulated genes in USLs tissues of POP compared with those of normal control, and 66 genes demonstrated the same expression pattern among the POP samples with different stages.

Figure 1

The heat map of significantly up- or down-regulated genes in the stage II, III, and IV in comparison with the normal USLs.

Functional annotation

To improve biological interpretation of the transcriptional profiles, GO and pathway enrichment analysis of DEGs was performed among 3 distinct POP stages. The result showed that the most significantly enriched GO term was Wnt receptor signaling pathway (P value=2.67E-05) for DEGs in POP stage II in terms of biological processes, while for DEGs in POP stage III, developmental growth was significantly enriched (P value=9.06E-05), and for DEGs in POP stage IV, DNA catabolic process was significantly enriched (P value=1.03E-04) (Figure 2). For the 66 DEGs that demonstrated the same expression pattern among the POP samples with different stages, canonical Wnt receptor signaling pathway was the most significantly enriched GO term (P value=3.33E-07). When KEGG pathway analysis was performed, we found that DEGs of different POP stages all enriched significantly in neuroactive ligand-receptor interaction (P value=1.24E-03), consisting of genes such as GABRQ, AGTR2, PRLR, PRLHR, ADRA2C, and TACR2.

Figure 2A,B

GO analysis of DEGs in terms of biological processes in POP stage II and stage III. (A) Stage II; (B) Stage III.

PPI network

The PPI networks of 81 dysregulated genes in USLs tissues of POP were established, including 249 nodes and 263 edges. The significant hub proteins contained TOP2A (Degree=54), KCNA5 (Degree=22) and PLA2G2A (Degree=19) (Figure 3), suggesting their important role in the development of POP.

Figure 3

The PPI networks of 81 dysregulated genes in USLs tissues of POP.

Discussion

A series of changes takes place in USLs tissues from the onset to the development of POP at a molecular level. Along with our present work, previous studies have demonstrated significant differences in genome-wide gene expression between POP patients and normal controls when working with pelvic floor muscle samples, such as pubococcygeus muscle, levator ani muscle, and vaginal muscularis by high-throughput technology, exploring the pathogenesis of POP on a more systemic level [1,20-22]. This is the first RNA-seq study to analyze the transcription profiles of USLs tissues from POP, updating our knowledge of proximal uterine support defects at a molecular level. More specifically, we examined gene expression profiles of POP in different stages.

As mentioned above, we identified a large group of DEGs in the USLs tissues between POP and controls, distinctly classifying the samples into normal and POP groups. In addition, there were also significant differences in gene expression of different POP stages, which was consistent with our hypothesis. However, the results do not appear to support our hypothesis that higher POP resulted in more differences in gene expression profile. On the contrary, the early-stage POP (Stage II) resulted in more DEGs compared with advanced-stage POP (Stage III–IV). The explanation for this result may be that the onset of POP involved the activation of driver genes to induce the transition from normal to pathologic status, which need more pathogenic genes to participate in the process.

Interestingly, 10 genes – CHRDL2, CPXM1, SALL1, KIAA1644, STAC2, CHODL, ZCCHC12, RPRM, PCP4, and LINC00890 – consistently displayed significant differential expression during different POP stages, emphasizing their important roles in the maintenance of POP phenotype. The function of CHRDL2 remains unclear, and is expressed in heart and liver and hormone-dependent tissues, such as ovary, testis, and prostate. Interestingly, CHRDL2 contains cysteine-rich pro-collagen repeats, a specific feature of several other extracellular matrix proteins, and is involved in Insulin-like 3 (INSL3) signaling. For orl fetal rats characterized by inherited cryptorchidism, aberrant muscle pattern formation was observed in the fetal gubernaculum. In vitro experiments revealed enhanced activation of INSL3/RXFP2 signaling in the orl rats [23]. In our study, CHRDL2 was dramatically up-regulated in USLs tissues of POP, suggesting an important role in the process of POP development, especially in formation of altered muscle phenotype.

Previous studies have shown that connective tissue remodeling occurred in POP through regulation of ECM components [24,25]. In our POP samples, altered expression of other ECM-related genes was observed, such as COMP, NDP, and SNAI2. As a noncollagenous ECM protein, COMP is predominantly expressed in cartilage and tendons [26]. COMP could function as an organizer of the dermal collagen I network in healthy human skin. A recent study showed that increased COMP levels altered supramolecular architecture of collagen matrix in fibrotic skin pathologies. In our study, we found that COMP was dramatically down-regulated in USLs tissues of early-stage POP, which was also observed in a previous study [1], suggesting that down-regulation of COMP may promote initiation of POP through aberrant organization deposition of ECM protein [27].

To investigate the biological role of 66 DEGs displaying the same expression pattern among all POP samples, GO enrichment was performed. As a result, the canonical Wnt receptor signaling pathway was significantly enriched, including SOX4, WNT16, FZD5, SNAI2, and NDP. The canonical Wnt receptor signaling pathway plays a key role in a variety of biological processes, and is linked to many human diseases. Fibrotic-like diseases, arising from a defect in connective tissue regulation, are mediated by impairing the Wnt pathway in terms of the abnormal expression of β-catenin [28]. The biological process of immune response was also significantly enriched, including CXCL6, IGKV4-1, CXCL1, IGHG1, CCL19, IGHM, IGHG2, and IGHG4, which was in accordance with the findings of Brizzolara [22].

When KEGG pathway analysis was performed, we found that DEGs of different POP stages were all significantly enriched in neuroactive ligand-receptor interaction, consisting of genes such as GABRQ, AGTR2, PRLR, PRLHR, ADRA2C, and TACR2. Additionally, we found some genes involved in the regulation of uterine smooth muscle contraction, such as TACR2 and ADRA2C. TACR2 functions as receptors for tachykinins, and mediates tachykinin-induced uterine contractions at late pregnancy [29]. Compared with that of controls, there was a decrease in TACR1 (P<0.05), and TACR2 (P<0.0001) in muscle samples of acute diverticular disease, accompanied by increased collagen fibers between muscle bundles [30]. In our study, TACR2 was highly expressed in POP, stimulating intense contractions in uterine smooth muscle to rescue the prolapsed uterus. Similarly, ADRA2C, a subtype of Alpha-2-adrenergic receptors involved in the suppression of neurotransmitter and hormone release, and the contraction of smooth muscles, have a critical role in the regulation of cervical resistance in late pregnancy rat [31].

Conclusions

We used the RNA-seq approach to identify a POP signature of 81genes, and some genes, including COMP, NDP, and SNAI2, might participate in the pathology of POP and be applied as potential therapeutic targets. Components in neuroactive ligand-receptor interaction and canonical Wnt receptor signaling pathway also contribute to the pathogenesis of POP.