Integrative analysis of microRNA and mRNA expression profiles in fetal rat model with anorectal malformation

Introduction

Anorectal malformations (ARMs) are among the most common congenital abnormalities with a wide clinical spectrum, occuring in approximately 1∕5,000–1∕1,500 live births (Endo et al., 1999; Van der Putte, 1986; Tsuda et al., 2011; Levitt & Peña, 2005). Recent studies demonstrated that ARMs are a group of complex malformations affected by many factors, and genetic factors have a great influence on the pathogenesis of ARM, such as Shh, Wnt, BMP-4, Gli2, Hox, Fgf, PCSK5, P63, and others (Wang, Li & Cheng, 2015; Wong et al., 2013; Khanna et al., 2018). Although ARM has been extensively studied, the detailed pathogenesis of ARM is still unclear; in particular, the molecular mechanisms that are involved in the ARMs remain uncharacteristic. It is generally known that non-coding RNAs (ncRNAs), are a family of RNA molecules that regulate gene expression but do not typically code proteins; ncRNAs participate in a variety of cellular processes including cell metabolic function and development (Mattick & Makunin, 2006; Thum & Condorelli, 2015; Godwin et al., 2010; Wang & Chang, 2011). NcRNAs can be subdivided into long ncRNAs (lncRNAs) with a length range from 0.2 and 2 Kb and small ncRNAs (<200 nucleotides long), which encompass microRNAs (miRNAs). In eukaryotic cells, miRNAs can regulate 31% of the target genes and participate in multiple processes of biological development and diseases (Lewis et al., 2013). However, there is a paucity of literature concerning the analysis of miRNA and mRNA co-expression in the development of ARMs. Hence, the purpose of our study was to integrate miRNA-mRNA expression profiles in the terminal hindgut tissue of ARM fetal rats at the twentieth day of the embryo (E20) with using Agilents miRNA and mRNA microarrays. In addition, the relevant data were used for bioinformatics analysis, which includes predicted target genes, GO enrichment analysis, and KEGG pathway analysis.

Materials and Methods

Results

miRNA and mRNA differential expression profiles

The miRNA expression profiles were significantly different in terminal hindgut tissues between the ARM and control groups. The process of the microarray and data analysis was shown as follows (Fig. 1). Of all the miRNA measured, 33 miRNAs were identified. Among the differently expressed miRNAs, eight miRNAs were up-regulated in the terminal hindgut tissue of the ARM group, while 25 miRNAs were down-regulated (>2-fold changes, P < 0.05; Table 1). The cluster analysis was shown in Fig. 2. In order to investigate the potential targets of the altered miRNAs, we also determined the mRNA expression profiles of the terminal hindgut tissues by using the Agilent microarray. Greater than 20,000 rat genes and transcripts were investigated and the ARM group and control group were compared. Among the 772 differentially expressed mRNAs, 350 of the genes were significantly up-regulated and 422 were down-regulated (≥2.0 fold changes, P < 0.05) in the ARM group compared with control group. The top 20 upregulated genes and the top 20 downregulated genes were listed in Table 2.

Table 1:

List of differentially expressed miRNAs between ARM and normal rats.

Upregulated miRNAs	P value	FC (abs)	Downregulated miRNAs	P value	FC (abs)
rno-miR-221-3p	0.0044	2.836	rno-miR-133a-5p	1.51E–05	128.189
rno-miR-10b-5p	0.0118	2.653	rno-miR-30c-2-3p	6.21E–04	45.897
rno-miR-326-3p	0.0129	40.909	rno-miR-144-5p	0.0024	43.322
rno-miR-183-5p	0.0201	2.305	rno-miR-543-3p	0.0034	67.572
rno-miR-3084a-3p	0.0259	13.513	rno-miR-136-5p	0.0037	40.949
rno-miR-200a-3p	0.0266	2.153	rno-miR-136-3p	0.0041	46.235
rno-miR-3559-5p	0.0352	2.116	rno-miR-376b-3p	0.0042	2.160
rno-miR-598-3p	0.0495	7.129	rno-miR-486	0.0053	3.118
			rno-miR-542-5p	0.0054	2.070
			rno-miR-1-3p	0.0093	10.819
			rno-miR-133a-3p	0.0143	3.407
			rno-miR-133b-3p	0.0146	6.663
			rno-miR-451-5p	0.0172	2.869
			rno-miR-347	0.0183	2.557
			rno-miR-1b	0.0192	2.778
			rno-miR-381-3p	0.0194	2.415
			rno-miR-503-5p	0.0199	2.719
			rno-miR-99a-5p	0.0209	2.147
			rno-miR-495	0.0248	2.758
			rno-miR-206-3p	0.0300	17.344
			rno-miR-411-5p	0.0353	2.031
			rno-miR-128-3p	0.0413	2.567
			rno-miR-431	0.0424	2.378
			rno-miR-122-5p	0.0475	212.683
			rno-miR-99b-3p	0.0490	12.287

DOI: 10.7717/peerj.5774/table-1

GO and KEGG pathway analysis of the predicted miRNA targets

The Gene Ontology (GO) system was used for functional enrichment analysis of the differentially expressed miRNAs predicted target genes. GO functional annotation of the predicted target genes of the differentially expressed miRNAs was shown in Fig. 3A. According to molecular function analysis, there were 166 GO terms significantly enriched (P < 0.05), and they were predominantly related to protein binding, identical protein binding, etc. in ARM development process. Analysis of the biological processes revealed the 622 significantly enriched GO terms (P < 0.05) to be predominantly associated with transcription, DNA-templated in ARM development process. Cellular component analysis identified 143 significantly enriched GO terms (P < 0.05) to be primarily related to cytoplasm, extracellular exosome in ARM.

Table 2:

Top 20 up- and downregulated mRNAs between ARM and normal fetal Rats.

Downregulated mRNAs	Fold change	P value	Upregulated mRNAs	Fold change	P value
Ky	87.915	0.0178149	Rpap1	40.312	0.0187
Fitm1	86.903	0.0158851	Alg2	21.131	6.17E–05
Col9a2	69.592	0.0101831	P2rx4	19.245	0.0133
Cited1	66.896	0.0132214	Duoxa2	14.719	0.0034
Cryab	59.710	0.0125412	Cyp4v3	13.758	0.0404
Serpinf1	52.046	0.0191088	RT1-A1	11.290	0.0312
Ablim3	50.569	0.0230989	Dhx32	11.101	9.78E-05
Gatm	46.246	0.0237851	Fgf9	9.546	0.0044
LOC685612	44.668	0.0173458	Gucy1b3	8.766	3.08E-04
Col7a1	43.956	0.0096485	Snrnp35	8.671	0.0021
Sgcg	41.283	0.0067426	RGD1562127	8.643	0.0016
Ttn	39.483	0.0055567	Kcnk2	8.621	0.0116
Myoz1	37.311	0.0189420	Aldh1b1	8.390	0.0175
Gfra4	33.696	0.0098202	Sntg1	8.082	8.35E-04
Csrp3	32.095	0.0056248	Mogat2	7.291	5.74E-04
Npepo	30.003	0.0205087	Slpi	7.255	0.0237
Clec4f	28.918	0.0107149	Adamdec1	6.929	3.22E-04
Epha3	28.519	0.0070704	Itgal	6.750	0.0342
Mybpc1	28.450	0.0113402	Madcam1	6.447	0.0352
Zbtb16	26.716	0.0063771	Tmem151a	6.351	0.0028

DOI: 10.7717/peerj.5774/table-2

The potential functional pathways of differentially expressed miRNA targets were then predicted by KEGG pathway system. The DAVID, and KEEG pathway system was used to analysis the microarray data. It showed that 64 pathways were predicted to be significantly related to ARM (P < 0.05). HIF-1 signaling pathway, Cytokine-cytokine receptor interaction, Pathways in cancer, GABAergic synapse, and cAMP signaling pathway were the main pathways involved in these predicted pathways (Fig. 3B).

GO and KEGG pathway analysis of the differentially expressed mRNAs

The Gene Ontology (GO) system was used for functional enrichment analysis of the differentially expressed mRNAs. The top 10 GO pathway enrichment terms for differentially expressed mRNAs were shown in Fig. 4A. According to molecular function analysis, there were 68 GO terms significantly enriched (P < 0.05), and they were closely related to heparin binding, actin binding, tropomyosin binding, etc. in ARM development process. According to biological processes analysis, it revealed that 331 significantly enriched GO terms (P < 0.05) to be predominantly associated with muscle contraction, muscle organ development, cardiac muscle contraction, and regulation of muscle contraction in ARM development process. Cellular component analysis identified 56 significantly enriched GO terms (P < 0.05) to be primarily related to extracellular space, myofibril, and troponin complex in ARM.

The KEGG pathways for the differentially expressed mRNAs were then performed by KEGG pathway system. The DAVID and KEEG pathway system were used to analysis microarray data. It revealed that 53 pathways were predicted to be significantly related to ARM (P < 0.05). Dilated cardiomyopathy, Hypertrophic cardiomyopathy (HCM), Cell adhesion molecules (CAMs), and ECM-receptor interaction were the main pathways involved in the development process of ARM (Fig. 4B).

Validation of miRNA and mRNA expression profile by RT-qPCR

To further investigate the credibility of our microarray based results, qRT-PCR was performed in 12 pairs of matched ARM/control fetal rat terminal hindguts. From the differentially expressed miRNAs and mRNAs observed by microarray, five miRNAs (down: rno-miR-1-3p, rno-miR-99b-3p, rno-miR-206-3p, Up: rno-miR-598-3p, rno-miR-3084a-3p) and five mRNAs (down: Fgf16, Wnt16, Cdh15, Up: Trim9, Ccl7) were chosen for RT-PCR validations. Relative expression levels of the selected miRNAs and mRNAs were depicted in Figs. 5A & 5B. The qRT-PCR results confirmed the accuracy of microarray findings. Compared with controls, rno-miR-598-3p and rno-miR-3084a-3p were significantly increased, while rno-miR-1-3p, rno-miR-99b-3p, and rno-miR-206-3p were significantly decreased in the terminal hindgut tissues of ARM fetal rat. Besides, qRT-PCR analysis revealed that the expressions of Trim9, Ccl7 were upregulated, whereas Fgf16, Wnt16, and Cdh15 were downregulated in the terminal hindgut tissues of ARM fetal rat, which is in line with the microarray data, too. Although the magnitude of changes differed between the two methods, these results demonstrated a high consistency between the microarray and RT-qPCR.

The regulatory network of miRNAs and target genes

The miRNA-target gene pairs via Cytoscape software were used to construct the miRNA-target genes regulatory network. Using the 142 miRNA-target gene pairs, a miRNA-target gene regulatory network was constructed (Fig. 6). In this network, rno-miR-128-3p, rno-miR-200a-3p, rno-miR-30c-2-3p, rno-miR-495, and rno-miR-543-3p, which regulate 15, 12, 13, 15 and 13 targets, respectively, demonstrated the highest connectivities. Whereas Igf1, Fut4, Kcnk2, Ngfr, and Nap1l2, which were regulated by eight, six, four, four and four miRNAs, respectively, were the mRNAs with the highest connectivities. According to the the network analysis, most of the identified transcripts were connected with several miRNAs, it suggested that multiple miRNA-mRNA interactions combinatorial effect in gene regulation by coexpressed endogenous miRNAs in ARMs.

Discussion

ARMs have been reported to be a multigene genetic disease. Both genetic and environmental factors are involved in the etiopathogenesis of ARMs. The spectrum of ARM phenotypes range from stenotic anus to cloacal malformation, and the etiology of ARMs are still unknown (Wang, Li & Cheng, 2015). The development of hindgut in the embryonic period is governed by multiple genes in the relevant signaling pathways. A result of genetic change in any stage of the hindgut development may lead to varieties of ARM phenotypes. Recent studies demonstrated that molecular and genetic facors were the contributing roles in its etiology. Only when the cascades of relevant genes expression are well orchestrated can the normal anorectum develop. Hence, further intensive study of ARMs may be benefical for illuminating the etiopathogenesis of this complex congenital malformation.

Studies of microRNAs (miRNAs) indicated that they are a series of noncoding RNAs with a length of 22 nucleotides, have the function of mediating silence and post-transcriptional regulation of gene expression (Ambros, 2004; Bartel, 2004), and are involved in many vital processes, including proliferation, cell differentiation, migration and apotosis (Hosako et al., 2009). The principle of miRNAs is that the miRNA usually interacted with the 3′ untranslated regions (3′-UTRs) of their targets mRNAs to negatively regulate gene expression at the post-transcriptional level (Saito et al., 2009). Up to now, there were more than 2,500 human miRNAs recorded in the miRBase (Kozomara & Griffiths-Jones, 2011; Kozomara & Griffiths-Jones, 2014), and research showed that more than 60% of protein-coding transcripts were predicted to be the targets of miRNAs and regulated by miRNAs (Fabian, Sonenberg & Filipowicz, 2010). However, there is paucity of literature related to the microRNA expression profiles and their biological effects in ARMs.

In the present study, we evaluated the terminal hindgut expression of miRNAs and mRNAs in both ETU-induced ARM rats and control rats using microarray analysis. We found 33 miRNAs and 772 mRNAs differentially expressed in hindgut tissues between the ARM fetal rats and the controls. GO terms and KEGG pathway analysis for the differentially expressed mRNA further revealed that the miRNA-mRNA contributed a lot to the development of ARMs. Seventy-six potential target genes for the 25 miRNAs were predicted by intersecting three bioinformatics databases and differentially expressed mRNAs, resulting in 146 miRNA-mRNA pairs. Morever, the expression of several miRNAs and mRNAs were validated by qRT-PCR and the outcomes were highly in accordance with the microarray data, which confirmed the credibility of our microarray data.

Among the 33 miRNAs that were differently expressed in the hindgut tissue of ARM fetal rats according to the miRNA microarray data, previously studies have reported that many of these miRNAs (rno-miR-221-3p, rno-miR-200a-3p, rno-miR-1-3p, rno-miR-133b-3p, rno-miR-133a-3p, rno-miR-486, rno-miR-451-5p) identified in our result were associated with many disease such as colorectal disease, prostate and pancreatic cancers (Sun et al., 2011; Qin & Luo, 2014; Kneitz et al., 2014; Sarkar et al., 2013). Kneitz et al. (2014) reported that miR-221 is an oncogenic miRNA which targeted CD117, and that it has the function of preventing cell proliferation and migration in endothelial cells. In this study, we also found that miRNA-200a-3p was significantly up-regulated; Pichler et al. (2014) reported that miRNA-200a-3p is able to exert its regulatory effect on EMT (epithelial to mesenchymal transition) and is involved in cancer stem cell properties in colorectal cancer.

miRNAs exert their function of regulating through inhibiting or degrading the translation of its targets. Hence, it is extremely important to identify miRNA target genes to clarify the development process of ARM. It is well known that the primary cause of ARMs was mainly attributed to the gene expression abnormalities. Gene mutations have been reported to be the most important contributing roles in the pathogenesis of ARMs. The normal development of hindgut will be disturbed by the mutations of the related genes, which thus lead to a series of ARM phenotypes.

Among the 772 mRNAs that showed to be significantly dys-regulated according to microarray datas, many of these mRNAs (Wnt10a, Wnt16, Fgf7, Fgf9, Fgf16, Fgf2, Bmp3, Hoxd1) identified in our result have been reported previously to be associated with the pathogenesis of ARMs. Based on GO analysis of biological processes and pathways, it revealed that the 772 significantly dys-regulated genes to be mainly relate to muscle contraction, muscle organ development, regulation of muscle contraction, skeletal muscle tissue development, response to drug, positive regulation of myotube differentiation, and these cellular events are closely related to extracellular space, myofibril, extracellular matrix, and contractile fiber in ARM development process. KEGG pathway enrichment analysis revealed that these differentially expressed mRNAs were potentially associated with the cell adhesion molecules (CAMs), ECM-receptor interaction, focal adhesion, PI3K-Akt signaling pathway, and etc. The PI3K-Akt signaling pathway is activated by many types of cellular stimuli or toxic insults and regulates fundamental cellular functions such as transcription, translation, proliferation, growth, and survival. The binding of growth factors to their receptor tryosine kinase (RTK) or G protein-couple receptors (GPCR) stimulates class Ia and Ib PI3K isoforms, respectively. PI3K-Akt signaling pathway containing many genes including Colla1, Fgf7, Fgf2, Fgf16, Fgf9, Igf7, Prkzc, etc., while some of those genes are also implicated in the development process of ARMs. These genes linked to the key pathways such as Shh, Fgf, and Wnt/β-catenin signaling pathways or were involved in the regulation of their downstream targets such as Bmp4, Fgf10, Gli2, Gli3, and Wnt5a to trigger the abnormal development of anorectum of ARMs.

In order to further investigate the role of miRNAs in ARMs, we integrated miRNA and mRNA expression profiles with miRNA target predictions by the miRWalk database to identify the positively and negatively correlated miRNA/mRNA pairs. Consequently, 25 miRNAs and 76 genes formed 146 miRNA-target gene pairs displaying reciprocal level of expression. The target prediction of eight miRNAs is not available in miRWalk databases. The negative correlation between miRNAs and their predicted target mRNAs expression identified in our study supports the hypothesis that miRNAs significantly regulate gene expression in the pathogenesis of ARMs. Although some miRNAs and their taget mRNAs did not display negative correlation, their regulatory effect might have been concealed by additional regulatory mechanisms or remain unclear. Hence, only the miRNAs and their predicted target mRNAs that displayed significant correlation with their cognate miRNAs expression were analysed acccording to our algorithm. In the present study, through detecting the predicted target genes of the differentially expressed miRNAs, we found that some important genes were involved in the pathogenesis of ARM, i.e., Fgf, Hox, Bmp, Wnt, and others. Recently, a study (Jin et al., 2017) showed that the miR-193 was significantly upregulated in the terminal hindgut tissues of ARM fetal rat, and miR-193 inhibited the expression of Hoxd 13 in ARM fetal rat by targeting Hox, while the normal expresion of Hoxd is absolutely essential to the development of the anorectum.

Fgf signaling is important during the early embryonic development in vertebrates (Böttcher & Niehrs, 2005; Gambarini et al., 1996); it plays a key role in patterning the gut tube through promoting posterior and inhibiting anterior endoderm cell fate (Kondoh, Kobayashi & Nishida, 2003). Fgf is expressed in epithelial cell, and has the function of patterning the gut tube by promoting posterior and inhibiting anterior endoderm cell fate (Gambarini et al., 1996). Fgf16 encoded a member of a family of proteins which possessed broad mitogenic and cell survival activities, and were associated with a variety of biological processes such as embryonic development, tumor growth and invasion, cell growth, morphogenesis, tissue repair. In addition, as predicted by bioinformatics databases, Fgf was the potential target gene of rno-miR-221-5p, rno-miR-133b-3p, rno-miR-381-3p, rno-miR-431, and rno-miR-495. The Wif1 (Wnt inhibitory factor 1) gene, which mediates the activity of Wnt (Wingless drosophila integration sites) signaling pathway, is a predicted target gene of the rno-miR-451-5p and rno-miR-495 according to TargetScan 7.1. Wif1 is an inhibitory Wnt (Wingless drosophila integration sites) that negatively regulates Wnt signaling. This gene is a member of the Wnt gene family, it plays a key role in the regulation of cell fate and patterning during embryogenesis (Jönsson & Andersson, 2001). The Wnt was reported to be a crucial signaling pathway and it has a paramount regulatory effect during the development of anorectum. Bambi (BMP and activin membrane-bound inhibitor) gene was significantly upregulated in terminal hindgut tissue of ARM fetal rats in our results, and is a predicted target gene of rno-miR-128-3p which was significantly downregulated. This gene encodes a transmembrane glycoprotein related to the type I receptors of the transforming growth factor-beta (TGF-β) family, whose members play important roles in signal transduction in many pathological processes, the encoded protein may function to limit the signaling range of the TGF-β family during early embryogenesis. Our results showed that there existed an obviously negative correlation between the expression of rno-miR-128-3p and Bambi during the development of ARMs.

The expression of rno-miR-221-3p, rno-miR-381-3p, rno-miR-128-3p, rno-miR-495, and rno-miR-221-3p were altered in the terminal hindgut in ARM fetal rats according to our results, and those miRNAs all involved in feed-forward loops that amplify or inhibit the Fgf, Bmp, Wnt signaling in the pathogenesis of ARMs. In the present study, a series of genes such as Fgf, Wnt that is necessary for embryonic development were low expressed, and they found to be related to the abnormal expression of rno-miR-221-3p, rno-miR-381-3p, rno-miR-128-3p, rno-miR-495, and rno-miR-221-3p in terminal hindgut tissue of ARM fetal rats. Therefore, we proposed that, in the process of ARMs, the dys-regulated expression of Fgf, Wnt mediated by relevant miRNAs, such as rno-miR-221-3p, rno-miR-381-3p, etc. would disturb the cytokine-cytokine receptor interaction, cell growth, cell differential during embryogenesis and promote abnormal development of anorectum.

It should be noted that because of the temporal and spatial expression pattern of mRNAs during the hindgut development in rat embryos with ETU-induced ARMs, there are some differences in our study in comparison with other studies. Also, some limitations existed in our study. On one hand, the expression of mRNAs are influenced by multiple factors. Both the upstream miRNA and the neighbour genes can affect the gene expression, and this condition is not taken into consideration in our work and can be integrated in future research. Also, the miRNA/mRNA that were predicted to have a potential effect on influencing the expression of pathogenic genes need further functional experiments for validation.

Conclusion

In this work, we identified the differentially expressed miRNAs and mRNAs in the terminal hindgut between ARM fetal rats and normal ones. We strongly suggest that miRNAs are paramount in regulating gene expression in the pathogenesis of ARM. The present work can be regarded as a new perspective and direction for future research on ARMs.

Supplemental Information