Comparison of Long Non-Coding RNA Expression Profiles of Cattle and Buffalo Differing in Muscle Characteristics

Buffalo meat consist good qualitative characteristics as it contains “thined tender” which is favorable for cardavascular system. However, the regulatory mechanisms of long non-coding RNA (lncRNA), differences in meat quality are not well known. The chemical-physical parameters revealed the muscle quality of buffalo that can be equivalent of cattle, but there are significant differences in shearing force and muscle fiber structure. Then, we examined lncRNA expression profiles of buffalo and cattle skeletal muscle that provide first insights into their potential roles in buffalo myogenesis. Here, we profiled the expression of lncRNA in cattle and buffalo skeletal muscle tissues, and 16,236 lncRNA candidates were detected with 865 up-regulated lncRNAs and 1,296 down-regulated lncRNAs when comparing buffalo to cattle muscle tissue. We constructed coexpression and ceRNA networks, and found lncRNA MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 could be as competitive endogenous RNA (ceRNA) containing potential binding sites for miR-1/206 and miR-133a. Tissue expression analysis showed that MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 were highly and specifically expressed in muscle tissue. Present study may be used as a reference tool for starting point investigations into the roles played by several of those lncRNAs during buffalo myogenesis.


INTRODUCTION
As a form of striated muscle tissue, skeletal muscle is an important object in the study of meat quality and plays key roles in regulating homeostasis and metabolism, accounting for 40-60% of the animal body (Li et al., 2018b). Skeletal muscle developmental inability via perantes lead to embryonic death, while the failure to repair or maintain skeletal muscle after birth that can lead to a decline in quality of life and even death. Myogenic progenitor cells from multi-potent mesodermal precursor cells are committed to the muscle fate, and express Pax3 and Pax7 Abbreviations: CeRNA, competing endogenous RNA, miRNA, microRNA; lncRNA, long non-coding RNA; qRT-PCR, realtime quantitative PCR. destined to become myoblasts (Buckingham and Relaix, 2015). Then, myoblasts experience proliferation, differentiation and fuse into myotubes through the regulation of myogenic regulatory factors: MyoD, MRF4, Myf5, and myogenin play key roles in the process of regeneration in adult muscle (Hernandez-Hernandez et al., 2017). In fact, various noncoding RNAs (ncRNAs) have been identifi ed and demonstrated to regulate myogenesis and muscle regeneration including microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) (Li et al., 2018b). MiRNAs profoundly influenced the physiology and pathology of skeletal muscle, such as miR-1 and miR-133 in vitro and in vivo have distinct roles in regulating skeletal muscle proliferation and differentiation (Li et al., 2018a).
LncRNAs have played versatile roles in regulating skeletal myogenesis and regeneration at multiple levels (Kallen et al., 2013;McHugh et al., 2015;Munschauer et al., 2018;Neumann et al., 2018). LncRNAs are well known for their involvement in transcriptional/ epigenetic regulation on chromatins through interacting with chromatin regulators, for example, acting as "molecular scaffold" or decoys to activate or repress transcription (Caretti et al., 2006;Korostowski et al., 2012). Other unique mechanisms have been found to explain the diverse modes of lncRNA action in myogenesis (Han et al., 2014;Zhou et al., 2015). For example, lncRNA could act as competing endogenous RNAs (ceRNAs) to sequester miRNAs from their target mRNAs (Cesana et al., 2011;Kallen et al., 2013). Some transcribed from antisense strand of protein coding genes could directly pair with the mRNA to modulate coding gene translation (Wang et al., 2013;Wang et al., 2016). Interestingly, Emerging research has shown that some lncRNAs could translate micropeptides (<100 amino acids) to perform micropeptidemediated functions (Anderson et al., 2015;Nelson et al., 2016). However, participation of miRNA and lncRNA in muscle development are still in their infancy, especially in livestock, for example, in buffalo studies.
Buffaloes usually used for labor purposes, but now optimized for meat or milk production (Pisano et al., 2016;Low et al., 2019). In developing countries, buffalo meat usually from old (>10 years old) period is eaten and therefore seemed tougher than beef (Sakaridis et al., 2013;Zhang et al., 2016). Compared to beef, buffalo meat indeed has less fat, lower calories, and less cholesterol, which is healthier and can confer significant cardiovascular benefits. "Buffalo meat is an amazing cure for diabetes" as is described in the Compendium of Materia Medica, which was considered to be the most comprehensive and complete medical work in the history of Chinese medicine. Here, the longissimus dorsi muscles of swamp buffaloes and Guangxi native cattle under the same feeding and management were selected, and analyzed its differences in physiological biochemical indexes. The histological staining and analytical chemistry methods were used to directly compare the differences in pH, water content, shear force, intramuscular fat content, ash content, and myofiber structure.
In the second part, the Ribo-Free RNA-seq method was selected to analyze the expression of lncRNA in longissimus dorsi muscles of swamp buffaloes and Guangxi native cattle. Differentially expressed genes and lncRNAs were identified in skeletal muscle samples, and the candidate lncRNAs were verified by Quantitative PCR(qPCR). We further constructed coexpression and ceRNA networks to select candidate lncRNA. Our research will be beneficial for the improvement of Chinese meat buffalo breeding and provide new insights into the genetic mechanism of Chinese swamp meat.

Sample Preparation
Chinese swamp buffalo (n = 3, 12 months old) and Guangxi native cattle (n = 3, 12 months old) under the same feeding and management were obtained from SIYE husbandry of Guangxi, China. The longissimus dorsalis muscle of adult buffalo and cattle were used for muscle quality analysis, transcriptome sequencing and qPCR analysis. The 4 months old buffalo and cattle fetal with a body length of 15 cm were selected from the local slaughterhouse in Nanning. Buffalo and cattle fetal tissues (skin, heart, liver, spleen, lung, kidney, small intestine, leg muscle, and longissimus dorsalis muscle) were used to extract RNA and analyze the expression of lncRNAs.

Meat Quality Evaluation
The longissimus thoracis muscles were taken between the 12th and 14th ribs from the left side of body, and performed the following analyses in triplicate: the pH was immediately measured using a pH meter (Thermo Orion, Hudson, NH, USA); water content was determined on drying at 100°C for 24 h; crude fat level was assessed by extracting for 12 h using petroleum ether; ash content was evaluated by ashing at 600°C for 10 h; shear force was measured using a C-LM3 digital display tenderness instrument (Northeast Agricultural University, Harbin, China). The amino acid composition in muscle samples was determined using an amino acid analyzer. Differences between the two groups were compared using a post hoc test.

Library Preparation
Total RNA of longissimus muscle was extracted and assessed by electrophoresis and quantified by NanoDrop spectrophotometer (NanoDrop, Wlinington, USA). Ribosomal RNA was removed by probe, and then the remaining RNA was used for library construction and sequencing (Ribo-Zero RNA-seq). cDNA library preparation and Illumina sequencing analysis were perofrmed as previous studies (Hui et al., 2018).

lncRNAs Identification
Potential lncRNAs were filtered through the following highly stringent criterion: (1) transcript length is not less than 200 nt; (2) transcript expression is more than 3 reads; (3) the transcripts were annotated as "i", "j", "o", "u", and "x" according to the cuffcompare classes; (4) the coding potential calculator (CPC) score less than -1, and coding-non-coding-index (CNCI) score less than 0 were kept; (5) the transcripts containing open reading frame (ORF) is greater than 100 aa were removed; (6) via aligning to the Swiss-Protein, Cpat, and Pfam database, the transcripts were removed with known protein-coding domain.

Coexpression Analysis
As a cis-regulator, lncRNA could regulate the expression of adjacent genes. The coexpression network of the candidate lncRNAs and their upstream or downstream 10 kb mRNAs was constructed. The connectivity and enrichment were performed due to Position Frequency Matrix.

Quantitative Real-Time PCR (qRT-PCR)
Total RNA were extracted using Trizol reagent (TaKaRa, Dalian, China), and reverse transcription was performed by HiScript R II One Step RT-PCR kit (Vazyme, Nanjing, China). qRT-PCR was performed with ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China) with the internal control of b-actin using 2 -DDCt method. All primers sequences were showed in Table S1.

Comparison of Meat Quality
In order to understand the difference of meat quality between buffalo and cattle, we analyzed the physiological and biochemical indexes of their longissimus dorsi samples. According to the indexes, significant differences were found in shear force between buffalo and cattle (P < 0.05), while no significant differences were found in pH, water content, intramuscular fat content and ash content (P > 0.05; Table 1). Moreover, the muscle samples were made into frozen sections and observed with HE staining by microscopy, and showed that the muscle fiber area, isometric diameter, circumference and density of buffalo were significantly smaller than those of cattle (P < 0.05; Table 2 and Figure 1).
Furthermore, the amino acid composition of muscle samples was also analyzed, and showed no difference in the ratio of essential amino acids or umami amino acids (P > 0.05), suggesting that under the same feeding and management conditions, the longest dorsal muscle of buffalo and cattle are basically similar in amino acid composition ( Table 3). These results revealed that under the same feeding and management, the muscle quality of buffalo can be equivalent to that of cattle, but there are significant differences in shearing force and muscle fiber structure. Therefore, whole transcriptome RNA-Seq was performed to analyze the differences of buffalo and cattle musculus longissimus.

Ribo-Zero RNA-Seq of Buffalo and Cattle Muscle
Three longissimus muscle samples were selected to perform Ribo-Zero RNA-Seq from cattle and buffalo at 12 M old, respectively. As shown in Figure 2A, a large number of lncRNAs were identified. On average, 83~137 and 116~131 million unique mapped clean reads were acquired from the buffalo and cattle libraries, respectively ( Table 4). We found 67.5% of the reads located in exon regions, while a significant reduction was observed in intergenic or intronic regions (32.5%; Figure 2B). Novel reliable lncRNAs were filtered by using Pfam and Cpat database and tested by CPC and CNCI, and a total of 16,236 potential lncRNA transcripts were identified to be expressed ( Figure 2C, Table S2). We found chromosome with longer length to be more likely to produce more lncRNAs, indicating that the number of reads distributed in the chromosome increased with chromosome length ( Figure 2D, E). According to the Cuffcompare classes, the lncRNAs aligned to intergenic regions (u) accounted for the largest proportion (8,605, 53%; Figure 2F).

Genomic Features of Identified LncRNAs
As shown in Figure 3A, B, the identified lncRNAs showed a low expression level, and the mean expression levels were 5.96 (FPKM). As illustrated in Figure 3C, the lncRNA data showed a good correlation between buffalo and cattle muscle samples. 2,161 lncRNAs were significantly (P < 0.05) differently expressed in buffalo and cattle muscle samples ( Figure 3D), and were listed in Table S3. Among all the differentially expressed lncRNAs, MSTRG.233222.1 and MSTRG.104517.1 showed the highest expression level of all up-regulated and down-regulated lncRNAs when comparing buffalo to cattle muscle tissue, respectively. To better understand potential functions of lncRNA, the scatter plot and volcano plot were displayed in Figures 3G-I. There were 865 lncRNAs were up-regulated, while  1,296 lncRNAs were down-regulated (P < 0.05; Table S3), and buffalo had a clear tendency for low expression of lncRNA ( Figures 3D-I).
LncRNA regulates the transcription of coding genes through cis and trans regulatory relationships: if the role of lncRNA is limited to the same chromosome (adjacent genes), cis regulation is exercised; trans works when it affects gene expression on other chromosomes (at long distances). The top 30 enriched KEGG pathways by cis and trans regulation were present in Figures 3E, F. Calcium signaling pathway and Valine, leucine and isoleucine degradation had the highest level in cis or trans' target genes cluster, respectively, indicating that these pathways may involve in regulating skeletal myogenesis.
Previous study has shown that lncRNA is shorter in length than protein-encoded transcripts. As illustrated in Figure 4A, the mean length of lncRNAs (1,087 nucleotides) was shorter than that of the mRNA (1,153 nucleotides). The average ORF length of lncRNA was 66.7 nt, and the mRNA was 309.7 nt, revealing that lncRNA has a lower coding potential than protein-coding genes ( Figure 4B). Moreover, lncRNAs had fewer exons (about 2.4) than protein-coding genes (about 3.6, Figure 4C). Interestingly, the average expression level of lncRNA was approximately 2.5-fold higher than that of protein-coding genes (6.0 vs 2.4; Figure 4D), revealing that lncRNAs in muscle do not act as transcriptional noise and may play important roles in regulating biological processes. Similar to mRNA, lncRNA has a similar number of isoforms, suggesting its important roles in transcriptional regulation ( Figure 4E). The expression of lncRNA and mRNA in different combinations of comparative analysis indirectly shows the expression relationship between lncRNA and mRNA in a certain biological period, thus volcano and MA interactive maps of differentially expressed lncRNA and mRNA were drawn ( Figures 4F, G). The positional relationship of genes on chromosomes is closely related to the functions of genes, and some lncRNAs may have regulatory functions on their neighboring genes. Therefore, we analyzed the chromosome distribution of differentially expressed lncRNAs and mRNAs ( Figure 4H).

Coexpression and CeRNAs Networks
In order to further investigate the cis-regulatory relationship of lncRNAs, the adjacent coding genes 10 kb upstream and downstream of the candidate lncRNAs were performed to construct coexpression network. The top 10 most significantly differentially expressed lncRNAs were chosen to hunt their neighboring coding genes ( Figure 5A). Each lncRNA had different number of adjacent coding genes. For example, MSTRG.266281.11 had maximum number of 11 neighboring coding genes, whereas MSTRG.203788.46 had only one nearby coding gene (MYH8) and was positively correlated with FIGURE 1 | Histological analysis of longissimus dorsi muscles. The HE staining of muscle results showed that the muscle fiber area, isometric diameter, circumference and density of buffalo (A) were significantly smaller than those of cattle (B). expression levels of MYH8. Interestingly, MSTRG.233222.1 were up-regulated in the buffalo muscle compared to the cattle muscle, and all its four neighboring coding genes (ICK, GSTA5, FBXO9, GSTA3) presented higher levels in the buffalo muscle, suggesting this lncRNA may has cis-regulatory relationship on its neighboring genes. The coexpression network could furnish valuable clue for these lncRNAs' potential function in regulating nearby coding genes.
LncRNAs to sequester miRNAs from their target mRNAs could be as a member of ceRNAs, and miRNAs act as common target of the lncRNAs and mRNAs. 15 muscle development related miRNAs were chosen with a total of 5 lncRNAs and 44 mRNAs to construct an ceRNA (mRNA-miRNA-lncRNA) network ( Figure 5B). For instance, MSTRG.30030.4 has multiple binding sites for muscle-related microRNAs, for example, miR-133a and miR-128 (Chen et al., 2006;Kim et al.,

2006
). This ceRNA network may provide valuable information for buffalo skeletal myogenesis.

Identification of Candidate LncRNA
To further validate lncRNAs expression profiles obtained from the RNA-Seq results, 14 lncRNAs that may be involved in muscle development regulation were selected and measured by qRT-PCR. The normalized read counts of the 14 lncRNAs were shown in Table S2. Overall, these randomly selected lncRNAs showed similar expression patterns between qRT-PCR and sequencing results, suggesting that the lncRNA-Seq data are highly accurate ( Figure 6A). Similarly, we also analyzed the expression of these 14 lncRNAs in the fetal dorsal longest muscle and found that the expression of these 14 lncRNAs varied more in the fetal period, indicating that muscle development was more complex in the fetal period and lncRNA was involved in this process ( Figure  6B). We also analyzed the changes of lncRNA in fetal leg muscles and found that the change trend was basically consistent with the expression of lncRNA in the longest dorsal muscle ( Figure 6C) As shown in Figure 6D, MSTRG.30030.4 and MSTRG.104517.1 were the lncRNAs with the most significant expression differences among the highly expressed lncRNAs, revealing that they may play important roles in muscle growth or differentiation. Therefore, we also examined the expression of 14 lncRNAs in different tissues in order to find potential lncRNAs that are specifically expressed and highly expressed in muscles. We examined the expression of these 14 lncRNAs in the heart, spleen, lung, liver, kidney, skin, small intestine, brain, leg muscle and dorsal longest muscle of cattle and buffalo. The buffalo tissue     .30030.4,MSTRG.104517.1, MSTRG.48330.7, MSTRG.58 818.1, MSTRG.71408.1, MSTRG.203788.46, and MSTRG.233222.1 was highly expressed in muscle tissue and low in other tissues, indicating potential roles in buffalo muscle development (Figure 7). Similarly, we found that MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 were highly and specifically expressed in cattle muscle tissue, and these lncRNAs could be chosen as candidates to analyze their real roles in vivo and in vitro in muscle development (Figure 8).

DISCUSSION
Muscle strength is a quantitative trait, which is related to a variety of physiological and biochemical indexes. The origin and evolution of buffalo are closely related to the cultivation of human beings (Low et al., 2019). Due to the rise of industrial revolution and the improvement of social productivity, few buffaloes are still retained for farming. At present, buffaloes were selected for meat or milk production, and their easement value is gradually abandoned (Pisano et al., 2016). By analyzing the longest dorsal muscle of buffalo and cattle in the same breeding and growing environment, it was found that their meat value was comparable. Therefore, high quality buffalo meat can be obtained through optimization of breeding management. Muscle freshness is related to the composition and content of umami amino acids in muscle. The flavor of beef is related to fatty acid content, and marbling level of beef is one of the important indicators in beef classification. Present study, focus on physiological and biochemical indexes, amino acid composition, and intramuscular fat content of the longest dorsal muscle of boar buffalo and local cattle were not significantly different under the same feeding and management conditions. Interestingly, there are significant differences in the  shear force and muscle fiber structure (muscle fiber area, diameter, and circumference) between buffalo and cattle, which may be due to genetic factors rather than the influence of breeding management on this trait. These results suggest that the strength trait of buffalo is positively selected, which is related to the role of buffalo in providing animal power under China's small-scale peasant economy for thousands of years. Similarly, candidate genes associated with strength trait were also positively selected.
Most of the studies on the molecular mechanism of skeletal muscle development in bovine are protein coding genes. However, the occurrence and potential functions of lncRNAs, which reflect the differences between the longissimus muscles of buffalo and cattle, are still unclear. Abundant lncRNAs were differentially expressed in the muscle tissue of buffalo and cattle, suggesting that lncRNAs have specific roles in muscle but not by-product of mRNA. In addition, some lncRNAs were specifically or mainly expressed in muscle tissue, such as MSTRG.30030.4, revealing that these lncRNAs are purposefully produced. LncRNA is more than just a by-product of protein coding genes, and many lncRNAs have been demonstrated to play a role in skeletal muscle development. Increasing studies show that lncRNA in cis or in trans is involved in the transcriptional or post-transcriptional regulation of gene expression (Caretti et al., 2006;Korostowski et al., 2012;Han et al., 2014;Zhou et al., 2015).
The top 10 most significantly differentially expressed lncRNAs were chosen with their neighboring coding genes to construct a co-express network. For example, MSTRG.233222.1 were up-regulated in the buffalo muscle, and all its four neighboring coding genes presented higher levels in the buffalo muscle, suggesting this lncRNA may has cis-regulatory relationship on its neighboring genes. AmRNA-miRNA-lncRNA network was constructed in buffalo muscle according to the common target miRNAs of the mRNAs and lncRNAs. MSTRG.30030.4 has multiple binding sites for muscle-related microRNAs, for example, miR-133a and miR-128. Tissue expression analysis showed that lncRNA MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 were mainly expressed in muscle tissue, that revealing its potential role in buffalo muscle development. Moreover, MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 had multiple binding sites for muscle-related microRNAs, for example, miR-1/206 and miR-133a which were the most representative muscle-associated miRNAs (Chen et al., 2006;Kim et al., 2006). Therefore, our next step is to explore the role of MSTRG.48330.7, MSTRG.30030.4, and MSTRG.203788.46 in the differentiation of cattle and buffalo myoblasts.

CONCLUSIONS
This study is the first to compare chemical-physical characteristics of muscle in cattle and buffalo, and provide an overview of lncRNA expression in buffalo muscle tissues. Thousands of lncRNAs were identified, and several of which were highly and specifically expressed in buffalo muscle tissues. We further constructed coexpression and ceRNA networks, and f o u n d M S T R G . 4 8 3 3 0 . 7 , M S T R G . 3 0 0 3 0 . 4 , a n d MSTRG.203788.46 could be as ceRNA which contained potential binding sites for miR-1/206 and miR-133a. This study may lay a foundation for in-depth investigations into the roles of those lncRNAs during buffalo muscle development.

DATA AVAILABILITY STATEMENT
The data generated and transcripts obtained were deposited at NCBI as the SRA accession SRP116252 and TSA accession GFWP00000000.1. The data and material are also provided as Supplementary Data.