A comprehensive evaluation of differentially expressed mRNAs and lncRNAs in cystitis glandularis with gene ontology, KEGG pathway, and ceRNA network analysis.

BACKGROUND: Cystitis glandularis (CG) is a proliferative disorder of the urinary bladder characterized by transitional cells that have undergone glandular metaplasia. The underlying mechanism associated with this transformation is poorly understood.
METHODS: The expression of messenger RNA (mRNA) and long non-coding RNA (lncRNA) from normal bladder mucosa and CG were compared using microarray analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to describe molecular interactions.
RESULTS: Microarray analysis identified 809 significantly dysregulated mRNAs in CG tissues; 606 were up-regulated and 203 were down-regulated (greater than 2-fold difference in expression from normal tissue, P<0.05). KEGG pathway analysis showed that the mRNAs that co-expressed with lncRNAs were enriched in the cell cycle regulation pathway. Four up-regulated lncRNAs (ENST00000596379, ENST00000463397, NR_001446 and NR_015395) were identified in the coding-non-coding co-expression (CNC) network analysis as being associated with the expression of four mRNAs (SMAD3, ORC1, CCNA2 and CCNB2). NR_015395 was revealed to be a competing endogenous RNA (ceRNA) of miR-133a-3p that targets SMAD3.
CONCLUSIONS: This is the first work to measure the expression of dysregulated lncRNA and ceRNA in CG and identify the crosstalk between mRNA and lncRNA expression patterns in the pathogenesis of CG. 2020 Translational Andrology and Urology. All rights reserved.

Introduction

Cystitis glandularis (CG) is a proliferative disorder of the bladder mucosa characterized by transitional cell glandular metaplasia. Typical and intestinal types of CG have been described (1). CG has been occasionally demonstrated to be a locally aggressive benign lesion (2). While it is reported that CG was associated with adenocarcinoma of the bladder, the causal relationship has yet to be fully established (3).

The pathogenesis of CG is not well understood. Through immunohistochemical techniques and in situ hybridization, expression of P16 in samples with CG were found to be abnormal, with expression in normal bladder mucosa defined as normal (4). After treatment, expression levels of p53 and Ki-67 were significantly lower in the bladder mucosa of patients with CG, suggesting that p53 may be involved in the development of CG (5). A recent comparative analysis of RNA sequencing of CG tissues and normal control tissues revealed differential expression (DE) of multiple genes including CX3CL1, CXCL6 and CXCL1, suggesting that these genes may be associated with the pathogenesis of CG (6). However, these limited studies do not effectively assess a correlation between molecular markers and the development of CG.

The understanding of the cellular control of genetic expression is rapidly expanding. Recently, gene expression has been shown to be inhibited by small, non-coding RNAs called microRNAs (miRNAs). miRNAs are RNAs of approximately 20 nucleotides in length that undergo base-pairing with complementary sequences in messenger RNA (mRNA), resulting in cleavage of the mRNA, destabilization of the mRNA, and less efficient translation of the mRNA. Individual miRNA can target many different genes. This miRNA-mediated post-transcriptional regulation of gene expression is believed to be present in half of all genes (7).

One type of long non-coding RNAs (lncRNAs) can bind to miRNAs and alter their inhibitory effect on mRNA expression, which adds another layer of control to mRNA expression (8). These lncRNAs are also called competing endogenous RNAs (ceRNAs) (9-12). The biological function of lncRNA has received great attention in recent years. Although it does not encode proteins, lncRNA does participate in many important biological processes (BP) (13-15). For example, Montes et al. reported that lncRNA MIR31HG targets p16 expression to regulate senescence (16). Additionally, Liu et al. found that lncRNA-UCA1 was regulated by p27 protein expression in primary cardiomyocytes to promote apoptosis (17).

We previously reported that several non-coding RNAs have been discovered, and we also found that lncRNA UCA1 plays an important role in the progression of CG (18-22). In this study, we aimed to investigate the role of novel lncRNAs expression in patients with CG. We used microarray screening to analyze the expression profiles of lncRNAs and mRNAs in normal and CG bladder tissues. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to predict the biological roles and potential signaling pathways of these differentially expressed lncRNAs. A coding-non-coding co-expression (CNC) network and ceRNAs analysis were performed to explore the potential roles of differentially expressed lncRNAs in CG pathogenesis.

Methods

Source of specimens

Tissue from three cases of typical CG and three normal bladders (N) were obtained from the Department of Urology, Xiangya Hospital. Diagnosis was made by pathologic examination of biopsied tissue. No patient had any other urinary disease. All specimens were stored in liquid nitrogen immediately after surgical resection and were transferred to a freezer maintained at –80 °C.

RNA extraction

Tissue was homogenized using TRIzol Reagent (Invitrogen, Grand Island, NY, USA), and total RNA was extracted in accordance with the manufacturer’s instructions. RNA quantity and quality were measured using a Nano Drop ND-1000 Spectrophotometer (ThermoFisher Scientific, Wilmington, DE, USA). RNA integrity was assessed using standard denaturing agarose gel electrophoresis.

Microarray analysis

Sample labelling and array hybridization were applied in the light of the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology, Santa Clara, CA, USA), with minor modifications. In brief, mRNA was purified from total RNA after the removal of rRNA using an mRNA-ONLYTM Eukaryotic mRNA Isolation Kit (Epicentre, Madison, WI, USA). Each sample was then amplified and transcribed along the full length of the transcript into fluorescent complimentary RNA (cRNA) without 3' bias using random priming methods (Arraystar Flash RNA Labelling Kit, Arraystar, Rockville, MD, USA). The labeled cRNAs were purified using a RNeasy Mini Kit (Qiagen, Hilden, Germany). The concentration (µg cRNA) and specific activity (pmol Cy3/µg cRNA) of the labeled cRNAs were determined using a NanoDrop ND-1000 Spectrophotometer.

Labeled cRNA from each tissue specimen was fragmented according to a modified protocol. In brief, cRNA was placed in 10× blocking agent and 25× fragmentation buffer and incubated at 60 °C for 30 min. Then 2× GE hybridization buffer was added to dilute the labeled cRNA to stop the fragmentation reaction. Subsequently, hybridization solution was dispensed into the gasket slide, which was then assembled with the lncRNA expression microarray slide. The slides were incubated for 17 h at 65 °C in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed and scanned using an Agilent Technologies G2505C SureScan High-Resolution DNA Microarray Scanner (Agilent Technologies, Santa Clara, CA, USA). Microarray hybridization and expression data collection were performed by KangChen Bio-tech, Shanghai, China.

A two-fold change (FC) in expression between CG and normal tissues was considered significant. This was calculated as log2 (CG/N). All values greater than 0.58 were up-regulated and all values less than –0.58 were considered down-regulated genes.

Quantitative real-time PCR (qRT-PCR) validation

Total RNA was isolated from tissues to generate cDNA using the SuperScript III Reverse Transcriptase (Invitrogen, Grand Island, NY, USA). Each cDNA sample was amplified using a Biosystems ViiA 7 Real-Time PCR System. Each cycle consisted of denaturation at 95 °C for 10 s and extension at 60 °C for 60 s. β-actin was used as an endogenous control to normalize each sample. The primers are listed in .

Table 1

Primers designed for qRT-PCR validation of candidate mRNAs and lncRNAs

Gene name	Primer	Product length (bp)
β-actin (H)	F: 5' GTGGCCGAGGACTTTGATTG 3'	73
	R: 5' CCTGTAACAACGCATCTCATATT 3'
SMAD3	F: 5' TGGAGCTGACACGGAGACAC 3'	121
	R: 5' CGCTGGTTACAGTTGGGAGAC 3'
ORC1	F: 5' CAGAATGAAGCGACCTCTACTC 3'	173
	R: 5' GCCTCTTCTTCGTCACTGCT 3'
CCNA2	F: 5' CATACCTTAGGGAAATGGAGG 3'	185
	R: 5' GGAAGACAGGAACCTATCAATG 3'
CCNB2	F: 5' CAAAATCGAGGACATTGATAACG 3'	206
	R: 5' GCAGAAGCCTAAACTTGGAGTG 3'
ENST00000596379	F: 5' CTGCTTTCAAACAAAGAGGC 3'	120
	R: 5' AGAGGATAGGAGGGCACAAC 3'
ENST00000463397	F: 5' GTTGCCTAAGGAGGGGTGAA 3'	113
	R: 5' TCACAGATTCACGTGGGTTTTA 3'
NR_001446	F: 5' TCAACAACCAACCAGGAGC 3'	145
	R: 5' CCTCTGGTCTTCTAACATTTGC 3'
NR_015395	F: 5' TGGAGTGCCACCTATGACCT 3'	86
	R: 5' TTGTAGACATCCTCAGGCTCG 3'

qRT-PCR, quantitative real-time PCR; lncRNA, long non-coding RNA; mRNA, messenger RNA.

GO and KEGG pathways analysis

The GO project provides a controlled vocabulary to describe gene and gene product attributes in any organism (http://www.geneontology.org). The ontology covers three domains: BP, cellular components (CC) and molecular functions (MF). A Fisher’s exact test was used to the compare expression of CG and normal tissues of items found in the DE and GO annotation lists. P values less than 0.05 were considered statistically significant. KEGG pathways (http://www.genome.jp/kegg/) were used to describe molecular interactions, reactions, and relation networks of identified genes and gene products.

CNC network analysis

To identify interactions among the differentially expressed mRNAs and lncRNAs, we constructed a CNC network. The network was constructed according to the normalized signal intensities of specific mRNA and lncRNA expression levels. A Pearson correlation coefficient >0.95 was used as a cutoff to identify coding and non-coding genes.

ceRNA network analysis

lncRNAs and mRNAs from CG specimens with significantly different expression than that found in normal tissues were subjected to ceRNA network nalysis in order to identify potential miRNA response elements (MREs). Overlapping of the same miRNA seed sequence binding site on both lncRNAs and mRNA sequences was identified in order to predict potential interactions. Potential targets of miRNAs were predicted using our miRNA target prediction software that is based on TargetScan (http://www.targetscan.org/) & miRanda (http://www.mircorna.org/).

Results

Microarray analysis identified 436 significantly dysregulated lncRNAs in CG tissues. Among them, 198 were up-regulated and 238 were down-regulated (greater than 2-fold difference as compared to normal tissue, P<0.05). The top 20 most significantly dysregulated lncRNAs are listed in .

Table 2

Top 20 aberrantly expressed lncRNAs in microarray analysis

Sequence name	P value	Log2 (CG/N)	Regulation	Chrom	Strand	Relationship
ENST00000426185	0.017	8.642	Up	Chr6	+	Intergenic
T184284	0.017	8.582	Up	Chr2	–	Intergenic
ENST00000555864	0.009	7.714	Up	Chr15	+	Intergenic
NR_003062	0.011	7.358	Up	Chr1	–	Exon sense-overlapping
ENST00000597260	0.019	7.316	Up	Chr19	+	Intergenic
ENST00000568243	0.010	7.134	Up	Chr18	–	Intergenic
TCONS_00024710	0.016	6.945	Up	Chr16	–	Intergenic
T117110	0.013	6.601	Up	Chr15	+	Intergenic
NR_038940	0.010	6.312	Up	Chr10	+	Intronic antisense
T106052	0.024	5.993	Up	Chr14	+	Intergenic
NR_121189	0.010	–7.944	Down	Chr9	+	Intergenic
NR_015423	0.003	–6.901	Down	Chr9	–	Intronic antisense
ENST00000608890	0.021	–6.584	Down	Chr18	–	Intergenic
T204932	0.023	–6.468	Down	Chr2	+	Intronic antisense
NR_121188	0.004	–6.250	Down	Chr9	+	Intergenic
NR_121192	0.021	–6.169	Down	Chr9	–	Intronic antisense
uc003jox.1	0.006	–6.151	Down	Chr5	+	Intergenic
uc003tcp.1	0.011	–6.149	Down	Chr7	–	Intron sense-overlapping
uc022baa.1	0.005	–6.090	Down	Chr8	–	Natural antisense
TCONS_00018898	0.015	–6.086	Down	Chr10	–	Intergenic

Chrom, chromosome number; lncRNA, long non-coding RNA; CG, cystitis glandularis; N, normal.

A total of 809 significantly dysregulated mRNAs were identified in CG tissues. Among them, 606 were up-regulated and 203 were down-regulated (greater than 2-fold difference as compared to normal tissue, P<0.05). The top 20 most significantly dysregulated mRNAs are listed in .

Table 3

Top 20 aberrantly expressed mRNAs in microarray analysis

Sequence name	Gene symbol	P value	Log2 (CG/N)	Regulation	Chrom	Strand
NM_144947	KLK11	0.015	9.662	Up	Chr19	–
NM_182502	TMPRSS11B	0.011	8.979	Up	Chr4	–
NM_015596	KLK13	0.010	8.940	Up	Chr19	–
NM_016190	CRNN	0.011	8.800	Up	Chr1	–
NM_014058	TMPRSS11E	0.013	8.771	Up	Chr4	+
NM_182606	TMPRSS11A	0.020	8.729	Up	Chr4	–
NM_001166034	SBSN	0.010	8.690	Up	Chr19	–
NM_002274	KRT13	0.012	8.558	Up	Chr17	–
NM_000359	TGM1	0.019	8.539	Up	Chr14	–
NM_005416	SPRR3	0.017	8.319	Up	Chr1	+
NM_002964	S100A8	0.042	8.264	Up	Chr1	–
NM_006783	GJB6	0.013	8.134	Up	Chr13	–
NM_001010909	MUC21	0.010	8.017	Up	Chr6	+
NM_198485	TPRG1	0.013	8.017	Up	Chr3	+
NM_002965	S100A9	0.028	7.918	Up	Chr1	+
NM_006061	CRISP3	0.023	7.859	Up	Chr6	–
NM_001944	DSG3	0.017	7.784	Up	Chr18	+
NM_020672	S100A14	0.009	7.770	Up	Chr1	–
NM_002974	SERPINB4	0.009	7.618	Up	Chr18	–
NM_006198	PCP4	0.003	7.658	Down	Chr21	+

Chrom, chromosome number; mRNA, messenger RNA; CG, cystitis glandularis; N, normal.

The variation of lncRNA and mRNA expression between CG and normal tissues is displayed in a Volcano plot (). Cluster analysis was performed to demonstrate relationships among lncRNA and mRNA expression patterns in normal bladder mucosa and CG ().

Figure 1

The expression microarray profiles of differentially lncRNAs (DE-lncRNAs) and coding RNAs (DE-mRNAs) in samples from three pairs of CG tissues and normal tissue. (A) DE-lncRNAs and DE-mRNAs are identified on a Volcano plot. Expression was plotted as P values of significance difference of expression (Y-axis) versus fold-difference in expression (X-axis). Horizontal green line marks a P value of 0.05. Vertical green lines mark a 2-fold increased or decreased expression in comparison to normal bladder mucosa; (B) heat map showing the expression profiles of lncRNAs and mRNAs. There were 436 differentially expressed lncRNAs and 809 differentially expressed mRNAs identified. The maps are based on the expression values of all expressed lncRNAs and mRNAs detected by microarray analysis and correspond to normalized expression values of significantly changed lncRNAs and mRNAs with 2-FCs in expression relative to normal tissue (P<0.05). Results from each tissue sample are presented in columns (normal 1–3 and CG 1–3). Red and green lines indicate high and low expression, respectively, of different RNAs. The intensity of expression is shown in the Color Key and Histogram. CG, cystitis glandularis; FC, fold change; DE, differential expression; lncRNA, long non-coding RNA; mRNA, messenger RNA.

GO and KEGG pathway analysis

Over- and under-expressed mRNAs that co-expressed with each of the differentially expressed lncRNAs were identified and listed as genes or gene products in the GO project (). The most enriched GO terms were epidermis development (ontology: BP, GO: 0008544), extracellular exosome (ontology: CC, GO: 0070062) and protein binding involved in cell-cell adhesion (ontology: MF, GO: 0098632). The most highly enriched GO sequences for the down-regulated transcripts were muscle system process (ontology: BP, GO: 0003012), contractile fiber part (ontology: CC, GO: 0044449) and carbohydrate derivative binding (ontology: MF, GO: 0097367).

Figure 2

GO and pathway analysis. The GO analysis included three domains: BP, CC, and MF. (A) The 10 most increased or decreased terms in each domain are presented; (B) KEGG enrichment analysis of top 10 enrichment scores for differentially expressed mRNAs. The vertical axis represents the pathway category and the horizontal axis represents the –log10(P value) of these significant pathways. The –log10(P value) denotes enrichment score representing the significance of GO term enrichment among differentially expressed genes. GO, gene ontology; BP, biological processes; CC, cellular components; MF, molecular functions; KEGG, Kyoto Encyclopedia of Genes and Genomes; mRNA, messenger RNA; DE, differential expression.

KEGG pathway analysis indicated that the mRNAs co-expressed with lncRNAs were involved in the regulation of cell cycle. The top 10 KEGG pathways are listed in .

Sequences that affect cell cycle progression were found in microarray analyses to be dysregulated in CG. Four related mRNAs, mothers against decapentaplegic homolog 3 (SMAD3), origin recognition complex subunit 1 (ORC1), cyclin A2 (CCNA2), and cyclin B2 (CCNB2), were selected for verification by qRT-PCR ().

Figure 3

CNC network analysis. (A) qRT-PCR was used to validate the up-regulated DE of four mRNAs found on microarray analysis; (B) the CNC network included sequences with a Pearson’s correlation coefficient ≥0.7. Green nodes represent dysregulated lncRNAs; yellow nodes represent four mRNAs: SMAD3, ORC1, CCNA2 and CCNB2; solid lines between lncRNAs and mRNAs indicate a positive correlation; dotted lines between lncRNAs and mRNAs indicate a negative correlation. The network consists of four mRNAs and 37 related genes; (C) qRT-PCR was used to validate the up-regulated DE of four lnRNAs found in the microarray analysis. The Y-axis indicates the FC (log2 transformed) computed from the qPCR and microarray data. Gene expression is normalized to β-actin and presented as mean ± SD. All the mRNAs and lncRNAs were up-regulated in the CG tissues (*, P<0.05). CNC, coding-non-coding co-expression; FC, fold change; qRT-PCR, quantitative real-time PCR; DE, differential expression; lncRNA, long non-coding RNA; mRNA, messenger RNA; CG, cystitis glandularis; N, normal.

Pearson coefficient analysis revealed that expression of these four genes was associated with expression of 37 other genes. Co-expression network profile analysis demonstrated 41 network nodes and 71 connections among the four mRNAs and 37 differentially expressed lncRNAs (). There were 24 negative and 47 positive correlations within the network. Four up-regulated lncRNAs were randomly chosen for qRT-PCR validation (ENST00000596379, ENST00000463397, NR_001446 and NR_015395). As expected, the four lncRNAs were up-regulated in the CG tissues compared to normal tissues (greater than 2-fold, P<0.05) (), which is consistent with the microarray chip data.

ceRNA network

Different ceRNAs have been shown to compete for the same MREs, regulating different levels of expression for each ceRNA (8). To investigate whether the four lncRNAs (ENST00000596379, ENST00000463397, NR_001446 and NR_015395) and four mRNAs (SMAD3, ORC1, CCNA2 and CCNB2) share any common MRE binding sites, a ceRNA network in CG was constructed with these elements from the microarray data. NR_015395 is a ceRNA of miR-133a-3p that targets SMAD3. Pearson coefficient analysis showed that expression of these two elements was correlated with the expression of 53 other genes ().

Figure 4

A competing endogenous RNA network was constructed based on lncRNA/miRNA and miRNA/mRNA interactions. The network consists of one mRNA, one lncRNA, and 53 related genes. The network included sequences with a Pearson’s correlation coefficient >0.7. In this network, the edges represent sequence matching lncRNAs connect expression correlated mRNAs via miRNAs. lncRNA, long non-coding RNA; mRNA, messenger RNA; miRNA, microRNA; ceRNA, competing endogenous RNA.

Discussion

Cellular proliferation and glandular metaplasia were observed in the transformation of normal bladder mucosa to CG. These changes were caused by alterations in the regulatory mechanisms that control normal cellular maintenance and growth-related functions. We examined the expression of mRNAs and lncRNAs in normal bladder mucosa and CG in order to identify abnormally expressed genetic material associated with the development of CG.

Using a second-generation lncRNA microarray, 809 mRNAs with abnormal expression in CG were identified. Of these, 606 mRNAs were found to be significantly up-regulated in CG as compared to normal bladder tissue (FC ≥2.0, P<0.05). Four mRNAs were related to cell cycle progression, SMAD3, ORC1, CCNA2 and CCNB2, and had increased mRNA expression in CG. This increased expression was verified with qRT-PCR.

Several studies have shown that lncRNAs play a significant role in regulating cell cycle activity (23-27). Microarray analysis identified 436 lncRNAs abnormally expressed in CG. Of these, 198 were found to be significantly up-regulated in CG (FC ≥2.0, P<0.05). GO and KEGG pathway analyses were performed to predict the potential functions of differentially expressed lncRNAs. GO analysis identified lncRNAs involved in epidermis development, extracellular exosome and protein binding related to cell-cell adhesion. KEGG pathway analysis indicated that the dysregulated lncRNAs were significantly associated with altered expression of cell cycle pathways in CG.

A CNC network analysis of SMAD3, ORC1, CCNA2 and CCNB2 mRNA expression identified 37 differentially expressed lncRNAs associated with their expression. SMAD3 expression was correlated with 15 lncRNAs expression, ORC1 expression was correlated with 18 lncRNAs expression, CCNA2 expression was correlated with 15 lncRNAs expression, and CCNB2 expression was correlated with 23 lncRNAs expression. Four up-regulated lncRNAs (ENST00000596379, ENST00000463397, NR_001446 and NR_015395) that were identified in the CNC network analysis as being associated with the expression of these four mRNAs were selected randomly for qRT-PCR verification. All four lncRNAs were also significantly up-regulated in CG.

We suggest a ceRNA mechanism maybe involved in this binding because of the consistent expression trends, and there may also be some shared MRE sites. The results of ceRNA network for CG showed that NR_015395 was a ceRNA of miR-133a-3p that targets SMAD3. These findings suggest that altered lncRNA expression (ceRNA) may be involved in CG-associated signaling pathways. In future work, we plan to establish more molecular mechanism at the clinical, cellular, animal, and molecular levels to provide new experimental evidence for the pathogenesis of CG. Solving these key problems will help us understand the molecular mechanism underlying the occurrence and development of CG and provide a theoretical basis as well as drug targets for new prevention and treatment strategies to inhibit dysuria through interruption of the inflammatory response.

This study has some limitations. A small number of tissue samples were examined. The analysis method identifies first-order miRNA-mediated lncRNA-mRNA associations but does not identify other influences. This method does not verify targets regulated by translational repression because it is based on changes in the transcript level. Further validation of the identified associations is needed.

Conclusions

Collectively, microarray data was used to identify mRNAs, lncRNAs, and ceRNAs with significantly altered expression profiles in CG. These findings were validated using qRT-PCR. Dysregulated mRNA from BP, CC and MF were unveiled. A significant dysregulation of lncRNAs was also observed. While these findings provide novel insights to the latent role of lncRNAs in CG, further research is required to fully explain their role in inducing transitional cell glandular metaplasia.

The article’s supplementary files as