The methods used in this study were performed in accordance with the guidelines of Good Experimental Practices adopted by the Institute of Animal Science (Qingdao Agricultural University, Qingdao, China). All surgical procedures involving cattle were performed according to the approved protocols of the Biological Studies Animal Care and Use Committee, Shandong Province, China (Cao et al., 2014).

Shandong black cattle(three)and Luxi cattle(three)were used in this study. Cattle were fed three times and received approximately equal amounts of green coarse feed and concentrate feed everyday according to standard NY5127—2002 pollution-free feeding management of beef cattle. They were raised in the same environment. At the age of 18 months, three healthy male beef cattle in each group were randomly selected for this study. The selected cattle had no scratches, scars, and scabs on their bodies, with no fat deposits in the internal organs or abdomen. No disease was found during examination, and all physiological and biochemical indexes were normal. Longissimus dorsi muscle samples were collected and immediately frozen in liquid nitrogen for RNA extraction. The cattle used in the experiment were euthanized as follows: first, a 1–3% sodium pentobarbital solution was prepared with physiological saline, and then intravenously injected. The injection dose is 90–135 mg/kg.

In order to better observe the histological morphology of muscle, we performed HE staining and fast/slow muscle staining. Paraffin sections were prepared from muscle tissue fixed with 4% paraformaldehyde. The HE staining protocol was performed as described previously (Liu et al., 2020). Briefly, dewaxing, covering with water, Hematoxylin staining, washing with water, 5% acetic acid differentiation, eosin staining, dehydration, natural drying, sealing, and image acquisition were performed.

Tissue sections were placed in a box filled with EDTA antigen repair buffer (Purchased from Shanghai Beyotime Biotechnology Co., Ltd) (ph 8.0) for antigen repair. After natural cooling, the slides were washed in PBS (pH7.4) with three ashes of 5 minutes each. BSA was added for blocking, and then the antibody was added, followed by DAPI to stain the nucleus. An autofluorescence quenching agent was added to the slices for 5 min and the samples were then washed with running water for 10 min. After natural drying, the film was sealed with an anti-fluorescence quenching sealing agent. Finally, the slices were observed under the fluorescence microscope and images were collected.

Image-Pro Plus software was used to count the images and measure the surface area. SPSS software was used for statistical analysis to determine significant differences.

Total RNA samples were isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a NanoPhotometer^® spectrophotometer (IMPLEN, Los Angeles, CA, United States). RNA concentration was measured using a Qubit^® RNA Assay Kit in a Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States). RNA integrity was assessed using a RNA Nano 6000 Assay Kit in a Bioanalyzer 2,100 system (Agilent Technologies, Santa Clara, CA, United States). Only samples with RNA Integrity Number (RIN) scores >8 were used for sequencing, which is different from the previous studies (Liu et al., 2020). The different RIN values are due to the different requirements for RNA quality when constructing different data sets. A total of 3 μg RNA per sample was used as the input material for RNA library preparation. This study used Illumina hiseq xten sequencing platform.

The original sequencing data were processed by Fastp software (v0.23.2) (Chen et al., 2018) with parameters “-Q 20 -P90” with disjointing sequence and low-quality sequence. Clean data were obtained by removing reads containing adapters, reads containing over 10% of poly(N), and low-quality reads (>50% of the bases had Phred quality scores ≤10) from the raw data. All downstream analyses were based on the high quality clean data. The Bos taurus genome reference genome and gene model annotation files were downloaded from the NCBI database (CHIR_1.0, NCBI) (Mahmoudi et al., 2020). An index of the reference genome was built using Bowtie v2.0.6 (Cai et al., 2015; Cai et al., 2018) and paired-end clean reads were aligned to the reference genome using HISAT v2-2.1.0 (Yuan et al., 2019). The mapped reads from each library were assembled with Cufflinks v2.2.1 (Andersson, 2009), using Cufflinks with ‘min-frags-per-transfrag = 0’ and’–library-typefr-firststrand’, and other parameters set as default.

Stringtie (Pertea et al., 2015; Kovaka et al., 2019) was used to sort reads into different classes and then generate a map for each class. Based on the length of LncRNA, and the characteristics of non-protein coding sequences, we established strict screening conditions to screen LncRNA as follows:

(1) Transcripts equal to or longer than 200 bp in length, and containing two or more exons;

To further assess if the identified transcripts are LncRNAs, a variety of coding potential analysis software was integrated, including CNCI analysis, CPC analysis, Pfam protein domain analysis, and CPAT analysis (only for animals) (Langfelder et al., 2008; Jing et al., 2010; Luo et al., 2017; Ito et al., 2018). The transcripts identified as non-coding by several methods were the final potential LncRNA dataset.

Comparison of raw counts data (Wang et al., 2010) for different genes is a very effective tool for quantitative estimation of gene expression based on RNA-Seq data. This method can eliminate the influence of gene length and sequencing quantity on the determination of gene expression levels. The calculated gene expression values allow direct comparison of gene expression differences among different varieties. Here, we used DEseq (Anders and Huber, 2010) to analyze the differential expression between the treatment group with the reference group, and selected |Log2ratio| ≥2 and q < 0.05 genes as indicative of significant differential expression. The numbers of up-regulated and down-regulated genes were also obtained.

LncRNA can act on target genes by Cis or Trans. Cis effect refers to LncRNA acting on adjacent target genes (Cabili et al., 2011). To predict the cis-regulatory genes of LncRNA, we screened the 50 kb of sequence upstream and downstream of LncRNA and looked for potential target genes. The coexpression relationship between LncRNA and mRNA was described by Trans. The basic principle of transaction is that the function of LncRNA has nothing to do with the position of the coding gene,but it is related to the protein-coding genes it coexpressed. According to the correlation coefficient between LncRNA and mRNA expression (correlation coefficient cor ≥ 0.9); WGCNA(Weighted correlation network analysis) was used to predict the target genes of lncRNA for cluster analysis and functional enrichment analysis of LncRNA target genes. p-value < 0.05 was set as the significance threshold.

GO (gene ontology) (Young et al., 2010) can be used for enrichment analysis of target genes with differential expression and Goseq (Young et al., 2010) was used to analyze the target genes of differentially expressed LncRNA. A value of p < 0.05 is considered as significant enrichment of differentially expressed genes.

The KEGG database (Kanehisa and Goto, 2000) describes advanced functions and utility of biological systems such as cells, organisms, and ecosystems (http://www.genome.jp/kegg/). We used KOBAS V3.0 software (Bu et al., 2021) to detect the enrichment of LncRNA target genes differentially expressed in the KEGG pathway.

As competing endogenous RNAs (ceRNA), LncRNA can act as a sponge to adsorb miRNA and affect mRNA expression. Therefore, to further understand the function of LncRNA in the growth and development of skeletal muscle, we used TargetScan (http://www.targetscan.org/vert_71/) miRanda ((http://www.microrna.org/microrna/home.do) and PITA qtar (https://genie.weizmann.ac.il/pubs/mir07/mir07_dyn_data.html) software to predict the differentially expressed miRNAs(Sequencing has been completed) adsorbed by LncRNA and downstream-regulated mRNA. Online miRNA binding site prediction software (RNA22: http://cbcsrv.watson.ibm.com/rna22.html and RNAhybird ()) predicted potential interaction of miRNA with LncRNA. We then constructed an LncRNA-miRNA-mRNA interaction network map using Cytoscape version 3.5.1 (Shannon et al., 2003).

Cells were seeded in 96-well plates at a density of 5 × 103 cells (HEK-293T) per well, 24 h before transfection. The cells were co-transfected with a mixture of 50 ng Firefly luciferase (FL) reporter vectors, 5 ng Renilla luciferase (RL) reporter vectors (pRL-TK), and miRNA mimics at the indicated concentration. The miRNA mimics were obtained from Life Technologies. After 48 h, the luciferase activity was measured with a dual luciferase reporter assay system using the psiCHECK-2 vector (Promega, Madison, WI). The LOC104975788 and LOC536229 (Pax7) sequences were separately cloned into the reporter gene vector (psiCHECK-2) to synthesize the predicted miRNA mimics and control. The potential binding target of each miRNA was cloned into the 3′UTR region of r-luciferase (hrluc), and then co-transfected with the miRNA to determine the activity of R-Luciferase. F-Luciferase (hluc +) was used as an internal reference to correct for differences in transfection efficiency between different samples. The miRNA mimics and psicheck-LOC104975788 or psicheck-LOC536229 (Pax7) were co-transfected into 293T cells. The expression level of reporter genes was detected using a multi-functional enzyme labeling instrument, and the miRNAs that exhibited down-regulated reporter gene expression were further screened.

The reaction system (20 μl) for the RT-PCR reaction consisted of the following: 1 μl of template cDNA, 10 μl each of the upstream and downstream primers, and 5 ml (5 × 1 ml vials) of RNase-free water. The thermal cycling procedure was as follows: 94°C for 10 min, 94°C for 30 s, 60°C for 30 s, and 72°C for 40 s, with 40 cycles12. The expression of GAPDH was calculated by the 2^-△△CT method. Primers used for qRT-PCR as shown in Supplementary Table S1.

First, we obtained samples from Shandong black cattle and Luxi cattle and performed fluorescence staining sections of fast and slow muscle fibers (Figure 1A) and HE staining (Figure 1B). The results showed significant differences in muscle fibers of longissimus dorsi muscles between Shandong black cattle and Luxi cattle. IPP(Image-Pro Plus) software analysis showed that the average area of muscle fibres of Shandong black cattle was significantly larger than that of Luxi cattle (p < 0.05), with significant differences in the muscle fibre density and the ratio of fast-twitch fibers to slow-twitch fibers and slow-twitch fibers to muscle fiber area (p < 0.05), but not in other muscle fibre properties (p > 0.05) (Table 1).

RNA was extracted from longissimus dorsi muscle samples of Shandong black cattle and Luxi cattle and sequenced. Totals of 105345396 and 126735736 reads were obtained from Shandong black cattle and Luxi cattle, respectively, with 101145946 and 122463284 mapped reads. In the data for Luxi cattle longissimus dorsi muscle, the percentages of exons, introns, and genes were 43.26%, 41.34%, and 15.4%, respectively, with 38.16%, 48.3%, and 13.54%, respectively, for Shandong black cattle longissimus dorsi muscle tissue, as shown in Figure 2.

A total of 20104 transcripts were obtained from the two libraries, including 18859 and 18976 mRNAs for Luxi cattle and Shandong black cattle, respectively. According to the screening criteria of |FC|≧2 and q < 0.05, 1,415 differentially expressed mRNA (Supplementary Table S2) were found, with 939 mRNAs significantly up-regulated and 476 down-regulated in Luxi cattle. The 92 mRNAs were specifically expressed in cattle, with 30 and 62 mRNAs specifically expressed in Shandong black cattle and Luxi Cattle, respectively (Supplementary Table S2). The top 40 genes with significant differential expression in the two breeds in cattle are listed in Supplementary Table S3, and suggest genes that could be used for the functional comparison and analysis of differences among varieties.

Gene Ontology (GO) is an international standardized classification system of gene function that includes three categories: biological process (BP), molecular function (MF), and cell composition (CC). We performed GO enrichment analysis on 1,415 differentially expressed genes (Figures 3A,B), and the results showed 2,662 items (p < 0.05), with 388 were enriched in cell composition (CC), 361 in molecular function (MF), and 1913 in biological process (BP) in up-regulated genes (Shandong black cattle as a reference group). The enrichment results of down-regulated genes showed 1883 items in total (p < 0.05), with 251 enriched in cell composition (CC), 247 in molecular function (MF), and 1,385 in biological process (BP). In the muscle regulation system, muscle contraction is significantly enriched in GO terms, indicating that some differential genes may participate in muscle development and function and other life activities. Muscle-related GO terms are shown in Table 2.

Different genes coordinate with each other to perform their biological functions in organisms. We next performed KEGG pathway analysis of the differentially expressed genes to identify the related signal transduction and metabolic pathways (Figures 3C,D). The results showed enrichment of 939 up-regulated genes in 243 pathways (p < 0.05), and enrichment of 476 down-regulated genes in 239 pathways, including metabolism, disease, cell communication and endocrine regulation (Supplementary Table S4). The enriched pathways (p < 0.05) included some signaling pathways related to skeletal muscle growth, development, and function (Table 3), such as the Wnt signaling pathway, the calcium signal transduction pathway, gap junction, the Hippo signaling pathway, the AMPK signaling pathway, and the thyroid hormone synthesis pathway.

A total of 8,427 LncRNAs were found in the small RNA Library of longissimus dorsi muscle of Shandong black cattle and Luxi cattle. Of these, 3,498 LncRNAs were annotated in non-coding regions and the remaining 4,929 LncRNAs were in regions that are not annotated. According to the screening criteria of |FC|≧2 and q < 0.05, 480 differentially expressed LncRNAs were found, with 245 and 235 LncRNAs significantly up-regulated and down-regulated in Luxi Cattle compared to the levels in Shandong black cattle (Supplementary Table S5). A total of 92 LncRNAs were only expressed in one cattle breed, with 30 specific LncRNAs in Shandong black cattle and 62 specific LncRNAs in Luxi cattle (Supplementary Table S5). These breed-specific longissimus dorsi muscle LncRNAs included eight annotated LncRNAs in Shandong black cattle and 13 annotated LncRNAs in Luxi cattle, as shown in Table 4. The LncRNAs specifically expressed in Shandong black cattle and Luxi cattle, the first 40 unexplained LncRNAs and annotated LncRNAs are listed in Supplementary Table S6.

We constructed a Circos diagram to visually display the distribution of differential LncRNAs and differential mRNAs (Figure 4A). The results showed that most differentially expressed LncRNAs mapped to chromosome 9, with others mapping to chromosomes 3, 1, 11, and 7. Pearson correlation coefficient analysis of differentially expressed LncRNA and differential mRNA genes (Figure 4B) showed 43844 pairs of co-expressed LncRNAs and differential mRNA genes with correlation coefficient >0.9 (Supplementary Table S7). Based on genomic location, 387 differentially expressed LncRNAs were selected for further study.

LncRNAs can regulate the expression of adjacent mRNAs. We analyzed the LncRNAs and protein-coding genes by analyzing mRNAs within 50 kb of LncRNAs. The results showed that 387 of the differentially expressed LncRNAs targeted 1,164 genes. In these predictions, one LncRNA targeted multiple mRNAs and many LncRNAs targeted the same mRNA (Supplementary Table S8). Of these, LOC112448071, LOC112444635, LOC112445963, LOC104975064, LOC101903261, LOC535280, and LOC521224 can target genes related to muscle development, including MYORG, Wnt4, PAK1, and ADCY7.

The potential target genes of LncRNA were next subjected to GO and KEGG analysis (Figure 5). We found 130 significantly enriched GO items (p < 0.05), with five related to the regulation of muscle development (Table 5). The 40 most enriched terms include multicellular organization development, single multicellular organization process, and multicellular organic process (Supplementary Table S9). The GO terms related to muscle development are listed in Table 5. KEGG analysis showed enrichment of 1,164 potential target genes in 299 pathways, with the 40 most significantly enriched pathways listed in Supplementary Table S9. Some of these are related to muscle development (Table 6), including the calcium signaling pathway, the AMPK signaling pathway, the cGMP-PKG signaling pathway, and the PPAR signaling pathway.

We used Targetscan (Agarwal et al., 2015), Miranda (Enright et al., 2003), and PITA qtar software (Li et al., 2018) to predict the target gene relationships between mRNA and miRNA, and used Miranda and PITA qtar to predict the target gene relationships between miRNA and LncRNA. LncRNA can function as ceRNA to regulate the expression of downstream genes at the post-transcriptional level by binding miRNAs. We selected several miRNAs related to muscle development (including miR-1, miR-133, miR-206, miR-181, miR-21, and miR-23a) to predict and screen downstream-regulated mRNA using Targetscan and Miranda software. We then constructed an LncRNA-miRNA-mRNA interaction network diagram using Cytoscape version 3.5.1 (Figure 6). LncRNAs that may be involved in the development of beef cattle muscle are summarized in Table 7.

We next looked for potential miRNA binding sites of LOC104975788 and found potential miR-133a and miR-103 binding sites. Online miRNA binding site prediction software (RNA22 and RNAhybird) predicted potential interaction of miR-133a with LOC104975788 and LOC536229(Pax7). The results showed that both LOC104975788 and LOC536229 (Pax7) contained miR-133a binding sites.

Nine differentially expressed LncRNAs and five differentially expressed mRNAs were verified by qRT-PCR. The qRT-PCR results (Figures 7D,E) confirmed differential expression of these LncRNAs and mRNAs.