Rapid Communication: MicroRNA co-expression network reveals apoptosis in the reproductive tract during molting in laying hens.

The aim of this study was to determine the regulatory mechanisms of molting and recrudescence via studying the micro-RNA (miRNA) expression in the oviduct of laying hens. We performed a cDNA microarray analysis in the magnum tissue from the oviduct to identify the whole miRNA profiles through the molting and recrudescence periods. A total of 35 laying hens (47-wk-old) were divided into 7 groups (0 d: a control group; 6 and 12 d: 2 molting-period groups fed on a high-zinc diet; and 20, 25, 30, and 35 d: 4 recrudescence-period groups fed on a normal diet after a 12-d period on a high-zinc diet). An miRNA co-expression network (miRCN) was generated using the differentially expressed miRNA (DEM) according to the entire data integration. The significantly co-expressed miRNA ( = 111) were highly differentially expressed from 12 to 20 d, which was a transition period between molting and recrudescence, while their expression patterns were contrary to the estrogen changes. The targets of highly connected miRNA ( = 12) indicated the significant biological pathways and gene ontology (GO) terms, such as MAPK and Wnt signaling and magnesium-ion binding, which are associated with apoptotic activities. These results suggest that the miRNA of the miRCN might play a role in the apoptotic progression of the reproductive tract during molting.

INTRODUCTION

Molting in female fowls is generally defined as the annual periodic phenomenon of feather shedding and replacement (Berry, 2003). In natural conditions, fowls undergo physiological changes at the end of each egg-laying cycle every year and then recover in preparation for a new egg-laying cycle. Otherwise, molting is artificially induced by feed withdrawal or modified diet to improve egg production rate and quality in commercial laying hens (al-Batshan et al., 1994; Webster, 2003). A molting process can be briefly divided into an induced molting period and a recrudescence period. The transition periods include dynamic morphologic and physiologic changes in the regression and regrowth of reproductive tissues following the dynamic fluctuations in hormones such as testosterone, progesterone, estradiol, and corticosterone (Jeong et al., 2013). Thus, we assumed that the avian female reproductive tissues exhibit molting-specific mechanisms in response to hormones.

Micro-RNA (miRNA) are short noncoding regulatory RNA that regulate the post-transcriptional expression of target mRNA (Vasudevan et al., 2007) and play important roles in immune surveillance and reproduction (Carletti and Christenson, 2009; Mehta and Baltimore, 2016). Like other organ systems (Harris et al., 2006; O'Rourke et al., 2007), miRNA are essential for the proper development and function of the female reproductive tract (Otsuka et al., 2008). Numerous miRNA have been implicated in the regulation of estrogens and their receptors by repressing the expression of their coregulators (Klinge, 2009). Therefore, we speculated that specific miRNA might be expressed in the avian reproductive tract during molting, even though a little is known about the network-based correlations among the miRNA expression levels during the artificially induced molting in laying hens.

Here, we determined the molting mechanism in laying hens by performing a cDNA microarray analysis to identify the whole miRNA profiles in their magnum part of the oviduct during molting and recrudescence and provided an miRNA co-expression network (miRCN) according to the whole data integration. Further, we report the miRNA that play important roles in regulating the progression of molting in laying hens.

MATERIALS AND METHODS

Molting and Recrudescence Induction

The experimental use of chickens for this study was approved by the Animal Care and Use Committee of Korea University. A total of 35 commercial White Leghorn laying hens aging 47 wk were provided ad libitum access to feed and water and subjected to standard poultry husbandry guidelines. Molting was artificially induced by adding 20 g/kg zinc to their diet as described previously (Jeong et al., 2013). Briefly, all hens were divided into 7 experimental groups. Except for the control group (normal feeding, n = 5), all hens were fed a high-zinc diet for 6 or 12 d; they completely ceased egg production within 12 d after feeding on the high-zinc diet. Recrudescence was induced in the groups that were treated for more than 12 d (20, 25, 30, and 35 d) by feeding a normal commercial diet for 8, 13, 18, and 23 d after a 12-d period of high-zinc diet.

Tissue Collection and Microarray Analysis

After euthanizing all hens by using 60%–70% carbon dioxide, we took the magnum tissue in the oviduct from each hen on each assigned day. A portion of the magnum from each oviduct was cut into 10- to 15-mm pieces, frozen immediately in liquid nitrogen, and stored at −80°C. For microarray analysis, the total RNA was extracted from the magnum from each hen at 7 different time points (0, 6, 12, 20, 25, 30, and 35 d after the onset of zinc feeding) and purified using an RNeasy Mini Kit (Qiagen, Valencia, CA). Subsequently, a microarray analysis was performed using the Affymetrix GeneChip Chicken Genome Arrays (Affymetrix, Santa Clara, CA). All experiments were performed using the individually independent RNA pools and 3 independent microarray chips. The signal intensity of each spot was calculated using the dChip software. Thereafter, we took all miRNA (n = 542) from the signal spots and determined the differential expression values of all miRNA (DEM) by calculating the net intensity ratios between the 6 treatment groups (6, 12, 20, 25, 30, and 35 d) and the control group (0 d). The DEM were weighted into the log2-based fold change (log2 FC) values.

miRCN Analysis and Visualization

We performed the miRCN analysis using the partial correlation coefficient with information theory (PCIT) algorithm (Reverter and Chan, 2008). After removing the miRNA that did not have a significant false discovery rate (FDR; q < 0.05) in any of the 6 groups, only those significant miRNA (n = 159) that had a stringent significant level in DEM (absolute log2 FC ≥ 1.0) on at least 1 of the 6 treatments days were employed to enhance the efficiency of network construction in the miRCN analysis. The PCIT algorithm was used to establish significant connections and to construct the network. A connection between 2 miRNA according to the significant correlations establishes an edge between 2 nodes in the construction of the miRCN. Then, the miRCN was visualized using the Cytoscape software, version 3.4.1 (Shannon et al., 2003). In the topological view of the network, the miRNA (nodes) that were closer together had more common neighbors than the others. The hierarchical heat map and expression patterns were analyzed on the basis of the log2 FC levels of the DEM detected in the 6 phases of the molting treatment. We used the customized code in the R package and the k-means clustering algorithm in the Multi Experiment Viewer software by narrowing down the optimal number of clusters with 1,000 iterations (Howe et al., 2011).

In Silico Functional Validation

The nodes (n = 12) with more than 30 connection degrees were further used to rebuild the core network using the CircleView. The 12 miRNA in the core-network were applied to investigate their target genes in the miRBase database (Griffiths-Jones et al., 2006). The genes that showed a target score value of greater than 90, as obtained from the 12 miRNA, were analyzed for the functional enrichment of pathways using the DAVID Bioinformatics Resources 6.7 (Huang et al., 2007). The significantly enriched pathways were represented by fold enrichment and a minus log10 P-value. The gene ontology (GO) annotations, including biological processes and molecular functions, were also enriched using the GO FAT option, which is a specific GO category filter in DAVID. Subsequently, the REVIGO visualization tool was applied to construct treemaps for the enriched GO terms (Supek et al., 2011).

RESULTS AND DISCUSSION

A previous study reported that the molting and recrudescence processes were successfully induced by the high-zinc feeding treatment (Jeong et al., 2013). The microarray analysis in the magnum tissue from the oviduct provided total miRNA profiles across the 6 time points in the molting period (Supplementary data 1). Subsequently, the significant DEM were identified in the network analysis. We created an miRCN with 111 nodes and 842 edges by selecting the significantly correlated miRNA with absolute co-expression correlations greater than 0.99 among those determined to be significant via the PCIT algorithm (Fig. 1A). The constructed miRCN was elaborately connected by 1 cluster. Among the 6 time points (6, 12, 20, 25, 30, and 35 d) in the molting period, most DEM showed a maximum level in 12 d (n = 94). Other nodes were observed in 20 d (n = 16) and 6 d (n = 1, miR-1607). The DEM were obviously clustered according to the molting progression. The hierarchical heat map showed a clearly separated expression pattern in 12- and 20-d clusters and in a 35-d cluster (Fig. 1B). Their density plots for the counted miRNA also represented the obviously clustered pattern. The relative DEM were revealed by a dynamic wave pattern with a peak visible in 12 to 20 d, which is a transition period between molting and recrudescence, followed by a dramatic decrease in the process of recrudescence (Fig. 1C).

Figure 1.

miRNA co-expression network (miRCN) through the molting period in laying hens. (A) A clustered miRCN of the significantly correlated miRNA across the 6 molting points. Node color and size represent the day of maximum differential expression of miRNA (d 6, orange; d 12, blue-green; and d 20, yellow) and log twofold-change level of the day. (B) Hierarchical heat map of miRCN and the color scale with density plot. (C) Differentially expressed pattern of the miRCN across the 6 molting points.

To further narrow down to more essential miRNA and to identify their biological significance, the highly connected nodes (> 30 edges, n = 12) in the core of miRCN were investigated by in silico functional analysis. The 12 miRNA were densely connected with each other (Fig. 2A). The target genes of these 12 miRNA (n = 745) were estimated from the miRBase database (Supplementary data 2). The genes were functionally enriched in the KEGG pathway terms (Fig. 2B). The top 3 significantly enriched pathway terms were the MAPK signaling pathway, Wnt signaling pathway, and ubiquitin-mediated proteolysis. The MAPK family members are crucial for the apoptotic regulation in cells (Wada and Penninger, 2004). In addition, the MAPK and Wnt signaling pathways have been reported as the core mechanisms of cell death in the development of vertebrates (Guardavaccaro and Clevers, 2012; Arya and White, 2015). The ubiquitin-mediated proteolysis is also naturally observed in the regulation of apoptosis (Wójcik, 2002). The GO term enrichment in the biological process showed a significant term for the negative regulation of gene expression; this term represents the interactions between miRNA and targets (Fig. 2C). Moreover, the magnesium-ion binding, the first significantly enriched GO term in a molecular function (Fig. 2D), is essential to regulate numerous cellular functions, including signaling pathways and metabolic cycles, and induces apoptotic cell death (Dribben et al., 2010; Romani, 2011).

Figure 2.

Core-clustered miRNA in the miRCN and functional term enrichment via in silico target gene identification. (A) Reconstruction of core-clustered miRNA in the miRCN. (B) Enriched KEGG pathways of the target genes for core-clustered miRNA. Gene ontology (GO) treemaps with regard to the biological process (C) and molecular function (D) terms of the target genes for core-clustered miRNA.

The molting of laying hens accompanies dynamic changes in the morphology, physiology, and function of the reproductive tract (Jeong et al., 2013). Particularly, the oviductal regression is achieved through apoptotic processes during the molting period in chicken (Heryanto et al., 1997). In the present study, the miRCN and their target genes also revealed significant biological pathway and GO terms that are associated with apoptotic activities. Moreover, our results showed that the DEM pattern dramatically progressed in the molting period and gradually regressed in the recrudescence period. The DEM pattern was antagonistic, resembling the changes in the levels of plasma estradiol (E2) in laying hens during a forced molting (Hoshino et al., 1988). Estrogen plays a major role in the regulation of the reproductive system as well as in the protection from apoptosis (Monroe et al., 2002; Cui et al., 2013). A previous report demonstrated the apoptotic processes induced by the withdrawal of estrogen regresses the chick oviduct (Berger and Sanders, 2000). Moreover, a decrease in estrogen is involved in the regulation of apoptotic processes and leads to oviduct regression during induced molting in hens (Monroe et al., 2002). As Klinge (2009, 2012) reviewed the role of the estrogen-regulated miRNA expression, E2 regulates selective miRNA expressions in female mammalian reproductive tissues. A recent study reported the E2-mediated regulation of miRNA in the chicken oviduct (In Lee et al., 2015). Therefore, the miRNA of the miRCN created in this study can be suggested to be the estrogen-mediated regulation factors in the induction of molting. Moreover, we supposed that the regression or deficiency of miRNA expression via estrogen restoration is related to the induction of recrudescence from molting. Further studies are required to functionally validate the interactions and regulations of miRNA mediated by estrogen secretions during the molting period.