Identification of microRNA-Associated-ceRNA Networks Regulating Crop Milk Production in Pigeon (Columba livia)

Pigeon belongs to altrices. Squab cannot forage independently. Nutrition can only be obtained from crop milk secreted by male and female pigeon. miRNA could regulate many biological events. However, the roles of miRNA and ceRNA in regulating crop milk production are still unknown. In this study, we investigated the miRNAs expression profile of female pigeon crop, explored the potential key genes, and found the regulatory mechanisms of crop milk production. A total of 71 miRNAs were identified differentially expressed significantly. Meanwhile, miR-20b-5p, miR-146b-5p, miR-21-5p, and miR-26b-5p were found to be the key miRNAs regulating lactation. Target genes of these miRNAs participated mainly in cell development; protein and lipid synthesis; and ion signaling processes, such as cell-cell adhesion, epithelial cell morphogenesis, calcium signaling pathway, protein digestion, and absorption. In the ceRNA network, miR-193-5p was located in the central position, and miR-193-5p/CREBRF/LOC110355588, miR-460b-5p/GRHL2/MSTRG.132954, and miR-193-5p/PIK3CD/LOC110355588 regulatory axes were believed to affect lactation. Collectively, our findings enriched the miRNA expression profile of pigeon and provided novel insights into the microRNA-associated-ceRNA networks regulating crop milk production in pigeon.


Introduction
Centuries ago, bird fanciers and naturalists had known that crops of parental pigeons could produce materials that nourished their springs [1]. Bernard (1859) first named it as crop milk and found that it consisted of clumps of epithelial cells shedding from the mucosal layer. As the only nutrition source of squab [2], crop milk has high nutritional value at the 1st, 2nd and 3rd weeks after hatching [3][4][5]. The lipid content of crop milk consisted mainly of triglycerides, along with phospholipids, cholesterol, free fatty acids, cholesterol esters, and diglycerides [6]. Young pigeons fed with crop milk showed remarkable growth rates, reaching 8-, 18-and 22-fold of the hatching weight. Feeding with granules including crop milk could improve the growth of domestic chickens and weaning rats significantly [7][8][9][10][11]. The crop has considerable expandability, particularly in granivores. Comparing with other poultry, except for temporary storage of ingesta and softening food, pigeon crop is more specialized for its milk production ability [12].
There were few studies on the lactating process of pigeon crop, and prolactin was the most widely investigated. As reported, prolactin promotes mitosis and early signal transduction, and growth of special epithelial cells lining the crop of pigeons could regulate the synthesis of crop milk proteins by activating the IRS1/Akt/TOR signaling pathway [13], which promotes the formation of crop milk [1,4]. Research on the molecular mechanism of

Ethics Statement
All experiments using animals were approved by the animal care and use committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS 2018-3). All procedures were conducted in accordance with institutional animal ethics guidelines set by the Ministry of Agriculture of the People's Republic of China.

Experimental Animals and Sample Collection
Ten healthy female White King pigeons were obtained from a pigeon breeding farm in Beijing, including five lactating (at lactating day 3) and five non-lactating birds. Pigeons were slaughtered by CO 2 asphyxiation. The crop tissues were removed aseptically. Then they were cut into small pieces of 1 cm thickness and frozen in liquid nitrogen immediately.

RNA Isolation and Quality Assessment
The total RNA was extracted from samples of the crop tissue using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The RNA purity was checked using a Nanodrop Spectrophotometer (Model kaiaoK5500 ® , Beijing, China). The integrity and concentration of the RNA were measured by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, San Diego, CA, USA).

miRNA Library Construction and Sequencing
Total RNA was separated by 15% agarose gels to extract the small RNA (18-30 nt). After being precipitated by ethanol and centrifugal enrichment of small RNA sample, the library was prepared according to the method and process of Small RNA Sample Preparation Kit (Illumina, RS-200-0048, San Diego, CA, USA). The main contents are as follows: 1. Connecting the 3 adaptor to the separated small RNA. 2. Connecting the 5 adaptor to the separated small RNA. 3. Quantitative real-time PCR. 4. Recycling strips of 145-160 bp (22)(23)(24)(25)(26)(27)(28)(29)(30). After library construction, miRNA sequencing was performed by using the Illumina HiSeq TM 2500 platform (Annoroad, Beijing, China) and generating 50 bp single-end reads.

Quality Control, Mapping and Acquisition of Sequencing Data
Raw Data were processed with perl scripts to ensure the quality of data used in the following analysis. The adopted filter criteria are: 1. Reads without 3 adapter were filtered out. 2. Reads with the number of bases whose phred quality value was no more than 19 accounting for more than 15% were regarded as low-quality reads and were filtered out. 3. Reads with the number of N bases accounting for more than 10% were filtered out.
The assembled transcripts were further subjected to a stringent filtering progress to obtain the candidate lncRNAs. First, we removed the transcripts with one exon and with lengths of less than 200 bp. Next, we calculated the reads' coverage of transcripts using String Tie [30] and removed transcripts with a reads' coverage of less than three. Subsequently, we eliminated known genes, pre-micro RNA and other non-coding RNAs (rRNA, tRNA, snoRNA and snRNA) using Cuffcompare.

Protein-Protein Interaction Analyze
To get insight into the potential relationships of mRNAs, differentially expressed between lactating and non-lactating pigeon crops, PPI network was constructed using STRING (version 11.0, https://string-db.org/cgi/input.pl) and visualized in the cytoscape (version 3.6.1).

Construction of ceRNA Network
The ceRNA network was visualized using cytoscape software. The topological features of this ceRNA network were calculated by a built-in network analyzer tool in cytoscape software.

Quantitative Real-Time PCR
The RNA samples were reverse transcribed into cDNA using the miRcute plus miRNA first-strand cDNA Kit (TIANGEN, Beijing, China) following the manufacturers instructions. The concentration of the RNA samples was adjusted to 200 ng/µL by adding RNase free ddH 2 O. The first step of the reverse transcription was conducted at 42 • C for 60 min. It contained 5 µL RNA samples, 10 µL 2× miRNA RT reaction buffer, 2 µL miRNA RT enzyme mix, and 3 µL RNase free dH 2 O. The second step was conducted at 95 • C for 3 min. qRT-PCR was performed using the miRcute plus miRNA qPCR Kit (TIANGEN, China) and the ABI Quant Studio 7 Flex Real-Time Detection System (Life Technologies Holdings Pte Ltd., USA). Each 10 µL qRT-PCR mixture contained 5 µL 2× miRcute plus miRNA premix (SYBR&ROX), 0.2 µL reverse primer (10 µm), 0.2 µL forward primer, 1µL cDNA (100 ng), and 3.6 µL RNase free ddH 2 O. After an initial denaturing at 95 • C for 15 min, there were 40 cycles of amplification (94 • C for 20 s and 60 • C for 34 s), followed by a thermal denaturing (95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s) to generate melting curves. The miRNA primer design software was provided by Vazyme Biotech Co., Ltd. Nanjing, China.
The RNA samples were reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) following the manufacturer's instructions. The concentration of the RNA samples was adjusted to 200 ng/µL by adding RNase free ddH 2 O. The first step of the reverse transcription was conducted at 42 • C for 2 min. It contained 5 µL RNA samples, 1 µL gDNA Eraser, 2 µL 5× gDNA Eraser Buffer, and 2 µL RNase free ddH 2 O. The second step was conducted at 37 • C for 15 min, then at 85 • C for 5s. It contained 10 µL last step reaction solution, 1 µL PrimeScript RT Enzyme Mix, 4 µL 5× PrimeScript Buffer, 1 µL RT Primer Mix, and 4 µL RNase free ddH 2 O. A qRT-PCR was performed using the PrimeScript One Step RT-PCR Kit (TaKaRa, Dalian, China) and the ABI Quant Studio 7 Flex Real-Time Detection System (Life Technologies Holdings Pte Ltd., Carlsbad, USA). Each 10 µL qRT-PCR mixture contained 5 µL SYBR Premix Ex TaqTM II, 0.5 µL (10 pM), and in each primer there was 0.2 µL ROX Reference Dye II (50×), 1.5 µL cDNA (100 ng) and 2.3 µL RNase free ddH 2 O. After an initial denaturing at 95 • C for 3 min, there were 40 cycles of amplification (95 • C for 30 s and 60 • C for 34 s), followed by a thermal denaturing (95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s) to generate melting curves. The primers of the genes were designed by the Primer Premier 5.
The primers in this study were summarized in Table S1. All mRNA expression data were normalized to hypoxanthine phosphoribosyltransferase 1 (HPRT1). All miRNA expression data were normalized to U6. The relative expression of the genes was calculated using the 2 −∆∆Ct method.

Data Accessibility
The sequencing data have been submitted to the SRA database (accession number: PRJNA612642). SRA records could be accessed via the following link (https://www.ncbi. nlm.nih.gov/sra/PRJNA612642).

Identification of Pigeon Small RNAs
As reported in Ma et al., the clean reads rate in individual samples was above 88%, and the clean Q30 bases rate was above 90% after data filtering [28].  Figure 1a), respectively. Known miRNA, rRNA, tRNA, snRNA, snoRNA, and novel miRNAs were identified after annotating for small RNAs. The number of known miRNAs accounted for more than 24% in all annotated small RNAs (Supplementary Table S2). In order of scale, miRNA > tRNA > rRNA > snoRNA in all the libraries sequenced. Others were less than 1%.     A total of 386 miRNAs (188 known miRNAs, 198 novel miRNAs) was found to be expressed in the crop of lactating and non-lactating pigeons (Supplementary Tables S3 and S4). A good correlation among individuals in the stage was observed in the miRNA cluster map and density distribution map (Supplementary Table S4, Figure 1b,c). A total of 50 known miRNAs and 21 novel miRNAs were differentially expressed (Supplementary Table S5, Figure 1d,e). Top 10 differential expressed miRNAs were shown in Table 2, including four up-regulated and six down-regulated ones.

miRNA Target Prediction and Analyze
mRNAs and lncRNAs targets of the differentially expressed miRNA were predicted using miRanda, PITA and Target scan databases. We identified 705 pairs of miRNA-mRNA (including 39 known miRNAs, 17 novel miRNAs and 620 mRNAs), 172 pairs of miRNA-lncRNA (including 28 known miRNAs, 12 novel miRNAs and 75 lncRNAs) and a pair of miRNA-circRNA. The characteristic descriptions of RNAs are shown in Figure 2.  GO and KEGG pathway analyses of the target mRNAs of differentially expressed miRNA were performed. The target genes were enriched in 52 GO biological process items and 14 KEGG pathways (Supplementary Tables S6 and S7), including development process, cell-cell adhesion, epithelial cell morphogenesis, calcium signaling pathway, pro- GO and KEGG pathway analyses of the target mRNAs of differentially expressed miRNA were performed. The target genes were enriched in 52 GO biological process items and 14 KEGG pathways (Supplementary Tables S6 and S7), including development process, cell-cell adhesion, epithelial cell morphogenesis, calcium signaling pathway, protein digestion and absorption, and glycosaminoglycan biosynthesis-of keratan sulfate ( Figure 3). GO and KEGG pathway analyses of the target mRNAs of differentially expressed miRNA were performed. The target genes were enriched in 52 GO biological process items and 14 KEGG pathways (Supplementary Tables S6 and S7), including development process, cell-cell adhesion, epithelial cell morphogenesis, calcium signaling pathway, protein digestion and absorption, and glycosaminoglycan biosynthesis-of keratan sulfate ( Figure 3).

PPI Network
The 620 target mRNAs of the differentially expressed miRNAs were submitted to string online tool where PPI network was constructed. Over 45% of the target mRNAs were nodes in the complex PPI network (Figure 4, Supplementary Table S8). Glutamylprolyl tRNA synthetase (EPRS), XP 005509099.1, sterol regulatory element-binding protein 1 (SREBF1), acetyl-CoA carboxylase α (ACACA), and protein phosphatase type 1 α (PPP1CA) were identified as the most highly connected nodes in the network.

PPI Network
The 620 target mRNAs of the differentially expressed miRNAs were submitted to string online tool where PPI network was constructed. Over 45% of the target mRNAs were nodes in the complex PPI network (Figure 4, Supplementary Table S8). Glutamyl-prolyl tRNA synthetase (EPRS), XP 005509099.1, sterol regulatory element-binding protein 1 (SREBF1), acetyl-CoA carboxylase α (ACACA), and protein phosphatase type 1 α (PPP1CA) were identified as the most highly connected nodes in the network.

miRNA-mRNA-lncRNA/circRNA Interaction Network
The interaction information of miRNA-mRNA, miRNA-lncRNA and miRNA-circRNA were predicted by PITA, miRanda and target scan. Based on these data, we constructed miRNA-mRNA-lncRNA/circRNA ceRNA regulatory networks and visualized them via cytoscape software. LncRNAs, circRNAs and mRNAs were connected by the same target miRNAs, and 1619 regulatory triplets (1024 pairs known miRNA-mRNA-lncRNA, 575 pairs novel miRNA-mRNA-lncRNA and 20 pairs miRNA-mRNA-circRNA) were obtained (Figure 5a-c). The topological features of ceRNA network were assessed by a built-in network analyzer tool in cytoscape software, including betweenness, network degree, and closeness centrality.
Statistics of four types of RNA degree value, betweenness centrality and closeness centrality in the network, retain higher degree value RNAs (Table 3) as the hub of ceRNA network. miR-193-5p, miR-92-2-5p, miR-7b-5p, and CREB3 regulatory factor (CREBRF) were at the core of the network, with the most lncRNAs and mRNAs associated with them. cli-let-7c-5p was the only miRNA connected to the lone target circRNA (circ_0003020). Genes 2021, 12, x FOR PEER REVIEW 8 of 19 Figure 4. Summary of protein-protein interaction network of target genes of the differentially expressed miRNAs between lactating and non-lactating pigeon crops. Different colors represent modules, and size of the node represents connectedness of the node.

miRNA-mRNA-lncRNA/circRNA Interaction Network
The interaction information of miRNA-mRNA, miRNA-lncRNA and miRNA-circRNA were predicted by PITA, miRanda and target scan. Based on these data, we constructed miRNA-mRNA-lncRNA/circRNA ceRNA regulatory networks and visualized them via cytoscape software. LncRNAs, circRNAs and mRNAs were connected by the same target miRNAs, and 1619 regulatory triplets (1024 pairs known miRNA-mRNA-lncRNA, 575 pairs novel miRNA-mRNA-lncRNA and 20 pairs miRNA-mRNA-circRNA) were obtained (Figure 5a-c). The topological features of ceRNA network were assessed by a built-in network analyzer tool in cytoscape software, including betweenness, network degree, and closeness centrality.
Statistics of four types of RNA degree value, betweenness centrality and closeness centrality in the network, retain higher degree value RNAs (Table 3) as the hub of ceRNA network. miR-193-5p, miR-92-2-5p, miR-7b-5p, and CREB3 regulatory factor (CREBRF) were at the core of the network, with the most lncRNAs and mRNAs associated with them. cli-let-7c-5p was the only miRNA connected to the lone target circRNA (circ_0003020).

Functional Annotation of the mRNAs in ceRNA Network
Functional enrichment analysis of mRNA targets of differentially expressed miRNAs was conducted to explore the function of mRNA in ceRNA network. These mRNAs were significantly enriched in 25 GO terms, including cell motility, cell migration, and calcium ion binding (Figure 6a). Only calcium signaling pathway and focal adhesion were significantly enriched in the KEGG pathway. The DEGs in the pathway are shown in the heat maps (Figure 6b,c).
As shown in Figure 7, these differentially expressed genes were involved in various steps of fatty acid elongation and synthesis. In particular, ACACA and elongation of very long-chain fatty acid proteins (ELOVL) have also been identified as differentially expressed. The synthesized products were also associated with signal pathways such as cutin, suberin and wax biosynthesis and fatty acid degradation.

Functional Annotation of the mRNAs in ceRNA Network
Functional enrichment analysis of mRNA targets of differentially expressed miRNAs was conducted to explore the function of mRNA in ceRNA network. These mRNAs were significantly enriched in 25 GO terms, including cell motility, cell migration, and calcium ion binding (Figure 6a). Only calcium signaling pathway and focal adhesion were significantly enriched in the KEGG pathway. The DEGs in the pathway are shown in the heat maps (Figure 6b,c).

Construct Fatty Acid Biosynthesis Pathway Function Network
Lipids were extremely important nutrients in crop milk, and we identified many genes that played an important role in the extension, synthesis and metabolism of fatty acids. So, we constructed a metabolic sub-network related to fatty acid biosynthesis (Figure 7). Metabolic networks were involved in insulin, PI3K-Akt, ECM-receptor interaction, GAP junction, fatty acid biosynthesis and extension, FoxO, and cell cycle signaling pathways. Purple ovals represent proteins involved in pathways. Pink ovals represent the differentially expressed genes identified in this study. The nine genes (hepatocyte growth factor (HGF), laminin subunit γ 2 (LAMC2), integrin subunit β 8 (ITGB8), platelet derived hydratase (trifunctional protein), and α subunit (HADHA) were involved in the ceRNA network.
As shown in Figure 7, these differentially expressed genes were involved in various steps of fatty acid elongation and synthesis. In particular, ACACA and elongation of very long-chain fatty acid proteins (ELOVL) have also been identified as differentially expressed. The synthesized products were also associated with signal pathways such as cutin, suberin and wax biosynthesis and fatty acid degradation.

The Exploration of ceRNA Regulates Axes
The key regulation axis of ceRNA regulating the crop milk production was screened by integrating the results of targeted relationship and functional regulation network. The results are shown in Table 4. miR-193-5p and many associated with pigeon milk production key genes were predicted as a target relationship. MSTRG.65211 and LOC110355588 might be in competition with miR-193-5p to regulate gene expression.

The Exploration of ceRNA Regulates Axes
The key regulation axis of ceRNA regulating the crop milk production was screened by integrating the results of targeted relationship and functional regulation network. The results are shown in Table 4. miR-193-5p and many associated with pigeon milk production key genes were predicted as a target relationship. MSTRG.65211 and LOC110355588 might be in competition with miR-193-5p to regulate gene expression.

Validation of Differentially Expressed RNAs
A total of seven miRNAs and seven mRNAs were chosen for qRT-PCR validation of the RNA sequencing results in lactating and non-lactating pigeon crop. ACACA, E-cadherin (CDH1), ELOVL6, DOT1 like histone lysine methyltransferase (DOT1L), epithelial splicing regulatory protein 2 (ESRP2), and miR-200a-5p were identified to be upregulated in lactating crop compared to non-lactating, whereas tetratricopeptide repeat domain 28 (TTC28), PIK3CD, let-7c-5p, let-7a-5p, miR-26-5p, miR-10a-5p, miR-99-5, and miR-125-5p were downregulated (Figure 8). These results were consistent with the sequencing results. Therefore, differential expression of these mRNAs might be involved in the process of crop milk production. Figure 8. An illustration of the qRT-PCR confirmation for the RNA-seq. The correlations of the gene expression level of the 14 differentially expressed genes between lactating and non-lactating crops using RNA-Seq and a qRT-PCR. The X-axis and Y-axis show the log2Fold change (T/C) measured by the RNA-seq and qRT-PCR, respectively. T represents lactating group and C represents non-lactating group. The hypoxanthine phosphoribosyltransferase 1 (HPRT1) gene was used as a housekeeping internal control.

Discussion
Male and female pigeons have the ability to produce nutrients in their crop for the nourishment of their young. The production of the nutrient has been likened to lactation in mammals, and hence, it has been called crop milk [39]; meanwhile, the pigeon crop undergoes significant changes to the tissue structure during lactation [17]. miRNAs have Figure 8. An illustration of the qRT-PCR confirmation for the RNA-seq. The correlations of the gene expression level of the 14 differentially expressed genes between lactating and non-lactating crops using RNA-Seq and a qRT-PCR. The X-axis and Y-axis show the log 2 Fold change (T/C) measured by the RNA-seq and qRT-PCR, respectively. T represents lactating group and C represents non-lactating group. The hypoxanthine phosphoribosyltransferase 1 (HPRT1) gene was used as a housekeeping internal control.

Discussion
Male and female pigeons have the ability to produce nutrients in their crop for the nourishment of their young. The production of the nutrient has been likened to lactation in mammals, and hence, it has been called crop milk [39]; meanwhile, the pigeon crop undergoes significant changes to the tissue structure during lactation [17]. miRNAs have demonstrated to play crucial functions in lipid metabolism, adipose tissue development and immune response [40][41][42]. In the present study, we found that over 40% of the small RNAs in pigeon crop were miRNAs; a total of 71 miRNAs were identified as significantly differentially expressed. We hypothesized that these differentially expressed miRNAs are involved in crop tissue morphological changes and pigeon milk synthesis. Histological studies showed that the layers of cells in the crop during lactation proliferated and that there were shed epithelial cells in the pigeon milk [43,44]. We noticed that many miRNAs involved in cell proliferation and migration have been identified, such as miR-20b-5p, which was proved to participate in the epithelial mesenchymal transition (EMT), migration and invasion process of prostate cancer [45]. MiR-146b-5p played a regulatory role in the proliferation and differentiation of human fibroblasts and chicken myoblasts [46,47]. MiR-21-5p is up-regulated in a variety of cancer cells such as non-small cell lung cancer, Hepatocellular carcinoma, melanoma, and colon adenocarcinoma cells, and has been shown to promote cell proliferation [48][49][50][51]. MiR-193-5p downregulation has significant angiogenic effect by inducing migration and proliferation in myocardial microvascular endothelial cells. Insulin-like growth factor 2 (IGF2) acted as a direct regulator to prevent this process from happening [52]. Gillespie suggested that the exfoliation of the crop epithelium due to inadequate blood supply leads to cellular keratinization [39], which may be an effect of miR-193-5p upregulation and inhibition of angiogenesis during lactation. Furthermore, Luo et al. found that miR-26b-5p serves as a direct negative regulator of TCF-4 expression within human adipose-derived mesenchymal stem cell, leading to inactivation of the Wnt/β-catenin pathway and thereby promoting the adipogenic differentiation of these cells in vitro [53]. This function is closely related to the characteristics of pigeon milk with high fat and nutrition. Based on the above evidence, we speculated that these miRNAs played an important role in the proliferation and exfoliation of the crop epithelial cells during the lactation period.
In order to further explore the functions of miRNAs, target gene prediction and functional enrichment analysis were performed. We identified 620 differentially expressed target genes and 75 differentially expressed target lncRNAs for differentially expressed miRNAs. By performing functional enrichment analysis of these genes, it was found that GO categories and pathways for DEGs such as development process, cell-cell adhesion, epithelial cell morphogenesis, calcium signaling pathway, protein digestion and absorption, and glycosaminoglycan biosynthesis of keratan sulfate signaling pathways were enriched between lactation and non-lactation periods. PPI analysis revealed EPRS, SREBF1 and ACACA closely associated with adipogenesis, degeneration, and lipid accumulation [54,55]. These results suggested that genes regulating substance synthesis remained active to support the synthesis of a large number of milk compositions in lactation.
Further, we found several GO terms related to cell adhesion and calcium signaling pathway. Classical cadherins are the key molecules that control cell-cell adhesion and are essential regulators of tissue homeostasis that govern cellular function and development, by transducing adhesive signals to a complex network of signaling effectors and transcriptional programs [56]. Specifically, CDH1 and M-cadherin (CDH15) are calciumdependent intercellular adhesion molecules. CDH1 played an important role in epithelial adhesion, geometrical packing and Spatial effort [57,58]. We suspect CDH15 has a similar function; Claudins as the cell adhesion protein are members of the membrane-associated guanylate kinases (MAGUK) family proteins and the main bridge proteins of connecting claudins with the underlying actin cytoskeleton [59]. Studies found that the formation of adhesion junctions could also promote tight junction assembly by altering the lipid composition of the plasma membrane. Low Ca 2+ led to internalization of claudins and reduced cholesterol in the plasma membrane [60]. This suggested that Ca 2+ concentration and the genes involved in cell adhesion might be related to the shedding of epithelial cells in the crop of lactating pigeon. Other genes involved in cell adhesion played a major role in the immune system and neural system, majorly participating in T cell receptor signaling pathway, inducing T cell activation and differentiation [61].
Regarding the constructed miRNA-associated-ceRNA networks, we desired to explore the roles of key RNA molecules in crop milk formation and epithelial cell morphogenesis. We found that many genes were identified as predicted target genes of miRNAs that might play critical roles in metabolism through several pathways. So HADHA, adrenoceptor β 3 (ADRB3), ACACA, ELOVL6, and 7-dehydrocholesterol reductase (DHCR7) were identified as key genes for lipid and fatty acid metabolism. The product of ACACA is the rate-limiting enzyme in the biosynthesis of long-chain fatty acids and was reported to be an important enzyme involved in the synthesis of saturated fatty acids. This enzyme was expressed ubiquitously, but the highest level was found in lipogenic tissues such as the liver, adipose tissue and lactating mammary gland [62][63][64]. Similarly, ADRB3 was a key gene in the process of adipose tissue browning and was reportedly down-regulated in adipocytes of obese people [65]. This allowed researchers to focus on adipose tissue browning in mice [66] and humans [67,68]. In this study, ADRB3 was downregulated in lactating pigeon crop, corresponding to the period of a high proportion of fat (33.8%) in crop milk [4]. Both HADHA and DHCR7 have been shown to be involved in lipid processes, including cardiolipin re-modeling, fatty acid β-Oxidation in human myocardium [69] and vitamin D-lipid interaction [70], respectively. Seven members of the ELOVL protein family have been described in mammals [71,72]. They are closely related to fatty acid metabolism and hepatic insulin sensitivity. ELOVL fatty acid elongase 4 (ELOVL4), ELOVL fatty acid elongase 5 (ELOVL5), ELOVL6, and ELOVL fatty acid elongase 7 (ELOVL7) were differentially expressed and identified in the crop tissues of lactating and non-lactating pigeons in our previous study [28]. However, only ELOVL6 appeared in the ceRNA network in the current study. ELOVL4 could be capable of elongating both saturated fatty acids and polyunsaturated fatty acid (PUFA) [71] and was responsible for the biosynthesis of saturated very long-chain fatty acids (VLC-FA) that are components of sphingolipids and ceramides and are important in the skin [73,74]. ELOVL6 and ELOVL7 can catalyze the elongation of saturated and monounsaturated fatty acids in mammals. ELOVL6 deficiency in mouse livers changes the composition of liver fatty acids, the length of the fatty acid chain, and the ratio of fatty acids, with a consequent reduction of SREBP-1 and PPARα in the liver. Decreased expression of SREBP-1 reduced fatty acid synthesis, but increased insulin sensitivity [75]. In this study, ELOVL6 was up-regulated in the crop of lactating pigeon, suggesting that it may be involved in fat synthesis.
It has been established that crop milk was formed mainly from proliferation and shedding of a large number of epithelial cells in the crop tissue. Indian Hedgehog (IHH) could regulate formation and proliferation of mesenchymal cells, which in turn affect epithelial proliferation and differentiation in intestinal epithelial cell [76]. Grainy head-like transcription factor 2 (GRHL2) was a member of grainy head-like transcription family, which played a fundamental role in epidermal integrity, embryonic neural tube closure, and wound healing processes [77,78]. ESRP2 protein was reported to regulate alternative splicing in epithelial cells [79,80]. Its loss disrupted the splicing of the transcripts from genes, thereby losing the ability to form sheets of cells with junctions between them [81]. In addition, ESRP2 and microtubule-associated serine/threonine kinase like (MASTL) clustered with epithelial markers and correlated inversely with the expression of mesenchymal markers [79,81]. EMT could activate a molecular program that drives epithelial cells to acquire mesenchymal phenotypes; it usually ends with the transition from mesenchymal to epithelial (MET). This state transition was very flexible [82]. Interestingly, this seems to be similar to the transformation between lactating and non-lactating crop tissues, and we infer that the ESRP2 gene was transformed between lactating and non-lactating crop tissues during the two periods. More importantly, in many cancers, the property of tumor dissemination and metastasis is closely associated with re-enabling developmental properties such as EMT [83]. These results indicated that the formation of crop milk is related to the transformation of cell morphology. Kinesin family member 26B (KIF26B) stabilized the cell-cell adhesion of mesenchymal cells in the developing kidney through enhancing the interaction of NMHC II and actin [84]. DOT1L between renal fibroblast activation and epithelial mesenchymal played a role in the process of transformation [85]. Phosphatidylinositol-3-kinase (PI3K) was a key signaling hub in immune cells. Either too little or too much PI3K activity was deleterious [86]. Mutations, activation, and low expression of PIK3CD were associated with autoimmune diseases and immunodeficiency [87,88].
In summary, a regulatory network was constructed. The regulatory relationship of these key genes has been listed and could be used as a key research object to further explain the production mechanism of crop milk. However, a limitation of this study is that most of the ceRNA networks identified were largely based on correlation analysis. Further strong evidence should be obtained to fully elucidate the identified ceRNA regulatory networks. Collectively, our study offered new insight for the ceRNA regulatory in crop milk production.

Conclusions
MiRNAs were abundant in the pigeon crop. Substance synthesis and cell morphogenesis genes remained active in lactation under the miRNA-associated-ceRNA regulation. Our findings enriched the pigeon miRNA expression profile and provide novel insights into the miRNA-associated ceRNA expression pattern and regulatory roles in crop milk production.

Conflicts of Interest:
The authors declare no conflict of interest.