Microarray analysis of miRNAs during hindgut development in rat embryos with ethylenethiourea‑induced anorectal malformations.

Anorectal malformations (ARMs) are one of the most common congenital malformations of the digestive tract; however, the pathogenesis of this disease remains to be fully elucidated. MicroRNAs (miRNAs) are important in gastrointestinal development and may be involved in the pathogenesis of ARMs. The present study aimed to profile miRNAs and examine their potential functions in rats with ethylenethiourea (ETU)‑induced ARMs. Pregnant Wistar rats (n=36) were divided randomly into ETU‑treated and control groups. The rats in the ETU‑treated group were gavage‑fed 1% ETU (125 mg/kg) on gestational day 10 (GD10), whereas the control group rats received a corresponding dose of saline. Embryos were harvested by cesarean section on GD14, GD15 and GD16. Hindgut tissue was isolated from the fetuses for RNA extraction and microarray analysis, followed by bioinformatics analysis and reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR) validation. Overall, 38 miRNAs were differentially expressed (all upregulated) on GD14, 49 (32 upregulated and 17 downregulated) on GD15, and 42 (all upregulated) on GD16 in the ARM group compared with the normal group. The top 18 miRNAs with |log2(fold change)| >4.25 were selected for further bioinformatics analysis. Among these miRNAs, five were differentially expressed at two time-points and were involved in ARM‑associated signaling pathways. The RT‑qPCR analysis revealed that three miRNA (miR), miR‑125b‑2‑3p, miR‑92a‑2‑5p and miR‑99a‑5p, were significantly differentially expressed in rats with ARMs compared with the normal group. In conclusion, the results suggested that the differential expression of miR‑125b‑2‑3p, miR‑92a‑2‑5p and miR‑99a‑5p during key time-points of anorectal formation in rats may have functions in the pathogenesis of ARM.

Introduction

Anorectal malformations (ARMs) are common congenital malformations of the digestive tract manifesting as ectopic anus, anal stenosis and/or an abnormal fistula between the rectum and urinary tract, with an incidence of 2-6/10,000 (1-4). ARMs are usually combined with malformations of other systems, and although the majority of patients with ARMs require surgery, their quality of life is likely to be impaired postoperatively (5). Although earlier studies have revealed that multiple genes and signaling pathways are important in the pathophysiology of ARMs (3,6), the specific mechanisms remain to be fully elucidated. Therefore, it is necessary to investigate the pathogenesis of ARMs.

MicroRNAs (miRNAs) regulate gene expression at the post-transcriptional level by binding to the 3′-untranslated regions (3′-UTRs) of target mRNAs, and have been reported to be important in various congenital diseases and in the embryonic development of various systems (7-10). Previous studies have investigated the functions of miRNAs in gastrointestinal development in various species. For example, miRNA (miR)-143 was reported to regulate genes involved in the differentiation of connective tissue cells and to induce the proliferation and differentiation of smooth muscle cells, thus influencing bovine gut development (8). As the main effector of miR-17-92 cluster components, a high expression of miR-17-5p in the crypt progenitor compartment is potentially involved in controlling cell differentiation and proliferation in human colon development (11). However, current understanding of the functions of miRNAs in the development of ARMs remains limited. An increased understanding of miRNA expression profiles may reveal their functional roles in ARM development, thus providing a valuable basis for investigating the mechanisms responsible for this condition.

Ethylenethiourea (ETU)-induced ARM rat fetuses are the most commonly used animal model for investigating this disease, and the key abnormalities in the experimental embryos include the following: i) no fusion between the urorectal septum and cloacal membrane; ii) delay of tailgut regression; iii) abnormal apoptosis in the cloacal wall; iv) underdevelopment of the dorsal cloaca and its membrane (12,13). In the present study, miRNA expression patterns in the hindgut of rat fetuses with ETU-induced ARMs were profiled between gestational day (GD)14 and GD16, which are critical time-points during anorectal development (4,12-14). Additionally, advanced bioinformatics analysis was performed, including target gene prediction, construction of regulatory networks and functional enrichment analyses to clarify the molecular mechanisms involved in the development of ARMs.

Materials and methods

Ethics statement

The present study was approved by the China Medical University Ethics Committee (Shenyang, China; no. 2015 PS213K) and all procedures involving animals were performed in accordance with the guidelines for the care and use of laboratory animals (15).

Animal models and tissue preparation

Mature female Wistar rats (n=36; age, 7-9 weeks-old; body weight, 250-300 g) were provided by the Experimental Animal Center of Shengjing Hospital of China Medical University (Shenyang, China) and housed in a specific-pathogen-free animal laboratory (room temperature, 22±2°C; humidity, 55±5%; 12 h light/dark cycle and ad libitum access to water and food) at The Key Laboratory of Health Ministry for Congenital Malformations (Shenyang, China). The surgical procedures were performed following sacrifice by intraperitoneal injection of sodium pentobarbital, and all efforts were made to minimize animal suffering.

ARM was induced in fetal rats as described in earlier reports (4). Briefly, 18 pregnant rats were gavage-fed a single dose of 125 mg/kg of 1% ETU (Sigma-Aldrich; Merck Millipore, Darmstadt, Germany) on GD10, whereas the remaining 18 pregnant rats received corresponding doses of ETU-free saline as a control. The rat fetuses were obtained by cesarean delivery on GD14-16, and the presence of ARMs was determined by light microscopy.

No malformations were observed in the 258 embryos of the normal rats. Among the ETU-treated embryos, all 236 embryos had a short or absent tail, and 14 died in utero. The incidence of ARMs in the ETU-treated embryos was 85.6% (202/236).

miRNA profiling and microarray analysis

Assessment of RNA quality

Total RNA was isolated using an miRNeasy Mini kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's protocol. RNA quality was assessed using a K5500 micro-spectrophotometer (Shanghai Drawell Scientific Instrument Co., Ltd., Shanghai, China) and formaldehyde agarose gel electrophoresis. Values of absorbance (A)260/A280 ≥1.5 and A260/A230 ≥1 indicated acceptable RNA purity, and an RNA integrity number ≥7 from the Agilent 2200 RNA assay (Agilent Technologies, Inc., Santa Clara, CA, USA) indicated acceptable RNA integrity.

miRNA microarray assay

miRNA profiling was performed using the RiboArray platform (Ribobio Co., Ltd., Guangzhou, China). In brief, l μg of total RNA was labeled with Cy3 using the ULS™ microRNA labeling kit (Kreatech; Leica Microsystems, Inc., Buffalo Grove, IL, USA) and hybridized on the microarray. The microarray contained 761 specific oligos for rat miRNAs, based on the Sanger miRBase 21.0 database (www.mirbase.org/).

Data normalization and analysis

Firstly, chip images were examined to exclude impurity, scratch marks and high background intensity. Following background adjustment by subtraction of background from foreground, the microarray data were normalized using a quantile (median) method and box-plots were generated based on the relative logarithm expression using the limma package in R software (version 3.3.0; www.r-project.org). Three-dimensional principal component analysis was performed using the scatter-plot3d package in R. Analysis of the differentially expressed miRNAs was performed using the limma package in R. Standard selection criteria for identifying the differentially expressed miRNAs were established as a fold change (FC) >2 or <−2, and P<0.05. Hierarchical clustering analysis of the miRNAs was performed and graphs were generated using the gplots package in R software.

Target gene prediction, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses

miRNAs with |log2(FC)| >4.25 were selected for further analysis. Target genes of these miRNAs were predicted using the target gene-prediction platforms miRWalk 3.0 (mirwalk.umm.uni-heidelberg.de/), TargetScan 7.1 (www.targetscan.org/vert_71/), miRBase 22 (www.mirbase.org/), miRDB (mirdb.org/miRDB/), and DIANA v5.0 (diana.imis.athena-innovation. gr/DianaTools/index.php?r=site/page&view=software). Those miRNAs predicted by at least three databases were considered as predicted target genes. Experimentally validated target genes were identified using miRWalk 3.0.

The biological functions of the key miRNAs were determined by GO clustering, including the three domains 'biological process' (BP), 'cellular component' (CC) and 'molecular function' (MF), and KEGG pathway analysis of the target genes using the Database for Annotation, Visualization and Integrated Discovery Bioinformatics Resources 6.7 (david-d.ncifcrf.gov/); corresponding results were displayed using the ggplot2 and gcookbook packages in R software.

Regulatory networks of the miRNA-target genes and miRNA-target gene signaling pathways were constructed and visualized using Cytoscape 3.4.0 (www.cytoscape.org/) in the Java Environment (16).

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

RT was performed using a PrimeScript RT Reagent kit (Takara Biotechnology Co., Ltd., Dalian, China). Stem-loop RT primers and PCR primers for each miRNA were purchased from Ribobio Co., Ltd. (cat. nos. U6 MQP-0202, rno-mir-125b-2-3p miRQ0004731-1-1, rno-mir-187-5p miRQ0017144-1-1, rno-mir-3542 miRQ0017796-1-1, rno-mir-92a-2-5p miRQ0017108-1-1 and rno-mir-99a-5p miRQ0000820-1-1). Fluorescence qPCR was performed on an ABI 7500 detection system (Applied Biosystems; Thermo Fisher Scientific, Inc., Waltham, MA, USA). Each 20 μl reaction mixture contained 2 μl template cDNA, 0.4 μM of each gene-specific primer, 10 μl SYBR Premix Ex Taq II (Tli RNaseH Plus, 2X), 0.4 μl ROX Reference Dye II (50X), and 6 μl sterilized Rnase-free water. The reactions were incubated at 50°C for 2 min, 95°C for 2 min, followed by 40 cycles of 95°C for 15 sec, and 60°C for 60 sec. The reactions were incubated at 50°C for 2 min, 95°C for 2 min, followed by 40 cycles of 95°C for 15 sec, and 60°C for 60 sec. With stable expression in the normal and ARM groups, U6 was selected as an endogenous reference gene for qPCR normalization purposes and the relative expression level of each miRNA was calculated using the 2−ΔΔCq method (17). RT-qPCR analysis was performed with three replicates.

Statistical analysis

Statistical analysis was performed using SPSS 21.0 software (IBM SPSS, Armonk, NY, USA). Statistically significant differences were identified using Student's t-test (two-tailed). All values are presented as the mean ± standard deviation. P<0.05 was considered to indicate a statistically significant difference.

Results

Screening of differentially expressed miRNAs and clustering analysis

The miRNAs potentially involved in ARM development were investigated using microarray screening of differentially expressed miRNAs at key time-points during anus formation (GD14, GD15 and GD16). A heatmap of the 761 miRNAs at the three time-points is shown in Fig. 1A. The results of the microarray screening revealed different expression patterns between the normal and ARM groups, and the majority of miRNAs exhibited increased expression in the ARM group. In the normal group, a higher expression level of miRNAs was detected at GD15 compared with GD14 and GD16, and the expression trend was opposite to that in the ARM group.

Figure 1

miRNA expression and miRNA-target gene regulatory network. In the heatmap, green and red indicate a decrease and increase of miRNA expression, respectively. The intensity of the color corresponds to the relative signal level on a logarithmic scale. (A) Hierarchical heatmap of the overall microarray screening. (B) Venn diagram of differentially expressed miRNAs at three time-points. (C) Hierarchical heatmap of differentially expressed miRNAs.

A total of 38 miRNAs were upregulated on GD14 in the ARM group compared with the normal group; 32 were upregulated and 17 were downregulated at GD15 in the ARM group compared with the normal group; and 42 were upregulated at GD16 in the ARM group compared with the normal group. The overlap of these miRNAs between the three time-points is shown in Fig. 1B. The detailed information and expression patterns of the differentially expressed miRNAs at the three time-points are shown in Table I and Fig. 1C.

Table I

Differentially expressed miRNAs at GD14, GD15 and GD16.

miRNA	lgFC	FDR
G14; N vs. A (38)
miR-125b-2-3p	5.6523	0.0181
miR-92a-2-5p	4.9863	0.0033
miR-99a-5p	4.7777	0.0207
miR-331-3p	4.7257	0.0171
miR-410-3p	4.5907	0.0124
miR-26b-3p	4.2617	0.0215
miR-190a-3p	4.2564	0.0052
miR-212-5p	4.1900	0.0495
miR-802-5p	4.1618	0.0230
miR-3586-3p	4.1185	0.0394
miR-143-5p	4.1014	0.0133
miR-132-3p	4.0487	0.0068
miR-200a-3p	3.9410	0.0020
miR-6319	3.8708	0.0209
miR-497-3p	3.8449	0.0105
miR-29a-3p	3.7200	0.0008
miR-133a-5p	3.6877	0.0338
miR-501-5p	3.4335	0.0379
miR-369-5p	3.4053	0.0480
miR-322-5p	3.3726	0.0006
miR-146a-3p	3.1359	0.0068
miR-9b-3p	3.0811	0.0077
miR-31b	3.0555	0.0216
miR-301b-3p	3.0376	0.0278
let-7i-3p	3.0366	0.0087
miR-291b	3.0339	0.0077
miR-490-3p	2.9198	0.0040
miR-30e-3p	2.8834	0.0112
miR-421-3p	2.8269	0.0162
miR-6329	2.6684	0.0479
miR-411-5p	2.4638	0.0261
miR-344i	2.2869	0.0288
miR-326-3p	1.8954	0.0491
miR-376b-5p	1.8833	0.0128
miR-666-3p	1.8752	0.0217
miR-192-3p	1.7802	0.0307
miR-3557-5p	1.6589	0.0399
miR-376b-3p	1.0623	0.0007
G15; N vs. A (49)
miR-741-3p	4.9521	0.0013
miR-935	4.9112	0.0077
miR-542-5p	3.3498	0.0025
miR-1249	3.2639	0.0007
miR-5132-3p	2.1773	0.0133
miR-129-2-3p	2.0751	0.0011
miR-409b	1.9867	0.0165
miR-7a-1-3p	1.8004	0.0432
miR-3550	1.7798	0.0096
miR-349	1.6717	0.0172
miR-874-5p	1.6689	0.0232
miR-3547	1.6623	0.0007
miR-551b-5p	1.6278	0.0027
miR-361-3p	1.5764	0.0205
miR-3546	1.5650	0.0199
miR-455-3p	1.5239	0.0484
miR-212-3p	1.5230	0.0169
miR-743a-3p	1.4840	0.0116
miR-196a-3p	1.3656	0.0000
miR-28-5p	1.3441	0.0328
miR-194-3p	1.3365	0.0083
miR-15b-5p	1.3145	0.0170
miR-191b	1.2964	0.0156
miR-540-3p	1.2545	0.0328
let-7b-3p	1.2488	0.0363
miR-191a-5p	1.2430	0.0180
miR-6315	1.1330	0.0307
miR-125b-5p	1.1254	0.0038
miR-125a-5p	1.1246	0.0073
miR-210-3p	1.1245	0.0308
miR-652-5p	1.0856	0.0302
miR-20b-5p	1.0384	0.0139
miR-144-5p	−6.6610	0.0130
miR-539-3p	−5.9386	0.0043
miR-708-5p	−5.4602	0.0468
miR-412-5p	−5.1433	0.0053
miR-187-5p	−4.9449	0.0021
miR-3542	−4.7610	0.0261
miR-224-3p	−4.2339	0.0112
miR-101a-5p	−4.2198	0.0133
miR-200c-5p	−4.0638	0.0208
miR-19b-2-5p	−3.4133	0.0476
miR-384-5p	−3.2683	0.0244
miR-135a-5p	−2.8032	0.0037
miR-204-3p	−2.7952	0.0009
miR-741-5p	−2.7242	0.0060
let-7f-5p	−2.2251	0.0124
miR-702-5p	−1.4132	0.0039
miR-449c-3p	−1.2088	0.0303
G16; N vs. A (42)
miR-409a-5p	6.3869	0.0271
miR-92a-2-5p	4.8504	0.0128
miR-484	4.7546	0.0084
miR-673-3p	4.5965	0.0110
miR-125b-2-3p	4.3396	0.0171
miR-3542	4.2864	0.0007
miR-127-5p	4.2316	0.0246
miR-6325	4.1901	0.0126
miR-23b-5p	4.1463	0.0001
miR-3068-3p	3.9677	0.0370
miR-99a-5pa	3.8542	0.0053
miR-143-5p	3.7453	0.0142
miR-6318	3.7171	0.0420
miR-9a-5p	3.6453	0.0070
miR-871-3p	3.5678	0.0002
miR-29b-1-5p	3.5150	0.0137
miR-219a-5p	3.4518	0.0295
miR-1912-5p	3.2806	0.0171
miR-187-5pa	3.2137	0.0074
miR-9a-3p	3.1949	0.0292
miR-296-5p	2.8674	0.0025
miR-483-3p	2.6111	0.0454
miR-500-3p	2.5249	0.0157
miR-881-3p	2.4781	0.0438
miR-665	2.4316	0.0108
miR-3552	2.3654	0.0418
miR-351-5p	2.1768	0.0025
miR-742-5p	2.0011	0.0080
miR-708-3p	1.9967	0.0259
miR-324-5p	1.9943	0.0130
miR-3573-3p	1.9430	0.0298
miR-29b-5p	1.8790	0.0485
miR-125a-5p	1.8771	0.0059
miR-485-5p	1.7577	0.0154
miR-129-2-3p	1.6690	0.0001
miR-191a-5p	1.6207	0.0372
miR-204-3p	1.6186	0.0182
let-7d-5p	1.6180	<0.0001
miR-3084b-5p	1.5984	0.0362
miR-125b-5p	1.2629	0.0247
miR-702-5p	1.2500	0.0419
miR-30c-1-3p	1.1605	0.0216

miRNAs marked with dashed lines are those with |log2(FC)| >4.25, and were selected for further analysis.

miR-99a-5p and miR-187-5p on GD16 had |log2(FC)| >4.25 on GD14 or GD15. miR/miRNA, microRNA; GD, gestational day; N, normal; A, anorectal malformation; lgFC, log2; FC, fold change.

Predicted and validated target genes of miRNAs

There were 18 differentially expressed miRNAs with |log2(FC)| >4.25 (miR-125b-2-3p, miR-144-5p, miR-187-5p, miR-190a-3p, miR-26b-3p, miR-331-3p, miR-3542, miR-409a-5p, miR-410-3p, miR-412-5p, miR-484, miR-539-3p, miR-673-3p, miR-708-5p, miR-741-3p, miR-92a-2-5p, miR-935 and miR-99a-5p), which were selected for further analysis. The numbers of predicted target genes for the above miRNAs were 229, 27, 104, 81, 177, 54, 219, 19, 57, 1, 112, 66, 17, 43, 143, 642, 148 and 5, respectively. Bone morphogenetic protein receptor type 2 (BMPR2) was the only experimentally validated target gene of miR-99a-5p in rats (18). The predicted and experimentally validated target genes of the 18 miRNAs were subjected to further analysis. The regulatory network between the miRNAs and their target genes is presented in Fig. 1D.

GO and KEGG enrichment analyses, and construction of the regulatory network

According to GO analysis of the target genes of the 18 miRNAs, the functions of DNA binding and gene regulation, metabolic and biosynthetic process of biomolecule, differentiation and projection of neuron were enriched. The top 15 enriched terms in the BP, CC and MF terms are presented in Fig. 2A.

Figure 2

GO and KEGG enrichment analyses ranked by enriched factor, which was calculated by -log(q-value). The q-value is a modified Fisher exact P-value provided by Database for Annotation, Visualization and Integrated Discovery enrichment analysis. (A) GO enrichment results; the top 15 GO categories are listed in the diagram. (B) Top 30 enriched signaling pathways from KEGG analysis. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function.

In terms of the KEGG analysis, a total of 49 signaling pathways were enriched, and the top 30 are listed in Fig. 2B. Among those pathways, 'axon guidance' was enriched the most, followed by 'phosphatidylinositol signaling system'. Notably, the enriched pathways included the Wingless Type MMTV integration site family (Wnt), mitogen-activated protein kinase (MAPK), and transforming growth factor-β (TGF-β) pathways, which are known to be closely associated with the process of embryonic development and ARM (6,19-21). Therefore, miRNA target gene-pathway networks were constructed for the miRNAs and target genes contributing to these three signaling pathways (Fig. 3A). A total of 15 core miRNAs associated with the signaling pathways of interest were screened out: miR-92a-2-5p, miR-125b-2-3p, miR-935, miR-190a-3p, miR-741-3p, miR-26b-3p and miR-539-3p were potentially associated with all three pathways; miR-484, miR-3542, miR-708-5p and miR-187-5p were potentially associated with two pathways; and miR-410-3p, miR-331-3p, miR-673-3p and miR-99a-5p were potentially associated with only one pathway. As demonstrated by the miRNA target gene-signaling pathway network, miR-92a-2-5p exhibited the most extensive connection with target genes in the relevant pathways.

Figure 3

miRNA-target gene-signaling pathway regulatory network and PCR validation. (A) miRNA-target gene-signaling pathway regulatory network. Diamonds and nodes represent miRNAs and target genes, among which the orange diamond and blue nodes represent an miRNA-mRNA pair validated by a dual luciferase experiment. Red polygons represent signaling pathways. Grey lines represent potential regulatory associations between miRNAs and their target genes. (B) Venn diagram of miRNAs differentially expressed at two time-points and included in the miRNA-target gene-signaling pathway regulatory network. Overlap between the two circles represent five miRNAs included in the network that are differentially expressed at two time-points and validated by PCR. (C) Relative expression levels of selected miRNAs in N and A groups. Relative miRNA expression is presented as the mean ± standard deviation. *P<0.05; **P<0.01. miR/miRNA, microRNA; PCR, polymerase chain reaction; N, normal; A, anorectal malformations; GD, gestational day.

Validation of important miRNAs

Among the 15 core differentially expressed miRNAs, five (miR-125b-2-3p, miR-92a-2-5p, miR-99a-5p, miR-187-5p and miR-3542) were differentially expressed at two time-points, and were selected for RT-qPCR validation (Fig. 3B). For all miRNAs, the relative expression levels were consistent with the microarray results (Fig. 3C). There were significant differences between the normal and ARM group for three miRNAs (miR-125-2-3p, miR-92a-2-5p and miR-99a-5p).

Discussion

The role of miRNAs in embryonic development and congenital deformities of gastrointestinal system have been described, however, the association between miRNAs and ARMs remains to be fully elucidated (8,10,11,22-24). Examining miRNA expression patterns may assist in identifying the pathogenesis of ARM. In the present study, miRNAs that may be involved in the course of ARM were profiled and bioinformatics analyses were performed to reveal the potential mechanism.

According to previous morphological investigations, rectourethal fistula and a common cloaca were present in all ARM rat fetuses, as a result of failed fusion of the urorectal septum and cloacal membrane at GD15 (4,13,25). Therefore, the present study detected aberrantly expressed miRNAs in ARM fetuses compared with normal rat fetuses at key time-points for anus formation (GD14, GD15 and GD16). The number of differentially expressed miRNAs at each time-point was 38, 49 and 42, respectively. Hierarchical analysis revealed that the increases of miRNAs in the ARM group were predominant, suggesting that increases in miRNAs may be one of the factors causing ARMs. In the normal group, there were differences in miRNA expression between the three time-points, which may be critical for embryonic development and organogenesis.

From the microarray results, 18 miRNAs with the highest significant FC values were selected for further bioinformatics analysis. By predicting target genes and searching for experimentally validated targets, an miRNA-target gene regulatory network was constructed. In a previous study by Jin et al (9), aberrantly expressed miRNAs in the terminal hindgut following the formation of ARM in Sprague-Dawley rat fetuses were profiled by the miRNA microarray, and miR-193 was screened out with |FC| >2. The present study focused on differentially expressed miRNAs prior to ARM having fully developed (GD14 and GD15), and the inclusion criteria of miRNAs selected for further bioinformatics analysis was |log2(FC)| >4.25; therefore, none of the miR-193 family members, including miR-193a-3p, miR-193a-5p, miR-193b-3p and miR-193b-5p, were included in the miRNAs selected for further analysis. Differences in the chip array and the samples used may also be factors in the discrepancy between these two studies.

Further functional analysis indicated that the target genes of the differentially expressed miRNAs were enriched in processes, including gene regulation and neuron development, and KEGG pathway analysis indicated that 49 significant signaling pathways were enriched, including the Wnt, MAPK and TGF-β signaling pathways, which have also previously been reported to be closely associated with ARM (6,19-21).

The neuron-associated terms that were enriched in the GO and KEGG functional analyses indicated that these miRNAs, were not only potentially associated with the development of ARMs, but may also be involved in neurode-velopmental abnormalities, which often appear concurrently with ARMs (26,27). In addition to the Wnt, MAPK and TGF-β pathways, the BMP signaling pathway is also important in ARMs (28-30). Wnt/β-catenin signal transduction dysregulation can lead to an ARMs phenotype (31). Studies have revealed that this may be caused by abnormal BMP signaling, as BMP4 and BMP7 were abnormally increased in the β-catenin-induced ARM phenotype, and furthermore, knockout of the BMP receptor, BMPR1A, partially rescued the malformation (32,33). However, in the previous study, only the BMRP1A receptor was knocked out, whereas the other BMP receptor, BMPR2, was not altered. Therefore, whether or not the knockout of BMPR2 and BMPR1A together is able to rescue the deformity to a greater degree, remains to be elucidated.

In the present study, the potential regulatory roles of the miRNAs and their corresponding target genes included in the enriched Wnt, MAPK and TGF-β signaling pathways were examined by constructing miRNA-target gene-pathway regulatory networks, which identified 15 core miRNAs indirectly associated with the these three signaling pathways. Among the 15 core miRNAs, five miRNAs (miR-125b-2-3p, miR-92a-2-5p, miR-99a-5p, miR-187-5p and miR-3542) were selected for RT-qPCR validation, as they were differentially expressed at two time-points in the microarray. The results revealed that the relative expression levels of miR-125b-2-3p, miR-92a-2-5p and miR-99a-5p were consistent with the microarray, with statistical significance.

miR-92 is a member of the miR-17-92 cluster, which can regulate fetal development at the early stages, and may promote proliferation and inhibit differentiation during embryogenesis (11,34). miR-125, which is a highly conserved miRNA family, can regulate early embryonic development by affecting cell fate and differentiation (35-37). Previously, GO analysis demonstrated that miR-125 controls axon guidance pathways, synaptic plasticity and catabolic processes (38). This GO and KEGG analysis was similar to the functional enrichment results in the present study, which suggests that miR-125 may be crucial in the pathogenesis of ARM. Regarding the miR-99a/100-125b cluster, different members can exert homogenous functions. RNAs of the miR-125 family can directly bind to the 3′-UTR of p53, B cell lymphoma-2 (Bcl-2), Bcl-2-like 12 and Mcl-1, and regulate cell proliferation and apoptosis (39,40). Similarly, miR-99 family members can reduce the expression of homeobox A1, which downregulates the expression of the downstream Bcl-2 gene, and leads to reduced cell survival and enhanced apoptosis (41). In hematopoietic stem and progenitor cells, miR-99a/100-125b can regulate homeostasis by shifting the balance between TGF-β and Wnt signaling (42); and in early chondrogenic differentiation in rats, miR-99a and miR-125b also have critical regulatory roles, with miR-99a being particularly important and negatively regulating chon-drogenic differentiation by directly targeting BMPR2 (18). These previous studies revealed that the dysregulation of apoptosis-associated genes, BMP family genes and homeobox genes are associated with ARM; however, the functions exerted by miRNAs in the disease remain to be fully elucidated, and further investigation is required.

In conclusion, the dynamic changes in differentially expressed miRNAs were investigated in normal and ARM rat fetuses at key time-points during anorectal development (GD14-16). The aberrant expression of miR-125b-2-3p, miR-92a-2-5p and miR-99a-5p during this period suggested that these miRNAs may be involved in ARM in rats. However, the result of this investigation into ARMs are preliminary, and further investigations are required to investigate the potential target genes which may contribute to ARM pathogenesis. In addition, normal and ARM fetuses were compared to screen miRNAs that may function in the course of this disease, however, the study did not evaluate miRNA expression in fetuses without the ARM phenotype following ETU treatment, which may assist in revealing the protective mechanisms that occur to avoid ARM following ETU stimulation. Further in vivo and in vitro investigations are also required to clarify the specific mechanisms responsible for the development of ARMs.