The experiment was executed in strict accordance with the guidelines for the care and use of experimental animals at Jilin University. All the experiments were approved by the Institutional Animal Care and Use Committee of Jilin University (license number: 201809041).

Three Kazakh cattle and three Xinjiang brown cattle were provided by the Xinjiang Yili Yixin Cattle and Sheep Breeding Cooperative. The beef cattle were slaughtered in accordance with the procedures of the slaughterhouse, and the longissimus muscle was collected at the slaughter line. Whole specimens were immediately snap-frozen in liquid nitrogen and stored at −80°C.

RNA was extracted from six pieces of longissimus muscle (three pieces from Kazakh cattle and three pieces from Xinjiang brown cattle) using the TRIzol reagent (Tiangen, Beijing, China) in accordance with the manufacturer’s guidelines. The RNA quantity was determined with a NanoDrop rdt-2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). After testing, the total RNA samples from the longissimus muscle tissue of three Kazakh cattle and three Xinjiang brown cattle were reverse transcribed and pooled for the construction of Kazakh cattle and Xinjiang brown cattle cDNA libraries. To construct these two libraries, we first extracted rRNA from the samples using an Epicenter Ribo-Zero kit and then added fragmentation buffer to the rRNA-depleted RNA. In addition, we synthesized first-strand cDNA using random hexamers with the rRNA-depleted RNA as the template, and buffer, dATP, dUTP, dCTP, dGTP, RNase H and DNA polymerase 1 were then added to the first-strand cDNA. AMPure XP beads were used to purify the cDNA. The purified double-stranded cDNA was subjected to end repair (poly-A addition and end ligation), and individual segment sizes were selected using AMPure XP beads. The U-containing cDNA libraries were obtained via PCR enrichment. In addition, the quality of the cDNA libraries was evaluated using the Qubit2.0 and Agilent 2100 instruments and quantitative reverse transcription polymerase chain reaction (RT-qPCR). The libraries were sequenced by BioMarker Technologies (Beijing, China) using an Illumina HiSeq X Ten platform.

The transcripts obtained from the longissimus muscle of the Kazakh cattle and Xinjiang brown cattle libraries were assembled, and several steps were used for the prediction of lncRNAs. Four computational tools, coding potential calculator (CPC), coding-non-coding index (CNCI), Pfam, and coding potential assessment tool (CPAT), were used to screen for protein-coding ability. Two-exon transcripts with a length greater than 200 nt were reserved as lncRNA candidates and subjected to further CPC/CNCI/Pfam/CPAT screening to ensure that each candidate transcript was a lncRNA. CPC2 is a tool for calculating the protein-coding potential based on sequence alignment. By aligning a transcript with known protein databases, CPC2 evaluates the coding potential of the transcript based on the biological sequence characteristics of each coding cassette of the transcript. CNCI allows distinguishing coding and noncoding transcripts based on adjacent nucleotide triplet features, and a score less than 0 indicates that a noncoding transcript. CPAT determines the coding and noncoding ability of a transcript by constructing a logistic regression model based on the ORF length and coverage and calculating the Fickett and Hexamer scores. The Pfam database is the most comprehensive classification system for protein domain annotations. Pfam divides the protein domain into different protein families and establishes an HMM statistical model of the amino acid sequence of each family through protein sequence alignment. An aligned transcript is defined as a transcript with a certain protein domain and is thus considered a coding transcript, and an unmatched transcript is a potential lncRNA. We downloaded the Bos taurus reference genome from the Ensembl genome browser (http://www.ensembl.org/Bos_taurus/Info/Index), and the lncRNA transcripts were classified through comparisons as lncRNAs, intronic lncRNAs, antisense lncRNAs or sense lncRNAs.

The expression of lncRNAs was analyzed using StringTie (1.3.1) in terms of the fragments per kilobase of transcript per million mapped reads values of the encoding genes. This program uses a scaling normalization factor for each sequence library. The differential expression of the transcripts in the two samples was analyzed using the DESeq R software package. The resulting false discovery rate (FDR) was controlled by the Benjamini and Hochberg’s approach, and an FDR<0.05 and a log2 (fold change)>1.5 were used as the criteria for significantly differentially expressed genes.

We found 100 kb of coding genes upstream and downstream of all the detected lncRNAs and predicted their effects. The names of the adjacent genes were subjected to gene ontology (GO) analysis using R software. KOBAS software [18] was utilized to measure the expression of different genes in Kyoto encyclopedia of genes and genomes (KEGG) pathways. The GO terms and KEGG pathways with p<0.05 were deemed to be significantly enriched.

The primers used for quantitative polymerase chain (PCR) were designed with Primer3 (http://bioinfo.UT.Ee/primer3-0.4.0/) and were synthesized by Changchun Kumei Bio (Supplementary Table S5). Total RNA was transformed into cDNA using an InRcute lncRNA First-Strand cDNA Synthesis Kit (with gDNase) and a FastQuant RT kit (with gDNase) (Tiangen, China), and RT-qPCR was performed using a Mastercycler ep realplex 2 system (Eppendorf, Germany). An INRCARE qPCR Test Kit (SYBR Green) and surrealistic premix + (SYBR Green) were used in accordance with the manufacturer’s instructions. Glyceraldehyde-3-phosphate dehydrogenase was used as an internal reference in the lncRNA analysis.

The results are shown as the means±standard deviations from three independent experiments and were analyzed using SPSS 19.0 [19,20]. The significance of the differences was examined by one-way analysis of variance [21], and p<0.05 was considered to indicate statistical significance.

To investigate the lncRNA expression profiles of the longissimus muscle, two cDNA libraries were constructed from longissimus muscle tissues of Kazakh cattle and Xinjiang brown cattle using the Illumina HiSeq platform. We performed paired-end sequencing and the read length was 150 bp. In this study, a total of 115.43 Gb of clean lncRNA sequence data was acquired from six samples. The clean data from each sample reached 18.36 Gb, and the Q30 base percentage was greater than 92.89%. The clean reads from each sample were aligned with the designated reference genome, and the efficiency of the alignment ranged from 93.35% to 94.79%. Circos was then employed to visualize the putative locations of the lncRNAs in the cattle genome. LncRNA transcripts were distributed in almost all chromosomes (Figure 1A). The differential expression of lncRNA and mRNA sequences on the chromosomes was analyzed. The outer loop in the ring diagram shows the chromosome of the species reference genome, the middle loop shows the distribution of differentially expressed mRNAs on the chromosome, and the inner loop shows the distribution of differentially expressed lncRNAs on the chromosome. In the Figure1A, the upregulated and downregulated genes are shown in red and green, respectively, and the upregulated and downregulated lncRNAs are presented in yellow and blue, respectively. A cross-analysis using the Pfam, CNCI, CPAT, and CPC tools was performed to obtain 11,964 lncRNAs (Figure 1B), including 7,856 lncRNAs (65.7%), 1,302 antisense lncRNAs (10.9%), 2,179 intronic lncRNAs (18.2%), and 627 sense lncRNAs (5.2%) (Figure 1C). We detected 279,170 protein-coding transcripts, and the lengths of most protein-coding and lncRNA transcripts were within the range of 200 to 800 bp (Figure 1D). Most of the protein-coding transcripts comprised one exon, whereas most of the lncRNA transcripts comprised two to three exons, and the proportion of lncRNA transcripts with two exons accounted for 73.93% of all lncRNAs (Figure 1E).

The DESeq R package was used to create volcano and MA plots of the differentially expressed transcripts. The volcano plot shows the symmetrical distribution of the main transcribed genes (Figure 2A). All abnormally expressed transcripts were screened to identify transcripts with a fold change ≥1.5. A total of 182 lncRNAs were differentially expressed in the longissimus muscle between Kazakh cattle and Xinjiang brown cattle (Supplementary Table S1), and 102 and 80 of these differentially expressed lncRNAs were upregulated and downregulated, respectively. A systematic cluster analysis was performed to elucidate the similarities and differences between the different libraries by examining the expression of differentially expressed lncRNAs. The DESeq R package was used for transcript expression analysis to detect lncRNAs that were differentially expressed between Kazakh cattle and Xinjiang brown cattle. Heatmaps were constructed to show the differentially expressed lncRNAs (Figure 2B). The results indicate the existence of particular lncRNAs that enable differentiation of longissimus dorsal tissues from Kazakh cattle and Xinjiang brown cattle.

We studied the potential roles of the lncRNAs by searching for protein-coding genes 100 kb upstream and downstream of all the detected lncRNAs to predict the underlying cis-regulatory targets of the lncRNAs. Overall, we found that some pathways associated with muscle metabolism and energy metabolism, such as the mitogen-activated protein kinase (MAPK) and phosphatidylinositol 3-kinase (PI3k)/Akt signaling pathways and signaling pathways related to spliceosomes, Ras, focal adhesion, and ribosome biogenesis in eukaryotes, were markedly enriched (Figure 3A and Supplementary Table S3). In addition, 60 GO terms were dramatically enriched (Supplementary Table S2), and these were mostly associated with single-organism processes (GO: 0044699), cellular processes (GO: 0009987), organelles (GO: 0043226), cell parts (GO: 0044464), biological regulation (GO:0065007), and other GO terms(Figure 3B).

The possible interactions among mRNAs and lncRNAs can be examined through lncRNA/mRNA coexpression networks. Hence, a coexpression network was constructed using the differentially expressed lncRNAs and mRNAs between Kazakh cattle and Xinjiang brown cattle. Because the target genes of the differentially expressed lncRNAs between Kazakh cattle and Xinjiang brown cattle were enriched in the PI3k/Akt signaling pathway, we focused on the mechanism through which lncRNAs regulate genes in this pathway (Figure 4). The network suggested that insulin-like growth factor 1 receptor (IGF1R), insulin receptor substrate 1 (IRS1), AKT serine/threonine kinase 1 (AKT1), phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1), mammalian target of rapamycin (mTOR), glycogen synthase kinase 3 alpha (GSK3A), glycogen synthase kinase 3b (GSK3B), forkhead box O1 (FOXO1), and nuclear factor kappa B subunit 1 (NFKB1) were correlated with 64 lncRNAs (Supplementary Table S4), which indicates that the lncRNAs and mRNAs regulate the PI3k/Akt signaling pathway during muscle development and differentiation.

To verify the RNA-seq results, eight genes were selected for validation by quantitative RT-PCR (RT-qPCR). Eight lncRNA transcripts with differential expression were selected at random, and their relative expression in longissimus muscle tissues from Kazakh cattle and Xinjiang brown cattle was verified (Figure 5A). The RT-qPCR results were highly consistent with the RNA-seq results. Based on sequence complementarity, we used the lncRNA Targets platform to establish an interaction network between the lncRNAs and IGF1R genes. The expression levels of all three lncRNAs (MSTRG.100266.1, MSTRG.103991.1, and MSTRG.119236.1) and IGF1R were measured by RT-qPCR (Figure 5B–5E). The MSTRG.100266.1, MSTRG.103991.1, and MSTRG.119236.1 expression results were exactly the opposite of the IGF1R expression results, which showed that MSTRG.100266.1, MSTRG.103991.1, and MSTRG.119236.1 might regulate muscle growth and development by targeting IGF1R.