Genome-Wide Identi ﬁ cation and Characterization of Long Non-Coding RNAs in Longissimus dorsi Skeletal Muscle of Shandong Black Cattle and Luxi Cattle

There is an increasing understanding of the possible regulatory role of long non-coding RNAs (LncRNA). Studies on livestock have mainly focused on the regulation of cell differentiation, fat synthesis, and embryonic development. However, there has been little study of skeletal muscle of domestic animals and the potential role of lncRNA. In this study, the transcriptome numbers of longissimus muscle of different beef cattle (Shandong black catle and Luxi catlle) were used to construct muscle related lncRNAs-miRNA-mRNA interaction network through bioinformatics analysis. This is helpful to clarify the molecular mechanism of bovine muscle development, and can be used to promote animal husbandry and improve animal husbandry production. According to the screening criteria of |FC| S 2 and q < 0.05, a total of 1,415 transcripts (of which 480 were LncRNAs) were differentially expressed ( q < 0.05) in the different breeds. Further, we found that the most differentially expressed LncRNAs were found on chromosome 9, in which the differentially expressed LncRNAs targeted 1,164 protein coding genes ( MYORG , Wnt4 , PAK1 , ADCY7,etc ) (upstream and downstream < 50 Kb). In addition, Pearson ’ s correlation coef ﬁ cients of co-expression levels indicated a potential trans regulatory relationship between the differentially expressed LncRNAs and 43844 mRNAs (r > 0.9). The identi ﬁ ed co-expressed mRNAs ( MYORG , Dll1 , EFNB2 , SOX6 , MYOCD , and MYLK 3) are related to the formation of muscle structure, and enriched in muscle system process, Interestingly, a candidate differential LncRNA ( LOC104975788 ) and a protein-coding gene ( Pax7 ) contain miR-133a binding sites and binding was con ﬁ rmed by luciferase reporter assay. LOC104975788 may combined miR-133a competitively with Pax7 , thus relieving the inhibitory effect of miR-133a on Pax7 to regulate skeletal muscle development. These results will provide the theoretical basis for further study of LncRNA regulation and activity in different cattle breeds. with LOC104975788 and LOC536229 (Pax7) target site (A) . Luciferase reporter assay showed signi ﬁ cantly lower luciferase activity of miR-133a co transfected with the psi-LOC112448162 vector compared to the control group ( p < 0.01). The luciferase activity of miR-133a co-transfected with psi-LOC536229 (Pax7) vector was also signi ﬁ cantly decreased relative to the control Q19 group ( p < 0.01) (B, C) . These results indicate that both LOC104975788 and LOC536229 (Pax7) contain miR-133a binding sites. Nine differentially expressed LncRNAs and ﬁ ve differentially expressed mRNAs were veri ﬁ ed by qRT-PCR (D, E) .

Interestingly, a candidate differential LncRNA (LOC104975788) and a protein-coding gene (Pax7) contain miR-133a binding sites and binding was confirmed by luciferase reporter assay. LOC104975788 may combined miR-133a competitively with Pax7, thus relieving the inhibitory effect of miR-133a on Pax7 to regulate skeletal muscle development. These results will provide the theoretical basis for further study of LncRNA regulation and activity in different cattle breeds.

INTRODUCTION
As an important economic trait that affects the production efficiency of beef cattle, meat production traits are a research focus in the field of beef cattle genetics and breeding. The analysis of the molecular regulation mechanisms of muscle growth can facilitate cattle breeding. Local cattle breeds in China typically present reduced growth, but increased meat quality when compared to imported cattle. These differences may be related to differences in skeletal muscle development among different breeds (Xu et al., 2018). The growth and development of beef cattle skeletal muscle can be affected by a variety of regulatory factors. Previous work mainly focused on contributions from DNA, mRNA, and miRNA, but regulatory effects of long non-coding RNA (LncRNA) on the growth and development of beef skeletal muscle remain poorly understood (Zhang et al., 2007;Dylan et al., 2008;Xu et al., 2020).
Previous studies reported contributions of LncRNA in the process of skeletal muscle proliferation and differentiation. Liu identified, mapped, and determined the global and skeletal muscle expression patterns of 7,188 LncRNAs, finding that these LncRNAs had similar open reading frames and expression levels as those of other mammalian LncRNAs (Sun et al., 2016;Yu-ying et al., 2017). Subsequently, Yue identified a highly expressed LncRNA-YYW in muscle tissue (Yue et al., 2017). Microarray analysis showed that LncRNA-YYW positively regulated the expression of growth hormone-1 and its downstream genes AKT1 and PIK3CD in bovine myoblasts, and promoted myoblast proliferation (Dylan et al., 2008;Billerey et al., 2014). Cai found that lnc-ORA interacted with IGF2BP2 to inhibit the PI3K/AKT signaling pathway, thereby inhibiting muscle production and inducing skeletal muscle atrophy (Cai et al., 2021). Studies had shown that LncMD could be used as a competitive endogenous RNA (ceRNA) to compete with IGF2 and bound miR-125b to weaken its inhibitory effect on IGF2, and promoted the differentiation of bovine skeletal muscle satellite cells (Sun et al., 2016). Linc-MD1 had been shown to be a competitive endogenous RNA for miR-133 and miR-135 targets during myoblast differentiation. These transcription factors activated muscle-specific gene expression to control the time of muscle differentiation (Marcella Davide et al., 2011). Lnc-31 maintained the proliferation of myoblasts and counteracts differentiation (Dacia et al., 2018). The active regulation of Rock1 translation by lnc-31 was of great significance to the control of myogenesis (Ballarino et al., 2015). These results suggest that LncRNA may regulate the growth and development of beef skeletal muscle and meat quality. Although LncRNAs lack coding capacity, many LncRNAs act in various biological processes, serving as an additional level of genomic regulation.
The first embryo transfer calf in China was obtained from a vitrified frozen somatic cell cloned embryo, which combined Japan black cattle, Luxi cattle, and Bohai black cattle. The new generation core breeding group of this new breed, Shandong black cattle, was established by the filial generation. The new generation core breeding group hybridized with Shandong black cattle, and the second generation group of Japan black cattle (3/4) was crossed with Luxi cattle (1/4) to produce the new generation (3/4) × Bohai black cattle (1/4). Ideal black cattle and Shandong black cattle were selected from the second-generation group for cross-cross fixation. In this stage, several good lines were crossbred with distant relatives, and then selected and the excellent lines were preserved. After four generations, this line was finally bred into a new Shandong black cattle variety matching line and used as a bull. In 2015, the Shandong black cattle variety was approved by the National Animal and Poultry Genetic Resources Committee as a new population and successfully established as a new cultivated variety in China (Lu, 2015;Liu et al., 2018).
In recent years, there has been growing research of LncRNA, but most work has mainly focused on the effects of neurologic diseases, tumors, embryonic development, and cell differentiation in human and mouse, with little research on domestic animals. The growth and development of beef cattle skeletal muscle are highly regulated processes. Although the contributions of DNA, mRNA, and miRNA to this regulation have been studied, the extent to which LncRNA regulates beef skeletal muscle growth and development remains poorly understood. To investigate the potential contribution of LncRNA, RNA samples were prepared from longissimus dorsi muscle tissues of Shandong black cattle (hybrid offspring) and Luxi Cattle (the first maternal generation). RNA-seq technology was used to identify LncRNA transcripts and their genomic characteristics, and LncRNAs differentially expressed in skeletal muscle tissues in Shandong black cattle and Luxi cattle were determined. Using these data, we constructed an interaction network of LncRNA-miRNA-mRNA and muscle development which can guide future breeding efforts. 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).

Animal and Tissue Preparation
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 FIGURE 1 | Fast/slow muscle fluorescence staining and paraffin section HE staining of muscle tissue. (A) Green represents slow muscle fiber, red represents fast muscle fiber, and blue represents the nucleus of the muscle cell. (B) Red represents a single muscle fibre, Blue represents a nucleus. 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.

Hematoxylin and Eosin Staining of Muscle Tissue and Fast/Slow Muscle Fiber Fluorescence Staining
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.

RNA Extraction, Library Construction, and Sequencing
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.

Transcriptome Assembly
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 The results showed that the genomic distribution of clean data of longissimus dorsi muscle tissue of Shandong black cattle and Luxi cattle was genome distribution; the abscissa was each sample, and the ordinate was the ratio of the sequence number of exon, intron and intergenic region between the reference genome.
Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 849399   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.

LncRNA Discovery
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 nonprotein 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; (2) Transcript read coverage of at least five reads; (3) No transcripts of known mRNA or other specific noncoding RNAs (rRNA, tRNA, snoRNA, or snRNA). This screening was based on gffcompare (http://ccb.jhu.edu/ software/stringtie/gff.shtml) using the same species annotation file; (4) According to the class in the comparison result Code information ("U", "I", "X") was used to screen potential lincrna and intronic LncRNA anti-sense LncRNA Coding potential prediction screening: 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.

Differential Expression Analysis
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  (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.

Analysis of Enrichment Pathway of GO and KEGG
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.

Luciferase Reporter Assay
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 cotransfected 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 downregulated reporter gene expression were further screened.

Validation of Sequencing Data by Quantitative Reverse Transcription PCR (qRT-PCR)
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.

RESULTS
To investigate the potential regulatory effects of LncRNA on the growth and development of beef skeletal muscle, we analyzed samples from muscles of different breeds of beef cattle.

Functional Analysis of Differential mRNA
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 upregulated 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.   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.

Prediction of Differentially Expressed LncRNAs Target Genes
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.
FIGURE 7 | Map of binding of bta-miR-133a with LOC104975788 and LOC536229 (Pax7) target site (A). Luciferase reporter assay showed significantly lower luciferase activity of miR-133a co transfected with the psi-LOC112448162 vector compared to the control group (p < 0.01). The luciferase activity of miR-133a cotransfected with psi-LOC536229 (Pax7) vector was also significantly decreased relative to the control Q19 group (p < 0.01) (B, C). These results indicate that both LOC104975788 and LOC536229 (Pax7) contain miR-133a binding sites. Nine differentially expressed LncRNAs and five differentially expressed mRNAs were verified by qRT-PCR (D, E).
Frontiers in Genetics | www.frontiersin.org May 2022 | Volume 13 | Article 849399 LncRNA-miRNA-mRNA Network Interaction Analysis We used Targetscan (Agarwal et al., 2015), Miranda (Enright et al., 2003), and PITA qtar software  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 posttranscriptional 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.

Luciferase Reporter Assay
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.

Verification of Sequencing Results
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.

DISCUSSION
The muscle fiber is the basic unit of muscle, and the type of muscle fiber greatly shapes muscle characteristics (Steinbacher et al., 2006). There were obvious differences in appearance in the muscle fibres isolated from longissimus dorsi muscles of Shandong black cattle and Luxi cattle. The average muscle fiber area of Shandong black cattle was significantly larger than that of Luxi cattle (p < 0.05), and the muscle fiber diameter was smaller than that of Luxi cattle. There was a significant positive correlation between muscle fiber area and carcass traits. Muscle fiber diameter is closely related to meat quality and taste and thicker muscle fiber diameter corresponds to decreased muscle tenderness. There is a negative correlation between muscle fiber diameter and muscle fiber density (Zhou, 2010). The amount of slow muscle fiber affects sarcomere length and has an important impact on meat quality (Hendricks et al., 1971). There were significant differences in the muscle fiber density, the ratio of fast-twitch fibers/slow-twitch fibers, and weight (p < 0.05) in this study. These differences may be key factors leading to the differences in meat production performance and meat quality of the two breeds of cattle after birth, which was also the research basis of this study to explore the underlying molecular regulatory mechanism. With low conservation of LncRNA sequences among species, bioinformatics methods are required to screen and identify LncRNAs. This type of analysis is based on transcript length, the number of exons, and the coding potential. In this study, transcriptome data of skeletal muscle of Shandong black cattle and Luxi cattle were generated and analyzed. In our studies, 1,415 transcripts were found to be differentially expressed, with 939 and 476 transcripts were significantly up-regulated and downregulated in Luxi cattle (Shandong black cattle as a reference group), respectively. Further, 19 GO items and 14 regulatory pathways related to muscle development were screened by hierarchical GO and KEGG cluster analysis. Many genes involved in the regulation of muscle development were identified, including myosin protein family genes (MYH1, MYH3, MYH4, MYH8, and MYL3), myogenic regulatory factor family genes (MY O 10 and MY O 5B), Homeobox family genes (HOXC10, HOXC11, and HOXD9), Troponin T family genes (TNNT2 and TNNT3) and some regulatory transcription factors (WNT4 and TNNT2). It is worth noting that three LncRNAs (LOC107133268, LOC112445823, LOC104969184) are directly enriched into muscle regulation items. Further target gene prediction analysis shows that LncRNAs may participate in the regulation of muscle development through bta-miR-2892 (LOC107133268), bta-miR-2360 and bta-miR-2449 (LOC112445823). The specific regulatory mechanism needs to be further verified and analyzed.
Compared with studies in human (Federica et al., 2019) and other model organisms (James et al., 2015), there has been limited identification and characterization of beef LncRNAs, especially ones related to skeletal muscle development, with most studies of cattly focused on the identification of genes and miRNAs (Xue et al., 2013;Peters et al., 2017). In this study, we identified 8,427 multiple exon LncRNAs in beef skeletal muscle, and 480 differentially expressed LncRNAs were identified. More LncRNAs were detected in this study than previously reported in goats (Sun et al., 2013;Ying-hui et al., 2019). Fifteen randomly selected differentially expressed transcripts were verified by qPCR, and the results were consistent with the results of RNA sequencing. In conclusion, these results confirm the reliability of the identification of LncRNAs (Langfelder et al., 2008).
Although LncRNAs can act on gene sites far from their chromosomal location (Cohen et al., 2000), genes in close proximity on a chromosome may participate in the same cellular metabolic pathways and have similar biological functions (Leonardo et al., 2011). Therefore, the distribution of differentially expressed LncRNAs on chromosomes and the linkage differential expression of nearby genes may have biological significance that can help us to determine the function of a gene. The most differentially expressed LncRNAs were found on chromosome 9, followed by chromosomes 3, 1, 11, and 7. In co-expression analysis, we detected 387 differentially expressed LncRNA transcripts related to protein-coding genes according to the expression correlation coefficient values (r > 0.9), and predicted 1,164 target genes. GO analysis showed 20 GO terms related to the regulation of genes involved in muscle development. KEGG analysis showed enrichment of 1,164 potential target genes in 299 pathways, with some related to muscle development, such as the calcium signaling pathway, the AMPK signaling pathway, the cGMP -PKG signaling pathway, and the PPAR signaling pathway. Interestingly, we also found MYORG, Dll1, EFNB2, SOX6, PROX1, MYOCD, NEBL and MYLK3 genes, annotated as related to muscle development. Overall, LncRNAs may play a regulatory role in skeletal muscle biological development through cis or trans mechanisms.
At the post transcriptional level, LncRNA binds miRNA competitively with mRNA for protein-coding genes, thus relieving the inhibitory effect of miRNA on protein-coding genes to promote expression of these genes (Kartha and Subramanian, 2014). We analyzed the expression patterns of LncRNAs differentially expressed among different breeds, and constructed an LncRNA-miRNA-mRNA interaction network related to skeletal muscle development. A total of 52 LncRNAs related to muscle development were identified, with nine upregulated and five down-regulated. All of these LncRNAs contain one or more putative miRNA binding sites with 48 LncRNAs predicted to interact with one to two miRNAs related to muscle development, and other LncRNAs predicted to interact with more than three miRNAs, such as LOC525506, LOC540051, LOC514189. With the miRNA seed sequences, LncRNAs can bind to miRNA to act as a sponge, preventing miRNA from binding to its target mRNA. As a classic example, M. Cesana (Marcella Davide et al., 2011) confirmed that linc-MD1, a long non-coding RNA specifically expressed during myoblast differentiation, regulates the expression of muscle specific transcription regulators MAML1 and MEF2C through miR-133 and miR-135. In particular, LOC525506 may be involved in the development of beef skeletal muscle by regulating miR-1, miR-23a, miR-378, let-7, miR-483, and miR-21. This is the most predicted LncRNA, but little is known about its expression and function. Both LOC112447073 and LOC104975788 are involved in the skeletal muscle development of beef cattle by interaction with miR-103. More interestingly, some target genes of LOC112447073, LOC104975788, and LOC101903367 directly regulate the development of muscle fiber and maintain the stability and development of muscle, such as miR-103, miR-133a, miR-145, MEF2C, myocardin, and Pax7 (Richard et al., 2008;Mathew et al., 2011). Choi identified 11 LncRNAs in bovine transcripts by studying skeletal muscle and adipose tissue of Korean cattle, with four related to muscle function (Jae et al., 2018). These data provide new insight into the role of LncRNA in muscle development and enrich the existing LncRNA mammalian skeletal muscle resources.
The second to eighth bases of the sequence at the 5 ′end of miRNA is called the seed sequence. Complete complementarity of this sequence and that of the target gene indicates the potential for binding of the target gene by miRNA (Lewis et al., 2005;Andrew et al., 2007). We predicted the miRNA binding sites of functional LncRNAs that may be related to skeletal muscle development. Using a double luciferase test, we found a recognition site of btamir-133a in the sequence of LOC104975788 and a binding site for bta-mir-133a in the downstream target gene Pax7. We speculate that LOC104975788 may be involved in the regulation of skeletal muscle development by competitively inhibiting the expression of the target gene Pax7 through binding of bta-miR-133a. Future work should test this proposed regulatory mechanism.

CONCLUSION
The expression patterns of LncRNAs in skeletal muscle of Shandong black cattle and Luxi cattle were elucidated by RNA sequence analysis, and the LncRNAs that may be involved in the skeletal muscle development of beef cattle were identified. The results allowed construction of interaction networks of LncRNAs-miRNA-mRNA regulated by muscle biology. We speculate that LOC104975788 may be involved in the regulation of skeletal muscle development by competitively inhibiting the expression of the target gene Pax7 through binding of bta-miR-133a.

DATA AVAILABILITY STATEMENT
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. The RNA sequencing data has been uploaded to the NCBI GEO database, and can be retrieved using the following accession numbers: GSM4904154, GSM4904155, 623 GSM4904156, GSM4904157, GSM4904158, GSM4904159.

ETHICS STATEMENT
The animal study was reviewed and approved by Biological Studies Animal Care and Use Committee, Shandong Province, China.

AUTHOR CONTRIBUTIONS
YD, XB, KY, RL and MH designed the study. RL and MH conducted the experiment. RL, MH, XL and XB performed and collected the data. RL analyzed the data and wrote the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was supported by Shandong modern agricultural industry technology system cattle industry innovation team (668-2216030). The team participated in the collection and analysis of data.