Animal care in this study was performed in accordance with the Animal Experiment Management Regulations (Ministry of Science and Technology of China, 2004) approved by the Animal Care and Use Committee of Henan Agricultural University, China.

The Gushi chickens used as the experimental animals in this study were obtained from the Animal Center of Henan Agricultural University. A total of 200 one-day-old female Gushi chickens were raised in cages in the same environment under standard conditions used for the pure breeding conservation of Gushi chickens. In this study, three healthy chickens were randomly selected at 6, 14, 22, and 30 weeks of age, and twelve chickens were thus used in this study. The feeding method is the same as our previous research (Li Y. et al., 2021). The chickens were anesthetized by an intravenous injection of sodium pentobarbital (40 mg/kg) at a concentration of 0.2% into the wing vein. Under deep anesthesia, the chickens were euthanized by an intravenous injection of KCl (1–2 mg/kg), after which their abdominal adipose tissues were harvested, immediately frozen in liquid nitrogen and stored at −80°C until total RNA extraction.

Total RNA was extracted from the abdominal adipose tissues of Gushi chickens aged 6, 14, 22, and 30 weeks. The quality of RNA was evaluated by Microspectrophotometer and Agarose gel electrophoresis, and then used to construct cDNA library. An initial total amount of 3 μg of RNA was used for the construction of each library. First, ribosomal RNA was removed with the Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, Madison, Wisconsin, United States), and sequencing libraries were then generated using rRNA-depleted RNA with the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, United States). The products were then purified using a TruSeq RNASample Prep Kit v2 (New England Biolabs, Ipswich, MA, United States), and library quality was assessed on an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Finally, the libraries were sequenced using an Illumina HiSeq 2500 platform, and paired-end reads were generated.

Quality control of the sequencing data was conducted using fastp, and the reads of the adapter, the reads of multimers and the low-quality reads were removed to obtain clean reads (Chen et al., 2018). At the same time, the Q20, Q30, and GC contents of the pure data were calculated. All downstream analyses were based on high-quality clean data. Reference genome and gene model annotation files were downloaded from a genome website (ftp://ftp.ensembl.org/pub/release-104/fasta/gallus_gallus/dna/). An index of the reference genome was constructed using Hast2, and paired-end clean reads were aligned to the reference genome using Hast2 (Pertea et al., 2016). LncRNA expression was quantified using FeatureCounts software (Liao et al., 2014), and StringTie was used for the identification and quantification of transcripts from the aligned RNA-seq reads (Pertea et al., 2015). The transcripts obtained from the splicing of each sample were merged by using StringTie merge software, and transcripts with uncertain chain directions whose length did not exceed 200 nt were removed. The transcripts were then compared with the annotation of Ensembl version 104 using gffcompare (Pertea and Pertea, 2020), and the known coding protein transcripts were filtered out. Finally, three coding potential software programs, CPC2, PfamScan, and CNCI, were employed to screen for coding potential and thereby identify novel lncRNAs (Li et al., 2015; Kang et al., 2017; Luo et al., 2017).

FeatureCounts software was used to calculate the number of lncRNAs in each sample, and StringTie was used to determine their expression levels based on FPKM values. Differentially expressed gene (DEG) analysis was performed using DESeq2 software V1.22.2. The thresholds of a log2 fold change (FC) value > 1.5 and a q-value < 0.05 indicated upregulated DE-lncRNAs, while a log2 (FC) value < −1.5 and a q-value < 0.05 indicated downregulated DE-lncRNAs (Love et al., 2014).

The cis- and trans-regulatory target genes of the DE-lncRNAs were predicted. Among them, the DEGs located approximately 100 kb upstream and downstream of the DE-lncRNAs were used as cis-acting target genes. In addition, a Pearson correlation coefficient (r) > 0.7 was used to predict the regulation of trans target genes by lncRNAs. Subsequently, gene ontology (GO) and KEGG pathway analyses of the predicted target genes of each group of DE-lncRNAs were performed using the Clusterprofiler package, and GO terms and KEGG pathways with a corrected p-value (q value) < 0.05 were considered to be significantly enriched (Yu et al., 2012).

The transcriptome (PRJNA551368) and small RNA (PRJNA528858) library were constructed using the same RNA samples collected from the abdominal adipose tissues of Gushi chickens at the four developmental stages in our previous research. The association analysis of DE-lncRNAs was performed using these data (Chen Y. et al., 2019). miRanda and TargetScan software were used to predict the target interaction relationships between DE-lncRNAs and miRNAs (Agarwal et al., 2015). The differentially expressed lncRNA-miRNA-mRNA pairs were identified by Pearson’s correlation analysis. Combined with functional annotation, the lncRNA-miRNA mRNA pairs associated with lipid metabolism were selected, and Cytoscape software was used to construct their interaction network (Shannon et al., 2003).

To verify the accuracy of the data obtained by high-throughput sequencing, four lncRNAs were randomly selected, and the results were confirmed by qRT-PCR for nine randomly selected lncRNAs. The RNAs were reverse transcribed into cDNA with HiScript II Q RT SuperMix for qRT-PCR (+gDNA wiper) (Vazyme, Nanjing, China). qRT-PCR was performed on a LightCycler^® 96 qRT-PCR system (Roche, Basel, Switzerland) with ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China). β-actin were used as internal reference genes for mRNA. The relative expression levels were analyzed with the 2^−ΔΔCt method. Experiments were repeated at least three times, and all data are presented as the mean ± standard error of the mean. Data were analyzed and plotted with Prism software (GraphPad, San Diego, CA, United States). The qRT-PCR primer sequences are listed in Supplementary Table S1.

The data were analyzed by single factor analysis of variance (ANOVA) with SPSS26.0. The difference significance test was carried out by Duncan multiple comparison. p < 0.05 was taken as the significant difference.

A total of 12 cDNA libraries were constructed using abdominal adipose tissues from Gushi chickens at four developmental stages. After sequencing and strict quality evaluation, each library contained 13.53–16.88 GB of clean bases, with the CG content ranging from 46.20 to 51.56%. Comparative analysis showed that 91.97–95.55% of the reads (83.08–90.16% of the number of species) were located in the reference genome. Details regarding the sequencing and quality evaluation of each library are presented in Supplementary Table S2.

StringTie merge software was used to merge the transcripts sequenced from each library, and the merged transcripts were then screened and annotated. Finally, a total of 5,466 known lncRNAs and 3,827 new lncRNAs were identified from the 12 libraries (Figure 1A). Among them, 9,427 lncRNA transcripts were located in the intergenic region, 1,061 lncRNA transcripts were located in introns, 1,200 lncRNA transcripts were antisense lncRNAs, and 79 lncRNAs were sense-strand overlapping transcripts (Figure 1B). In addition, structural comparison and analysis of the identified lncRNAs and mRNAs revealed that the lncRNAs were shorter and had fewer exons than the mRNAs (Figures 1C,D).

To gain insight into the key lncRNAs involved in chicken abdominal fat development, we analyzed the DE-lncRNAs (|fold change, FC| ≥ 1.5, q-value < 0.05) at four different developmental stages, namely, 6 weeks (W6), 14 weeks (W14), 22 weeks (W22), and 30 weeks (W30). Among the five different comparison groups, there were 62, 79, 140, 58, and 62 DE-lncRNAs in the W14 vs. W6, W22 vs. W6, W30 vs. W6, W30 vs. W14, and W30 vs. W22 comparison groups, respectively. However, DE-lncRNAs were not identified in the W22 vs. W14 comparison group. Venn diagram analysis showed that no DE-lncRNAs were common among the five comparison groups (Figure 2A).

LncRNAs MSTRG.16790 and MSTRG.13152 were the common DE-lncRNAs in four comparison groups, W22 vs. W6, W30 vs. W6, W30 vs. W14, and W30 vs. W22; lncRNA ENSGALG00000047342 was common among the W14 vs. W6, W22 vs. W6, W30 vs. W14, and W30 vs. W22 groups; and the other four combinations shared DE-lncRNAs. In particular, MSTRG.11907, MSTRG.2857, MSTRG.18130, MSTRG.9971, and ENSGALG00000053669 showed the same temporal expression characteristics in the W30 vs. W6, W30 vs. W14, and W30 vs. W22 comparisons. In addition, cluster analysis revealed the DE-lncRNAs in each group that clustered together, and the within-group differences were smaller than the between-group differences (Figure 2B). In particular, the 22- and 14-week-old clusters were closest, which was consistent with the result that no DE-lncRNAs were identified between 14 and 22 weeks of age. In addition, we randomly selected 4 DE-lncRNAs for qRT-PCR verification, and the results revealed an expression pattern similar that obtained by RNA-seq sequencing, which indicated that the RNA-seq sequencing data were authentic and reliable (Figure 3).

To further elucidate the potential roles of the DE-lncRNAs in the development of abdominal fat in Gushi chickens, their cis- and trans-regulatory target genes were predicted, and functional enrichment analysis was performed. GO enrichment analysis results showed that the cis-regulated target genes showed significant enrichment for GO terms such as translation regulation, posttranscriptional regulation of gene expression, negative regulation of cell proliferation, and cell adhesion (Figures 4A–C). The trans-regulated target genes were significantly enriched for GO terms such as peptidase inhibitor activity and endopeptidase inhibitor activity (Figures 4D–F). In addition, KEGG pathway analysis revealed that the cis-regulated target genes showed significant enrichment for pathways related to glycosaminoglycan biosynthesis-keratan sulfate, glycosphingolipid biosynthesis-lactose and neolacto series, and other types of O-glycan biosynthesis (Figures 5A–C). The trans-regulated target genes showed significant enrichment for the PPAR signaling, fatty acid metabolism, fatty acid degradation, fatty acid biosynthesis, and adipocytokine signaling pathways (Figures 5D–F).

Based on the prediction results for the DE-lncRNA target genes, 9,207 cis-interaction regulatory relationships were observed between 136 DE-lncRNAs and 1,182 mRNAs, and 13,970 trans-interaction regulatory relationships were observed between 276 DE-lncRNAs and 559 mRNAs. Based on this, interaction regulatory networks between the DE-lncRNAs and their target genes was constructed. Many of these networks were related to lipid metabolism. In particular, the PPAR signaling pathway, which includes 10 genes, was significantly enriched in the W22 vs. W6, W30 vs. W6, and W30 vs. W22 comparison groups. In total, 275 lncRNA-mRNA interaction pairs were formed with 150 DE-lncRNAs (Figure 6), yielding a complex regulatory network. Among them, lncRNA ENSGALG00000037616 was shown to target the mRNAs of ACOX1, ADIPOQ, CPT1A, FABP5 and other genes in the PPAR signaling pathway. These lncRNA-mRNA regulatory networks may play an important role in fatty acid metabolism and abdominal fat deposition in chickens.

Based on the miRNA transcriptome profile data of abdominal fat samples from Gushi chickens at the four developmental stages, 51 differentially expressed miRNAs (DE-miRNAs) were identified, and the target relationships between the DE-lncRNAs and DE-miRNAs were predicted. In total, 4,058 lncRNA-miRNA interaction pairs (Supplementary Table S3) were observed between 272 DE-lncRNAs and 50 DE-miRNAs. Among these miRNAs, gga-miR-200b-3p, gga-miR-130b-3p, gga-miR-215-5p, gga-miR-122-5p, gga-miR-223, and gga-miR-125b-5p were previously shown to be related to lipid metabolism (Kim et al., 2013; Tao et al., 2016; Wei et al., 2017; Du et al., 2018; Xu et al., 2018; Li W. D. et al., 2019; He et al., 2019; Rodríguez et al., 2020; Li G. et al., 2021). To this end, we built a miRNA-lncRNA (FPKM > 10) interaction regulatory network based on these miRNAs (Figure 7).

Based on the whole-transcriptome data association analysis of abdominal adipose tissue samples from Gushi chickens at the four developmental stages, 1,740 differentially expressed lncRNA-miRNA-mRNA pairs (Supplementary Table S4) were screened to construct a ceRNA regulatory network related to abdominal fat development. In particular, we found a ceRNA regulatory network that was related to lipid metabolism (Figure 8), which contained 82 lncRNA-miRNA-mRNA pairs formed by 51 DE-lncRNAs, 4 DE-miRNAs and 5 DE-mRNAs. Meanwhile, we used lncLocator to predict the subcellular localization of 51 lncRNAs involved in the construction of ceRNAs (Cao et al., 2018). The results showed that 48 lncRNA predictions were located in the cytoplasm and another 4 lncRNAs were predicted to be located in Exosome and Cytosol. These results show that most of the lncRNAs we used to construct the ceRNA networks are located in cytoplasm (Supplementary Table S5). Among the components of this network, ASBG2, LPL, and FABP5 are involved in fatty acid transport. They are related to three miRNAs (gga-miR-122-5p, ga-miR-125b-5p, and gga-miR-200b-3p), the latter two of which form 51 ceRNAs and 36 lncRNAs. In addition, PLIN2 was related to adipocyte differentiation and forms 11 ceRNAs with gga-miRNA-125b-5p and 11 lncRNAs. In addition, PPARD, gga-miR-130b-3p and 16 lncRNAs also formed 16 ceRNAs.