MicroRNAs and their regulatory networks in Chinese Gushi chicken abdominal adipose tissue during postnatal late development

Background Abdominal fat is the major adipose tissue in chickens. The growth status of abdominal fat during postnatal late development ultimately affects meat yield and quality in chickens. MicroRNAs (miRNAs) are endogenous small noncoding RNAs that regulate gene expression at the post-transcriptional level. Studies have shown that miRNAs play an important role in the biological processes involved in adipose tissue development. However, few studies have investigated miRNA expression profiles and their interaction networks associated with the postnatal late development of abdominal adipose tissue in chickens. Results We constructed four small RNA libraries from abdominal adipose tissue obtained from Chinese domestic Gushi chickens at 6, 14, 22, and 30 weeks. A total of 507 known miRNAs and 53 novel miRNAs were identified based on the four small RNA libraries. Fifty-one significant differentially expressed (SDE) miRNAs were identified from six combinations by comparative analysis, and the expression patterns of these SDE miRNAs were divided into six subclusters by cluster analysis. Gene ontology enrichment analysis showed that the SDE miRNAs were primarily involved in the regulation of fat cell differentiation, regulation of lipid metabolism, regulation of fatty acid metabolism, and unsaturated fatty acid metabolism in the lipid metabolism- or deposition-related biological process categories. In addition, we constructed differentially expressed miRNA–mRNA interaction networks related to abdominal adipose development. The results showed that miRNA families, such as mir-30, mir-34, mir-199, mir-8, and mir-146, may have key roles in lipid metabolism, adipocyte proliferation and differentiation, and cell junctions during abdominal adipose tissue development in chickens. Conclusions This study determined the dynamic miRNA transcriptome and characterized the miRNA–mRNA interaction networks in Gushi chicken abdominal adipose tissue for the first time. The results expanded the number of known miRNAs in abdominal adipose tissue and provide novel insights and a valuable resource to elucidate post-transcriptional regulation mechanisms during postnatal late development of abdominal adipose tissue in chicken.


Background
The chicken is one of the most important agricultural animals [1]. Abdominal fat is the major adipose tissue in chickens. In poultry production, the characteristics of abdominal adipose tissue are closely related to the economic traits of chickens [2]. Excessive accumulation of abdominal fat often affects meat quality. Moreover, chickens can also be used as a biomedical model for studies on human metabolic disorders, such as insulin resistance, diabetes, obesity and metabolic syndrome [3,4]. Therefore, understanding the molecular mechanisms of abdominal adipose tissue development can provide insight into biomedical research and benefit poultry genetic breeding and production in chickens. However, further research investigating the molecular mechanisms of chicken abdominal adipose tissue development is warranted.
The development of adipose tissue is a result of both adipocyte hyperplasia and hypertrophy. This complex process is regulated by a series of cellular and molecular events [5,6]. MicroRNAs (miRNAs) are a class of~22nucleotide (nt)-long endogenous small noncoding RNAs that regulate gene expression at the post-transcriptional level in organisms; thus, they play an important role in many biological processes, including development, metabolism, and differentiation. A number of studies have shown that miRNAs also participate in the biological processes involved in adipose tissue development, including adipocyte proliferation and differentiation, adipogenesis, and lipid metabolism [7]. For example, miRNAs such as miR-143, miR-17-92 cluster, miR-21, miR-204, miR-378, and miR-375 have been shown to promote adipogenesis [8][9][10][11][12][13], while miRNAs such as miR-27b, let-7, miR-22, and miR-130 have been reported to impair adipocyte differentiation [14][15][16][17]. Therefore, miRNAs have become novel targets for investigating the molecular mechanisms related to adipose tissue development in animals [18][19][20]. Over the past decade, miRNAs in adipose tissues have been extensively studied. These studies have focused on the identification, expression, and function of miRNAs in adipose tissues of different anatomical depots [21][22][23][24][25], and these findings have helped to elucidate the molecular mechanisms underlying adipose tissue development in animals. However, previous studies were limited to humans, mice and a few agricultural animals, as well as metabolic diseases, such as obesity [26,27].
Although 1232 miRNAs have been identified in the chicken genome (miRBase Version 22.1, October 2018), to date, there have been notably few reports on abdominal adipose tissue development in chickens [28]. In an earlier study, only 47 miRNAs were identified by molecular cloning experiments in abdominal adipose tissue from 28-day-old Arbor Acres commercial chickens [29]. More recently, Huang et al. (2015) identified 230 known miRNAs and 83 potentially novel miRNAs by deep sequencing in abdominal adipose tissues from a Chinese local breed and a commercial broiler line at 93 days of age and determined 62 differentially expressed miRNAs [30]. Two other studies detected 159 known miRNAs [31] and 33 differentially expressed miRNAs [32] in preadipocytes obtained from chicken abdominal adipose tissues, respectively. These studies helped us to understand miRNA functions in the development of chicken abdominal adipose tissues, but they also had the following limitations: i) they mainly concentrated on a few breeds and a specific stage of abdominal fat development; ii) they lacked an overall understanding of the regulatory mechanism of miRNAs during abdominal fat development in chickens as miRNAs function as an interactive network; and iii) they only identified a relatively small number of miRNAs using early experimental techniques. It is wellknown that the physiological characteristics of abdominal adipose tissue are distinct between different stages of development in chickens post-hatching. The growth status of abdominal fat during postnatal late development ultimately affects meat yield and quality in chickens. However, there is still a paucity of knowledge on the regulatory mechanisms of miRNAs underlying abdominal adipose tissue development in the post-hatch late stage in different breeds. For these reasons, the identification of more miRNAs, construction of miRNA dynamic expression profiles, and further research investigating the interaction regulatory network of miRNAs are necessary to reveal the miRNA regulatory mechanism in chicken abdominal adipose tissue during postnatal late development.
The Gushi chicken is a domestic Chinese breed, and it is often used for breeding and production. To reveal the miRNA regulatory mechanism in the postnatal late development of abdominal fat in this variety, we constructed dynamic expression profiles of miRNAs in abdominal adipose tissue from Gushi chicken at 6, 14, 22, and 30 weeks after hatching. In addition, we also established the miRNA regulatory network involved in abdominal fat development based on the integrated analysis results of miRNA and mRNA expression profiles in Gushi chicken. The results expanded the number of known miRNAs in abdominal adipose tissue in chickens and provided a better understanding and novel insights into the miRNA regulatory mechanism of abdominal adipose tissue development in the postnatal late development stage of chickens.

Small RNA library sequencing and sequence analysis
Using Illumina HiSeq2500 sequencing, we obtained 12, 316,628, 26,800,794, 13,427,210 and 14,073,376 raw reads from the four small RNA libraries of Gushi chicken abdominal adipose tissue at 6 (z06), 14 (z14), 22 (z22), and 30 (z30) weeks old, respectively ( Table 1). The raw  sequence reads were deposited in the NCBI database Sequence Read Archive with the accession number  PRJNA528858. After filtering, 11,220,507, 25,261,754, 12,  928,203 and 13,520,406 high-quality clean reads measuring 18-35 nt were finally obtained from the z06, z14, z22 and z30 libraries, respectively. The lengths of these small RNAs were mainly 21-24 nt. In particular, the 22-nt small RNAs accounted for 47.35-58.54% of the total number of small RNAs in the libraries (Additional file 1: Figure S1). In addition, the 18-35-nt small RNA sequences were aligned with the chicken genome sequence. The results showed that 80.95, 85.63, 88.29 and 87.02% of the sequences from the z06, z14, z22 and z30 libraries were perfectly mapped to the chicken genome sequence, respectively (Table 1). These mapped sequences were used for subsequent miRNA identification.

MiRNAs expressed in Gushi chicken abdominal adipose tissues
By category annotation, the tRNA, rRNA, snoRNA, and other snRNA sequences in the clean reads were eliminated (Additional file 7: Table S1). The remaining sequences were then compared to the mature sequences of miRNAs from chickens in miRBase (Release 22.0). The 11,460 unique sequences (24,059,167 reads) were fully matched to the miRNAs of the chicken. Combining structural predictions of precursor sequences, a total of 507 known miRNAs and 53 novel miRNAs were identified from the four small RNA libraries in Gushi chicken abdominal adipose tissue (Additional file 11). Detailed information on novel miRNAs is shown in Additional file 8: Table S2. The statistics for the miRNAs identified in each library are shown in Table 2.
In addition, we also analysed the expression levels of 560 identified miRNAs in different developmental stages of Gushi chicken abdominal adipose tissue. Among these miRNAs, 312 miRNAs were expressed in all four developmental stages, while 29, 54, 18, and 33 miRNAs were expressed specifically at 6, 14, 22, and 30 weeks of age, respectively ( Fig. 1 and Additional file 11). This finding indicated that there were differences in the miRNA expression profiles during different developmental stages of chicken abdominal fat. In addition, the abundance of these identified microRNAs was also different. Most of the novel miRNAs were relatively weakly expressed in the abdominal adipose tissue samples, while the abundance of miRNAs expressed specifically in different developmental stages of abdominal adipose tissue was also considerably lower (Additional file 11). The percentages of miRNAs with TPM > 60 in z06, z14, z22, and z30 were 18.38, 19.29, 19.64, and 18.57%, respectively, which suggests that only a few miRNAs were abundantly expressed during abdominal adipose tissue development in chicken. Of these abundant miRNAs, some dominated the miRNA libraries in different developmental stages of abdominal fat. For example, the expression level of miR-148a-3p was the highest in all four developmental stages, and its reads in the z06, z14, z22, and z30 libraries accounted for 7.53, 12.95, 14.08, and 9.95% of total clean reads, respectively.

Target gene predictions and functional analysis
Based on the SDE miRNAs identified in different combinations, 3489, 3212, 3071, 1105, 3067 and 3113 target genes were predicted in the z14 vs. z06, z22 vs. z06, z30 vs. z06, z22 vs. z14, z30 vs. z14, and z30 vs. z22 combination, respectively. Therefore, 843 predicted target genes were shared among all six combinations, while the number of unique target genes in the z14 vs. z06, Z22 vs. z06, Z30 vs. z06, Z22 vs. z14 and Z30 vs. z22 combination were 21, 12, 110, 17 and 43, respectively (Additional file 2: Figure S2). Subsequently, the predicted target genes of SDE miRNAs were analysed by GO enrichment. In the biological process category, 108, 106, 74 respectively. Most of these GO terms did not differ significantly between combinations. Their functions were mainly concentrated in the biological process related to localization, metabolic processes and transport (Additional file 3: Figure S3), for instance, cellular macromolecule localization, cellular metabolic processes, cellular protein localization, cytoplasmic transport, intracellular transport, macromolecule localization, and primary metabolic processes. However, the SDE miRNA biological processes were also different in the different developmental stages of Gushi chicken abdominal adipose tissue. The significantly enriched GO terms in the z14 vs. z06 combination were mainly involved in the cell cycle process, regulation of the apoptotic process, and regulation of cellular component organization, and the GO terms in the z22 vs. z14 combination were mainly related to Golgi vesicle transport and macromolecular complex disassembly, while the GO terms in the z30 vs. z22 combination were mainly focused on cellular protein complex assembly, fatty acid metabolism, and regulation of kinase activity ( Fig. 4 and Additional file 13). Moreover, we also found many GO terms related to lipid metabolism or fat deposition in the GO enrichment results, although these GO terms were not significantly enriched within combinations. In particular, some important biological processes involved in the regulation of fat cell differentiation, lipid metabolism, fatty acid metabolism and unsaturated fatty acid metabolism were the same in the three combinations, z14 vs. z06, z22 vs. z14 and z30 vs. z22 (Additional file 10: Table S4). Only a few biological processes were different, including acylglycerol homeostasis, medium-chain fatty acid metabolism, and lipid hydroperoxide transport (Additional file 4: Figure S4). These results suggested that the basic biological processes in which SDE miRNAs participated were similar in the development of Gushi chicken abdominal adipose tissue from 6 weeks to 30 weeks, but there were differences in key biological processes at different developmental stages.

Discussion
Approximately 25-70% of body fat distribution changes are influenced by genetic background in animals. In chickens, abdominal fat weight (AbFW) and abdominal fat percentage (AbFP) are often used as phenotypic indices of fat traits. The heritability coefficients (h 2 ) are 0.62 for AbFW and 0.24 for AbFP. However, these traits are genetically different in different chicken breeds and lines. In fact, the h 2 of abdominal fat mass varies from 0.5 to 0.8 in chicken [33]. This finding suggests that it should be feasible to modulate abdominal fat deposition in chickens using genetic control strategies. In poultry production, excessive abdominal fat is often removed as a poultry by-product because it affects meat quality. Today, genetic selection for low-abdominal fat has become the main goal of meat-type chicken breeding. Therefore, investigations into the molecular mechanisms of abdominal fat deposition in chickens are important for increasing productivity. In the last few decades, there have been extensive studies on the cellular basis, as well as regulatory gene identification and the expression profile of abdominal adipose tissue in chickens [2,[33][34][35][36][37][38]. However, knowledge of post-transcriptional regulatory mechanisms underlying abdominal adipose tissue development in chickens, particularly the role and function of miRNAs, is still very limited. In this study, we constructed dynamic expression profiles of miRNAs in the abdominal adipose tissue of Gushi chickens at 6, 14, 22, and 30 weeks and identified 507 known miRNAs and 53 novel miRNAs. These results expand the number of miRNAs expressed in abdominal adipose tissue and provide a valuable resource for the selective breeding of abdominal fat traits in chickens.
According to our previous studies, the abdominal fat rates (i.e., abdominal fat weight/live weight× 100%) of Gushi chickens at 6, 14, 22 and 30 weeks were 0.63 ± 0.42, 1.95 ± 1.08, 3.25 ± 1.26 and 2.96 ± 0.35, respectively. In particular, the period from 14 weeks to 22 weeks is the rapid deposition stage of abdominal fat in Gushi chicken. Functional enrichment analysis of differentially expressed miRNAs was completed in this study. The results showed that the enriched biological processes were Fig. 5 The miRNA-mRNA interaction networks containing eleven lipid metabolism or deposition related pathways. The dot indicates miRNA, and the box indicates a target gene that had a negative correlation with a given miRNA. Red indicates up-regulated, green indicates down-regulated, and yellow indicates a miRNA that was up-regulated in a given comparison combination and was down-regulated in other comparison combinations. The eleven pathways are fatty acid degradation, fatty acid metabolism, biosynthesis of unsaturated fatty acids, alpha-linolenic acid metabolism, fatty acid biosynthesis, arachidonic acid metabolism, fatty acid elongation, glycerophospholipid metabolism, glycerolipid metabolism, the adipocytokine signaling pathway and the PPAR signaling pathway. The miRNA-mRNA pairs in each pathway are shown in Additional file 9: Table S3 significantly different at different developmental stages of Gushi chicken abdominal fat. The biological processes comprising the regulation of apoptotic process, cell cycle process, and regulation of cellular component organization were mainly enriched between 6 weeks and 14 weeks, and the biological process of macromolecular complex disassembly was enriched between 14 weeks and 22 weeks, while the biological process of fatty acid metabolism was significantly enriched between 22 weeks and 30 weeks (Fig. 4). This finding suggests that 14 weeks and 22 weeks are important stages in abdominal fat deposition in Gushi chicken. Generally, the abdominal fat mass is determined by the number and volume of adipocytes in the depot [39]. It is speculated that in Gushi chicken, the abdominal fat pad and adipocyte hyperplasia dominate before 14 weeks, and adipocyte hypertrophy dominates between 14 weeks and 22 weeks, while growth of abdominal fat is primarily the filling of existing adipocytes with lipids after 22 weeks. The dynamic change in abdominal fat deposition in Gushi chicken is largely consistent with the cellular basis of abdominal adipose tissue development described in previous studies [2]. Therefore, the miRNA dynamic expression profiles obtained in this study reflect the mechanisms of molecular regulation of abdominal adipose tissue development at the post-transcriptional level in Gushi chicken.
Previous studies have confirmed that many miRNAs are expressed in a tissue-specific manner. In this study, 45 miRNA families with more than two members were detected, and most of the members in these miRNA families were expressed in high abundance in Gushi chicken abdominal adipose tissue. The most abundant miRNAs were the let-7 family, of which nine members were detected in Gushi chicken abdominal fat. In particular, six members, including let-7a-5p, let-7b, let-7f-5p, let-7 g-5p, let-7i, and let-7 k-5p, were expressed in high abundance in abdominal fat at four developmental stages. Let-7 has been reported to directly target the high-mobility group AT-hook 2 (HMGA2) and to inhibit adipocyte differentiation [15]. Of these abundant miR-NAs, the miRNA with the highest abundance was gga-miR-148a-3p. MiR-148a has previously been shown to promote adipogenesis by targeting WNT1 [40]. It is Fig. 6 The miRNA-mRNA interaction networks containing five cell proliferation and differentiation related pathways. The dot indicates miRNA, and the box indicates a target gene that had a negative correlation with a given miRNA. Red indicates up-regulated, green indicates downregulated, and yellow indicates a miRNA that was up-regulated in a given comparison combination and was down-regulated in other comparison combinations. The five pathways are FoxO signaling pathway, MAPK signaling pathway, TGF-beta signaling pathway, ErbB signaling pathway and Wnt signaling pathway. The miRNA-mRNA pairs in each pathway are shown in Additional file 9: Table S3 Fig. 7 The miRNA-mRNA interaction networks containing two matrix connectivity (a) and (b) represent the miRNA-mRNA interaction networks with three cell junction (focal adhesion, ECM-receptor interaction, gap junction) and two matrix connectivity (tight junction, and adherens junction) related pathways, respectively. The dot indicates miRNA, and the box indicates a target gene that had a negative correlation with a given miRNA. Red indicates up-regulated, green indicates down-regulated, and yellow indicates a miRNA that was up-regulated in a given comparison combination and was down-regulated in other comparison combinations. The five pathways are focal adhesion, ECM-receptor interaction, gap junction, tight junction, and adherens junction. The miRNA-mRNA pairs in each pathway are shown in Additional file 9: Table S3 Chen et al. BMC Genomics (2019) 20:778 known that the miR-17-92 miRNA cluster promotes adipocyte differentiation by targeting RB2/p130 [9]. As a member of the miR-17-92 cluster, miR-17-3p was previously found to enhance 3 T3-L1 differentiation by targeting fatty acyl desaturase genes [41]. Interestingly, six members of the mir-17 family were found in Gushi chicken abdominal adipose tissue. Similarly, other abundant miRNAs, miR-30a-5p, miR-146b-5p, miR-30d, miR-21, miR-101-3p, and miR-27b-3p, have been shown to regulate adipogenesis in other species [10,14,[42][43][44]. Moreover, some miRNAs that showed developmental stage-specificity were also found in Gushi chicken abdominal adipose tissue (Additional file 11). Although they were expressed in low abundance in Gushi chicken abdominal fat, some of them, such as miR-103, miR-31, and miR-135a, have previously been reported to participate in the regulation of adipocyte differentiation or adipogenesis [45][46][47][48][49]. These results suggest that developmental stage-specificity and high-abundance miRNAs may play an important role during abdominal adipose tissue development in chickens. Many miRNAs are expressed in a spatiotemporal-specific manner in organisms. In this study, 51 significantly differentially expressed (SDE) miRNAs were screened from Gushi chicken abdominal adipose tissue at four developmental stages. These SDE miRNAs showed obvious temporal expression characteristics, and their expression patterns were divided into six types (Fig. 2). Specifically, three miRNAs, miR-215-5p, miR-122-5p and miR-499-5p, were differentially and abundantly expressed across the four stages during abdominal adipose tissue development. MiR-122 is a liver tissue-specific miRNA and plays a role in lipid metabolism [50]. Recent research has shown that miR-122 overexpression induces hepatic differentiation of adipose tissue-derived stem cells [51]. In addition, miR-122 is also the most significant signature between visceral fat and subcutaneous fat, and its abundance affects PPAR-γ signalling and adipocyte differentiation in vitro and in human adipose tissues [52]. Although the function of miR-499 in lipid metabolism has not been reported, a recent study has shown that miR-499 can regulate PRDM16 to affect insulin-induced skeletal muscle satellite cell (SMSC) adipogenic differentiation [53]. Moreover, previous studies have suggested that miR-215-5p is a negative regulator during early adipogenesis in 3 T3-L1 cells [54]. For the above three miRNAs, further studies are required to clarify the relationship between their expression and fat-related phenotype and between the modulations of their expression and fat content, which will ultimately determine whether they can serve as molecular markers of chicken abdominal adipose tissue development. In addition to the above three miRNAs, other SDE miRNAs, miR-29b-3p, miR-1a-3p, miR-10a, miR-206, miR-429, miR-200a, miR-454, miR-200b, miR-31-5p, miR-204, miR-375, miR-155, miR-194, miR-130, and miR-365, have previously been shown to be related to adipogenesis [11,13,17,[55][56][57][58][59]. miR-375 has been reported to promote adipogenesis by extracellular signal regulated kinase (ERK1/2) signalling [13]; miR-204 has been found to improve adipogenesis by targeting runt-related transcription factor 2 (RUNX2) [11]. Conversely, miR-130 has been reported to directly target PPAR-γ and to inhibit the differentiation of human preadipocytes [17]. Similarly, miR-155 has been demonstrated to inhibit adipogenesis by targeting the cAMP response element binding protein (CREB) and C/EBPβ [59]. The above results indicate that these temporally and differentially expressed miRNAs are closely related to abdominal fat development in chickens.
In organisms, miRNAs construct a sophisticated control network and participate in the regulation of various biological processes through complex interactions with their target genes. To reveal the potential regulatory relationships of miRNAs underlying abdominal adipose tissue development in Gushi chicken, we carried out correlation analysis between miRNA and mRNA and performed KEGG pathway enrichment analysis of differentially expressed miRNA-mRNA pairs. In addition, based on the results of KEGG pathway analysis, we constructed miRNA-mRNA interaction networks involved in fatty acid metabolism, glycerolipid metabolism, cell junctions, and cell proliferation and differentiation during abdominal adipose tissue development in Gushi chicken (Figs. 5, 6 and 7). In these miRNA-mRNA interaction networks, some key genes as components of multiple pathways link multiple pathways to each other. For instance, ACSBG2, ACOX1 and FADS2 are found across five pathways, including fatty acid metabolism, fatty acid degradation, PPAR signalling pathway, biosynthesis of unsaturated fatty acids and alpha-linolenic acid metabolism, and ADIPOQ is found in the adipocytokine signalling pathway and PPAR signalling pathway. Similarly, THBS1 also connects three pathways: the focal adhesion, transforming growth factor (TGF)-beta signalling pathway and ECM-receptor interaction. In these interaction networks, some miRNAs regulate the functions of multiple pathways by targeting multiple genes. Thus, a complex regulatory network formed by the interactions between miRNAs and their target gene and between pathways regulates abdominal adipose tissue development in Gushi chicken. After exhaustive literature mining of the miRNAs constituting the abovementioned interaction network, some miRNAs, including miR-148a, miR-103, miR-31, miR-204, miR-375, miR-155, miR-130, miR-142, miR-101, miR-200, miR-29b, miR-9, miR-32, miR-222, miR-206, miR-1a, miR-146b, miR-181, miR-30, miR-22, miR-27, and miR-194, as well as let-7 family members and the miR-17-92 cluster, have previously been reported to be related to lipid metabolism and adipogenesis [9, 11, 13-17, 40, 42-45, 47, 55-57, 59-64].
Therefore, these interactive networks reflect the complexity of molecular regulation of abdominal adipose tissue development at the post-transcriptional level in chickens.

Conclusion
We have described the miRNA dynamic expression profiles and the miRNA-mRNA interaction networks involved in abdominal fat development in Gushi chicken. It was found that miRNAs regulate abdominal adipose tissue development, primarily by affecting lipid metabolism or deposition, adipocyte proliferation and differentiation, matrix interaction, and cell junctions, through a complex interaction network in chicken. These results provide novel insights and a valuable resource for a better understanding of the post-transcriptional regulatory mechanisms of abdominal fat development in chickens. The findings reported in this study will be further analysed and validated.

Animal feeding and sample collection
In this study, the experimental animals were Gushi chickens from the Animal Center of Henan Agricultural University. Female chickens were selected and raised in cages in the same environment with standard conditions. The diet comprised 18.5% crude protein and 12.35 MJ/ kg energy before 14 weeks, and 15.6% crude protein and 12.75 MJ/kg energy after 14 weeks. Three healthy individuals were randomly selected at the ages of 6 weeks (z06), 14 weeks (z14), 22 weeks (z22) and 30 weeks (z30). These individuals were euthanized by intravenous injection of KCl (1-2 mg/kg) under deep anaesthesia. The abdominal fat pad and adipose tissue around the stomach were isolated and frozen immediately in liquid nitrogen. The frozen sample was stored at − 80°C. Total RNA was isolated from abdominal adipose tissue using Trizol reagent (TaKaRa, China), and the mixed RNA samples from three individuals at each developmental stage were used for library construction.

Small RNA library construction
Four small RNA libraries were constructed from abdominal adipose tissue of Gushi chicken at 6, 14, 22 and 30 weeks old using the Small RNA Sample Prep Kit (Illumina). First, 16-30-nt small RNAs were added to 3′ and 5′ adapters, and cDNA was reverse transcribed. Subsequently, PCR amplification was carried out, and the target DNA fragment was separated by PAGE gel electrophoresis, and the final cDNA library was obtained after enrichment. The library was diluted to 1 ng/μl, and then the library's insert size was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Finally, the qualified cDNA library, i.e., The insert size met our expectation, and the effective concentration was more than 2 nM in the libraries and sequenced on the Illumina HiSeq 2500 sequencing platform.

Post-sequencing analysis
The raw reads were processed to obtain clean reads by removing reads containing poly-N with 5′ adapter contaminants without 3′ adapter or the insert tag containing poly-A or T or G or C and low-quality reads. Then, the 18-35-nt clean reads were mapped to the chicken genome by Bowtie software [82]. The reads mapped onto the chicken genome were selected and blasted against the GenBank and Ensembl ncRNA databases, as well as the RNA families in Rfam, to annotate and remove the non-coding RNA sequences. The remaining sequences were then used to identify known miR-NAs in the abdominal adipose tissue of Gushi chickens by comparing them with the mature chicken miRNAs in miR-Base (Release 22.0). Finally, to identify potentially novel miR-NAs in Gushi chicken abdominal fat, the remaining unannotated sRNA sequences were aligned against the chicken genome assembly (Release 4.0) from Ensemble to obtain genomic sequences containing sRNA. The hairpin structures and folding energy for the surrounding 300 bases flanking each small RNA sequence were predicted using miREvo [83] and miRDeep2 [84] software. The sequences with a typical stem-loop hairpin structure were considered potential novel miRNAs. These novel miRNAs were named as described by Ambros et al. (2003) [85].

Identification and GO enrichment analysis of SDE miRNAs
The miRNA expression levels at different developmental stages of Gushi chicken abdominal fat were estimated using the transcripts per million clean reads (TPM). The following formula was used to calculate the normalized expression values: (read count × 1,000,000)/total miRNA read count in four libraries. If the normalized expression of a given miRNA is zero in all samples, this miRNA is removed in future differential expression analysis. The differentially expressed miRNAs were analysed by DEGseq [86]. The fold change for each miRNA between two discretionary developmental stages was calculated based on the comparison of six combinations, i.e., z06/z14, z06/z22, z06/z30, z14/z22, z14/z30, and z22/z30. Fisher's test was used to determine the P-values, and the false discovery rate (FDR) method was introduced to determine the threshold p value in multiple tests [87]. MiR-NAs with a q value < 0.01 and | log 2 (fold change) | > 1 were considered SDE miRNAs.
To reveal the potential functions of SDE miRNAs during the development of Gushi chicken abdominal fat, miRanda [88] and TargetScan [89] were used to predict the potential target genes of the SDE miR-NAs. For the miRanda software, all detected targets with scores and energies less than the threshold parameters of S > 140 and ΔG < -10 kcal/mol and strict 5′ seed pairing were selected as potential targets. For the TargetScan software, the 2-8-nt sequences that start from 5′ small RNA were chosen as seed sequences to predict the 3'UTR of mRNA, and the predicted targets were ranked based on the predicted repression or aggregate P CT score of the longest 3′-UTR isoform. Finally, only the sites belonging to the top quartile of ranked predictions and present in several species, including Gallus gallus, were retained as true binding sites. To make the identification of target genes robust, miRNAs that were simultaneously predicted by two programs were selected for this study. Gene ontology (GO) enrichment analysis of the predicted targets was performed using the GOseq method [90]. To obtain a corrected P-value, the FDR method was used to adjust the P-value threshold of multiple tests [87]. When the corrected P-values were less than 0.05, GO terms were considered significantly enriched.

Interaction analysis between miRNAs and mRNAs
To analyse the interactions between miRNAs and mRNAs, we also determined the dynamic transcriptome profiles of abdominal fat using the same tissue samples and the small RNA libraries. A total of 12 cDNA libraries were constructed from Gushi chicken abdominal fat tissues at 6, 14, 22 and 30 weeks using NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) and were sequenced on an Illumina Hiseq 2500 platform. The raw sequences were deposited into the NCBI Sequence Read Archive under BioProject number PRJNA551368. Post-sequencing analyses, such as sequence analysis and differential expression analysis, were performed using the methods described by Li et al. (2019) [91].
Based on the abovementioned transcriptome profile data, the potential target genes of differentially expressed miRNAs were predicted using miRanda [88] and TargetScan [89] software according to the above described method, and these target genes were matched with the miRNAs. According to the expression levels of miRNAs and genes in Gushi chicken abdominal fat, the differentially expressed miRNA-mRNA pairs were identified by Pearson's correlation analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of all target genes in each combination was performed using KOBAS software [92] to identify the pathways associated with lipid metabolism and adipogenesis. Given the inverse correlation between miRNA and target expression patterns, the negatively correlated miRNA-mRNA pairs were selected from the above pathways, and their interaction networks were determined by Cytoscape 3.4 (http://www.cytoscape.org/).

Quantitative real-time PCR (qRT-PCR) analysis
To verify the accuracy of the sequencing data, a primer-Script™ RT reagent Kit (TaKaRa, Dalian, China) was used for reverse transcription, and BμLge-Loop™ miRNA qRT-PCR Primer specific for the differentially expressed miRNAs was designed by GenePharma (GenePharma, Shanghai, China). The qRT-PCR analysis was performed using the LightCycler® 96 instrument (Roche Applied Science). All reactions were repeated three times. With the U6 small nuclear RNA as an endogenous control, the relative expression levels were determined using the study design or in any aspect of the collection, analysis and interpretation of data or in writing the manuscript.

Availability of data and materials
All data generated or analysed in this study are provided in additional files. The raw sequences from this study have been deposited to the NCBI Sequence Read Archive under accession number PRJNA528858 (BioProject ID of miRNA) and PRJNA551368 (BioProject ID of mRNA).
Ethics approval and consent to participate Experimental and animal care is executed in accordance with the programme approved by the Experimental Animal Management Ordinance (Ministry of Science and Technology, 2004), and were approved by the Institutional Animal Care and Use Committee (Use Committee of Henan Agricultural University, China; Permit Number: 17-0118).

Consent for publication
Not applicable