Analysis of potential regulatory LncRNAs and CircRNAs in the oxidative myofiber and glycolytic myofiber of chickens

SART and PMM are mainly composed of oxidative myofibers and glycolytic myofibers, respectively, and myofiber types profoundly influence postnatal muscle growth and meat quality. SART and PMM are composed of lncRNAs and circRNAs that participate in myofiber type regulation. To elucidate the regulatory mechanism of myofiber type, lncRNA and circRNA sequencing was used to systematically compare the transcriptomes of the SART and PMM of Chinese female Qingyuan partridge chickens at their marketing age. The luminance value (L*), redness value (a*), average diameter, cross-sectional area, and density difference between the PMM and SART were significant (p < 0.05). ATPase staining results showed that PMMs were all darkly stained and belonged to the glycolytic type, and the proportion of oxidative myofibers in SART was 81.7%. A total of 5 420 lncRNAs were identified, of which 365 were differentially expressed in the SART compared with the PMM (p < 0.05). The cis-regulatory analysis identified target genes that were enriched for specific GO terms and KEGG pathways (p < 0.05), including striated muscle cell differentiation, regulation of cell proliferation, regulation of muscle cell differentiation, myoblast differentiation, regulation of myoblast differentiation, and MAPK signaling pathway. Pathways and coexpression network analyses suggested that XR_003077811.1, XR_003072304.1, XR_001465942.2, XR_001465741.2, XR_001470487.1, XR_003077673.1 and XR_003074785.1 played important roles in regulating oxidative myofibers by TBX3, QKI, MYBPC1, CALM2, and PPARGC1A expression. A total of 10 487 circRNAs were identified, of which 305 circRNAs were differentially expressed in the SART compared with the PMM (p < 0.05). Functional enrichment analysis showed that differentially expressed circRNAs were involved in host gene expression and were enriched in the AMPK, calcium signaling pathway, FoxO signaling pathway, p53 signaling pathway, and cellular senescence. Novel_circ_004282 and novel_circ_002121 played important roles in regulating oxidative myofibers by PPP3CA and NFATC1 expression. Using lncRNA-miRNA/circRNA-miRNA integrated analysis, we identified many candidate interaction networks that might affect muscle fiber performance. Important lncRNA-miRNA-mRNA networks, such as lncRNA-XR_003074785.1/miR-193-3p/PPARGC1A, regulate oxidative myofibers. This study reveals that lncXR_003077811.1, lncXR_003072304.1, lncXR_001465942.2, lncXR_001465741.2, lncXR_001470487.1, lncXR_003077673.1, XR_003074785.1, novel_circ_004282 and novel_circ_002121 might regulate oxidative myofibers. The lncRNA-XR_003074785.1/miR-193-3p/PPARGC1A pathway might regulate oxidative myofibers. All these findings provide rich resources for further in-depth research on the regulatory mechanism of lncRNAs and circRNAs in myofibers.

). Heatmaps generated from the expression of DE-mRNA, DE-lncRNA, and DE-circRNA were used to show the expression patterns of these RNAs between PMM and SART ( Fig. 3B-D). These results proved the high reproducibility and reliability of transcriptome profiling performed in the present study.

Functional enrichment analysis of neighboring target genes of differentially expressed lncR-NAs.
Recent studies suggested that lncRNAs may act in cis-regulation and affect the gene expression of their chromosomal neighborhood 10 kb upstream and downstream 36,37 . In this study, 49 significantly differentially expressed lncRNAs in PMM vs. SART were transcribed close to (< 10 kb) 70 mRNAs. Gene Ontology (GO) analysis of the cis lncRNA target genes was performed to explore their functions (GO, https:// www. omics hare. com/ tools/ home/ report/ goenr ich. html) (Fig. 4A). There were 178 extremely significantly enriched GO terms (p < 0.05) in PMM vs. SART. Many of the significantly enriched biological processes were associated with muscle cell differentiation and myofiber regulation, such as striated muscle cell differentiation, regulation of cell proliferation, regulation of muscle cell differentiation, myoblast differentiation, and regulation of myoblast differentiation. Interestingly, T-box 3 (TBX3) was annotated in many myoblast differentiation-related GO terms. T-box3 (TBX3) could regulate the cell cycle, inhibit cell apoptosis and promote proliferation 38,39 and was regulated by XR_003077811.1 in this study, which suggested that lncRNAs played important roles in regulating myofiber type.
To further understand the DE-lncRNA regulatory roles in myofiber type, Kyoto Encyclopedia of Genes and Genomes (KEGG, https:// www. omics hare. com/ tools/ home/ report/ koenr ich. html) pathway analysis was performed according to the cis-target genes (Fig. 4B). The MAPK signaling pathway was identified as a significantly enriched pathway for oxidative and glycolytic myofibers. Additionally, the MAPK family plays crucial roles in complex cellular processes, such as proliferation, differentiation, and development, by regulating the cell cycle and other cell proliferation-associated proteins 40,41 . Interestingly, only CALM2 was found to be differentially expressed in the MAPK signaling pathway and was regulated by XR_001465741.2. CALM2, as one subunit of CaM 42 , is a key redox-sensitive modulator of muscle physiology. Therefore, XR_001465741.2 was speculated to play a key role in myofiber growth. Functional enrichment analysis of host genes targeted by differentially expressed circR-NAs. Many studies have shown that circRNA expression can promote the transcription of host genes 46-48 . Therefore, to explore the mechanism of DE-circRNAs in myofibers, the host genes of DE-circRNA were used to perform a functional enrichment analysis (Fig. 5). There were 112 significantly enriched GO terms (p < 0.05) identified in the PMM and SART. The significantly overrepresented gene ontology (GO) terms were related to positive regulation of muscle hypertrophy, mitosis, positive regulation of developmental growth, and mononuclear cell proliferation (Fig. 5A). KEGG pathways, calcium signaling pathway, AMPK signaling pathway, FoxO signaling pathway, p53 signaling pathway, cellular senescence, etc., were closely related to myofiber development, resulting in the transformation of myofiber types (Fig. 5B). Ca 2+ signaling pathways were implicated at every step of myogenesis 49 , which contributed to the transformation of glycolytic myofibers to oxidative myofiber type 50 . AMPK has three subunits: α, β and γ; α is the catalytic subunit; and β and γ are the regulatory subunits. The distribution of the three subunits was related to muscle fiber type, and the β2 subunit was highly expressed in skeletal muscle and mainly distributed in glycolytic muscle 51 . The AMPK signaling pathway regulates the differentiation directions of myoblasts and changes the types of myofiber 52 . The FoxO signaling pathway is regulated www.nature.com/scientificreports/ by a variety of phosphorylated kinases and plays an important role in the proliferation and differentiation of myoblasts and the transformation of myofiber type 53 . p53-related genes promoted the cell cycle by upregulating p21 and enhancing muscle differentiation in MSTN knockout QM7 cells 54 . In this study, the host genes CAB39L, FOXO3, BMPR1B, SESN1, UBE2D1, PPP1CB, XPO1, NFATC1, BMPR1B, and PPP3CA of novel_circ_001395, novel_circ_003490, novel_circ_004266, novel_circ_003485, novel_circ_005184, novel_circ_003035, novel_ circ_002810, novel_circ_002121, novel_circ_004266, and novel_circ_004282 were involved in these pathways. These results suggest that these circRNAs might play important roles in regulating myofiber type. Interestingly, differentially expressed novel_circ_004282 and novel_circ_002121 were transcribed from the PPP3CA and NFATC1 genes, respectively, which are CaN substrates involved in maintaining and regulating muscle phenotypes, promoting the growth of different types of myofibers, and maintaining muscle mass in the skeletal muscle of mammals after birth [55][56][57] . The CaN-NFAT pathway was the best candidate for an activity-dependent signaling pathway responsible for the maintenance of the oxidative gene program in adult oxidative muscles and the induction of this program in regenerating oxidative muscle 58 . These results suggested that novel_circ_004282 and novel_circ_002121 might directly regulate fiber type heterogeneity in muscles.  Over the past decades, QKI has important functions in neural progenitors, myelin formation, smooth muscle differentiation, and monocyte to macrophage differentiation [59][60][61] . QKI could promote the myoblast differentiation of C2C12 cells in mice 31 . HSPA2 is abundantly expressed in skeletal muscle 62 and is induced by exercise-associated oxidative stress 63 . It prevents apoptosis by interacting with apoptosis-inducing factors (AIFs) 64 . MYBPC1 is known to encode the oxidative myofiber isoform of the major myosin-binding proteins in muscles [65][66][67] and acts as an adaptor to connect the ATP consumer (myosin) and the regenerator (MM-CK) 68 . The higher expression level of the MYBPC1 gene could result in more IMF deposition in skeletal muscle by controlling the energy www.nature.com/scientificreports/ metabolism and homeostasis of oxidative myofiber 69 . MYBPC1 could contribute to better meat quality of red muscle and was more highly expressed in oxidative myofibers than in glycolytic myofiber 70 . PPARGC1A is a coactivator of transcription involved in several aspects of skeletal muscle physiology, such as mitochondrial biogenesis, glucose utilization, fatty acid oxidation, thermogenesis, gluconeogenesis and insulin signaling 71 .

Validation of differentially expressed lncRNAs and circRNAs by qRT-PCR.
To validate the differential expression results of lncRNAs and circRNAs, the relative expression of three randomly selected lncRNAs and three randomly selected circRNAs was quantified by qRT-PCR. In Fig. 9A,B, all selected DE-lncRNAs and DE-circRNAs showed concordant expression patterns between the RNA-seq and qRT-PCR results.

Materials and methods
Ethics statement. All

Animal materials, tissue collection. A total of 15 female Qingyuan partridge chickens, provided by
Guangdong Tiannong Food Ltd, Guangdong, China, with similar body weights were euthanized by stunning followed by exsanguination at 140 days of age (marketing age). Pectoralis major and sartorius major on the left side were sampled immediately after slaughter, one part of the samples was stored in liquid nitrogen for frozen sectioning, and the other part was stored in liquid nitrogen for RNA extraction. Pectoralis major and sartorius major on the right side were sampled and stored at 4 ℃ for meat color measurement.

Meat color measurement.
After 24 h of storage, meat color was measured by a color difference meter (Shenzhen 3NH Technology CO., Ltd, Guangdong, China). The average of triplicate measurements was recorded, and the results were expressed as lightness (L*) and redness (a*) values.
Frozen section analysis. A 1 × 1 cm 2 section in the middle of the right SART and PMM was selected.
Measurement of the myofiber characteristics, including density, cross-sectional area, average diameter, and myofiber ratios, was carried out using ATPase staining. Myosin ATPase staining was used to identify myofiber type and to measure myofiber size 83,84 . The lightly dyed fibers were type I, the darker dyed fibers were type IIB, and the middle of the dyeing was type IIA. Oxidative myofibers included type I and IIA myofibers, and glycolytic myofibers were type IIB myofibers.
RNA extraction. Total RNA was extracted by using TRIzol reagent (Invitrogen, CA, USA) following the manufacturer's protocol. The RNA quantity of each sample was examined using a NanoDrop ND-2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) at 260/280 nm (ratio > 2.0). The integrity of total RNA was analyzed with the Agilent Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent Technologies) with RNA integrity number (RIN) ≥ 7.
Transcriptome library construction and sequencing. A total of eight cDNA libraries were constructed with four PMMs and four SART muscle tissues. A total of 3 µg RNA per sample was used as input material for RNA sample preparation. After the mRNAs and noncoding RNAs were enriched by removing rRNAs from the total RNA, the enriched RNAs were fragmented into short fragments and reverse transcribed into cDNAs. Buffer, dNTPs, RNase H, and DNA polymerase I were added to synthesize the second-strand cDNA. The resulting double-stranded cDNAs were ligated to adaptors after being end-repaired and A-tailed. Then, uracil-Nglycosylase (UNG) was used to digest the second-strand cDNAs. The digested products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced by Gene Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq 4000 85 .
Identification of lncRNAs and circRNA. The Illumina sequencing raw reads were obtained by removing adapter sequences, reads with poly-N, and low-quality reads, in which the number of bases with a quality value Q ≤ 20 was > 50%. All downstream analyses were based on high-quality clean data. The clean reads of each sample were mapped to the reference genome (Galgal 6.0) using TopHat (version 2.  www.nature.com/scientificreports/ Next, new transcripts were screened based on the position of the assembled transcript on the reference genome and the following screening criteria: transcript length ≥ 200 bp and the number of exons ≥ 2 (plant sample ≥ 1) to obtain known and new transcripts from the sample. Next, the Coding-Non-Coding Index (CNCI) 89 and Coding Potential Calculator (CPC) 90 were used to remove potential protein-coding transcripts. In this study, the resulting transcripts with no protein-coding potential in the two software analyses resulted in the lncRNA dataset.
In this study, find_circ software was used to identify circRNAs. Briefly, the unmapped back-spliced junction reads (default 20 bp) were used to extend the anchor sequences by find_circ software with default parameters 91 . In addition, identified circRNAs that were expressed in at least two samples were used for further analysis.  Functional enrichment analysis. Differentially expressed lncRNAs were selected for cis-target gene predictions and were integrated with differentially expressed gene data to improve the veracity of target prediction. In the present study, DE-mRNAs located ∼10 kb upstream and downstream of DE-lncRNAs were classified as cis-acting target genes; then, their functional roles were predicted as follows. The host genes of DE-circRNA were used to perform a functional enrichment analysis.
All DE-mRNAs and target genes of DE-lncRNAs and DE-circRNAs were annotated and classified by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis with OmicShare tools (HTTP ://www. omics hare.com/tools/) [43][44][45] . The results with a p value < 0.05 were considered significantly enriched. DE-mRNA and DE-lncRNAs were used to construct lncRNA-mRNA interaction networks by Cytoscape 3.5.1.
Construction of lncRNA-miRNA-mRNA network. The lncRNA-miRNA-mRNA network was constructed based on lncRNA-miRNA-mRNA theory as follows: (1) Expression correlation between mRNA-miRNA or lncRNA-miRNA was evaluated using the Spearman rank correlation coefficient (SCC). Pairs with SCC < − 0.7 were selected as coexpressed negative lncRNA-miRNA pairs or mRNA-miRNA pairs, both mRNA and lncRNA were miRNA target genes, and all RNAs were differentially expressed. (2) The expression correlation between lncRNAs and mRNAs was evaluated using the Pearson correlation coefficient (PCC). Pairs with PCC > 0.9 were selected as coexpressed lncRNA-mRNA pairs, and both mRNA and lncRNA in this pair were targeted and coexpressed negatively with a common miRNA. (3) A hypergeometric cumulative distribution function test was used to test whether the common miRNA sponges between the two genes were significant. As a result, only the gene pairs with a p value less than 0.05 were selected [93][94][95] .
Validation by qRT-PCR. To validate the differential expression results from sequencing, three lncRNAs and three circRNAs were selected for qRT-PCR. Total RNA for sequencing was reverse transcribed into cDNA using the PrimeScript RT reagent kit (TaKaRa, Dalian, China). Then, qRT-PCR was conducted using a KAPA SYBR Fast universal qPCR kit (Kapa Biosystems, USA). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) genes were used as the internal reference. All primers are shown in Table 3.

Statistical analysis.
Comparisons of the two myofiber types were analyzed using an independent sample T-test procedure in SPSS (Version 20.0, SPSS, Inc., Chicago, IL, USA). Differences between PMM and SART muscle samples were considered statistically significant at p < 0.05.

Data availability
The datasets supporting the results presented here are available in the Sequence Read Archive (SRA) repository under accession number PRJNA578179. Table 3. Primers for quantitative real-time PCR.