Comparison of transcriptome between high- and low-marbling fineness in longissimus thoracis muscle of Korean cattle

Objective This study compared differentially expressed genes (DEGs) between groups with high and low numbers of fine marbling particles (NFMP) in the longissimus thoracis muscle (LT) of Korean cattle to understand the molecular events associated with fine marbling particle formation. Methods The size and distribution of marbling particles in the LT were assessed with a computer image analysis method. Based on the NFMP, 10 LT samples were selected and assigned to either high- (n = 5) or low- (n = 5) NFMP groups. Using RNA sequencing, LT transcriptomic profiles were compared between the high- and low-NFMP groups. DEGs were selected at p<0.05 and |fold change| >2 and subjected to functional annotation. Results In total, 328 DEGs were identified, with 207 up-regulated and 121 down-regulated genes in the high-NFMP group. Pathway analysis of these DEGs revealed five significant (p<0.05) Kyoto encyclopedia of genes and genomes pathways; the significant terms included endocytosis (p = 0.023), protein processing in endoplasmic reticulum (p = 0.019), and adipocytokine signaling pathway (p = 0.024), which are thought to regulate adipocyte hypertrophy and hyperplasia. The expression of sirtuin4 (p<0.001) and insulin receptor substrate 2 (p = 0.043), which are associated with glucose uptake and adipocyte differentiation, was higher in the high-NFMP group than in the low-NFMP group. Conclusion Transcriptome differences between the high- and low-NFMP groups suggest that pathways regulating adipocyte hyperplasia and hypertrophy are involved in the marbling fineness of the LT.


INTRODUCTION
In beef, intramuscular fat (IMF) or marbling is associated with flavor, juiciness, and tenderness and is a critical factor determining the grade of beef in some countries, including Korea and Japan [1,2]. The IMF content is positively correlated with the beef price [3], and highly marbled beef is a premium consumer product, especially in Korea and Japan [4].
Computer image analysis (CIA) systems have been developed to assess various characteristics of marbling particles (MPs), such as their number, size, and distribution (fineness or coarseness) [5]. We recently reported that the number of fine MPs (NFMP) in the longissimus thoracis muscle (LT) was positively correlated with the auction price of Korean cattle beef [6], indicating that Korean consumers prefer beef with fine MPs (each covering an area of 0.01 to 0.5 cm 2 ). The marbling fineness is determined by the marbling fineness index, which is obtained by dividing the NFMP by the ribeye area (cm 2 ) and represents the distribution of fine MPs in the LT. A CIA method has been established for measuring fine MPs, but the mechanism of fine MP formation is unknown.
High-throughput RNA sequencing is a powerful tool for quantifying transcripts [7]. Several RNA sequencing studies have examined the molecular mechanisms of IMF deposition and identified differentially expressed genes (DEGs) between different IMF levels in the LT of cattle [8,9]. The molecular mechanisms responsible for marbling fineness are unknown. Therefore, this study identified DEGs between cattle with high and low NFMP in the LT to understand molecular events associated with formation of fine MPs.

Animals
All experimental procedures followed the Animal Experimental Guidelines of the Seoul National University Institutional Animal Use and Care Committee and were approved by the committee (SNUIACUC: SNU-160408-4). We used LT samples from Korean cattle steers sampled in our previous study [6]. The details of the feeding trial, sampling, and CIA methods can be found in our previous study [6]. Briefly, the size and distribution of MPs were assessed using a CIA method in 43 Korean cattle steers. Of all LT samples, 10 were selected based on the NFMP and assigned to either a high-(n = 5) or low-(n = 5) NFMP group ( Table 1). The NFMP was used as the selection criterion, because it is highly correlated with auction price of beef and is positively correlated with quality grade (R 2 = 0.41, p<0.001) [6]. Three samples from each group were randomly selected for RNA sequencing analysis, and all 10 samples were used for quantitative real-time polymerase chain reaction (qPCR) analysis. (Table 1). Each three steers from the high and low group were randomly selected for RNAseq, and total ten steers were used for qPCR to validate the result of RNA-seq.

Carcass grading and the computer image analysis method
The carcass traits (back fat thickness, ribeye area, marbling score, and quality grade) were graded by trained meat graders according to the Korean carcass grading system of the Korea Institute for Animal Products Quality Evaluation [10]. The carcass auction price was determined by wholesalers after grading in the slaughter house. Detailed methods have been described previously [6].
Cross-sectional photographs of the ribeye area in the LT were taken alongside the 13th thoracic vertebra using HK-333 beef carcass photography equipment (mirror type) (Hayasaka Ricoh, Sapporo, Japan) developed by Kuchida et al [11], and the photographs were analyzed using Beef Analyzer II (Hayasaka Ricoh, Japan), as described previously [6,12]. CIA traits were obtained as described previously [6]. Briefly, the total number of MPs was counted from the binarized fat image. Fine MPs were defined as those covering an area between 0.01 and 0.5 cm 2 . The fineness index was obtained by dividing the NFMP by the ribeye area (cm 2 ).

RNA sequencing
Total RNA was extracted from LT tissue using TRIzol reagent (Sigma-Aldrich, Milwaukee, WI, USA) in accordance with the manufacturer's protocol. Total RNA was quantified with a ND-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) by recording the absorbance at 260 nm. The RNA quality was checked via agarose gel electrophoresis and ethidium bromide staining of the 28S and 18S bands, using an Agilent 2100 Bioanalyzer with the RNA  1) These data were reported in our previous study [6]. CIA, computer image analysis; qPCR, quantitative real-time polymerase chain reaction; NFMP, number of fine marbling particles; SEM, standard error of mean; MP, marbling particle. 2 6000 Nano Chip (Agilent Technologies, Santa Clara, CA, USA). RNA samples were stored at -70°C until RNA sequencing and qPCR analyses. The cDNA library was constructed and sequenced at ebiogen (Seoul, Korea). The library was constructed using a SENSE 3′ mRNA-Sequencing Library Prep kit (Lexogen, Vienna, Austria) [13] according to the manufacturer's instructions. Briefly, an oligo-dT primer containing an Illumina-compatible sequence at its 5′ end was hybridized with 500 ng of total RNA, and reverse transcription was initiated. After degrading the RNA template, the second strand was synthesized using a random primer containing an Illumina-compatible linker sequence at its 5′ end.
To remove all reaction components, the double-stranded library was purified using magnetic beads. The library was then amplified to add the complete adapter sequences required for cluster generation. High-throughput sequencing was conducted as single-end 75-bp sequencing using a NextSeq 500 system (Illumina, San Diego, CA, USA). The quality of the sequencing reads from LT tissues was checked using FastQC. The reads that passed quality control were mapped to the Bos taurus genome (UMD3.1) from the UCSC database using Tophat2 (v2.0.2) and counted using HTSeq (v0.5.3p3).

Identifying and annotating differentially expressed genes
For data analysis, the SENSE 3′ mRNA sequencing reads were aligned using Bowtie2 ver. 2.1.0. [14]. Bowtie2 indices were generated from either the genome assembly sequence or representative transcript sequences for alignment to the genome and transcriptome. The alignment file was used to assemble transcripts, estimate transcript abundances, and detect DEGs. DEGs between the high-and low-NFMP groups were determined based on counts from the unique and multiple alignments. The read count data were processed based on the quantile normalization method using Genowiz ver. 4.0.5.6 (Ocimum Biosolutions, Hyderabad, India  [17]. qPCR data were normalized relative to the level of the β-actin housekeeping gene, as validated in our previous study [16].

Statistical analysis
Statistical differences between the two groups were determined using the independent t-test or Mann-Whitney test depending on data normality with IBM SPSS ver. 22 (SPSS, Chicago, IL, USA).

Carcass and computer image analysis characteristics of the animals
We compared the carcass and CIA characteristics between the high-and low-NFMP samples used RNA sequencing and qPCR ( Table 1). The high-NFMP group had a higher marbling score, NFMP, and fineness index than did the low-NFMP group (all p<0.001). The meat from the high-NFMP group also fetched a higher auction price (p<0.001), confirming the result in our previous report [6].

RNA sequencing and identification of differentially expressed genes
Six cDNA libraries were constructed, and an average of 20,585,914 single-end 75-bp sequence reads was generated. At least 95.1% of the reads in each sample were mapped to the bovine reference genome (UMD3.1) using Tophat2 (v2.0.2) (Supplementary Table S2). In total, 328 DEGs were identified (p<0.05 and |FC| >2), with 207 up-regulated and 121 downregulated genes in the high-NFMP group (Figure 1). The upand down-regulated genes included 97 and 122 proteincoding genes, respectively. Table 2 lists the top 10 up-and down-regulated genes.

Functional annotation of differentially expressed genes
The DAVID tool was used for functional annotation of the 328 DEGs, which were categorized as related to biological process, molecular function, and cellular component (Supplementary Tables S3 to S5). The number of genes associated with each pathway was small (Table 3). Nevertheless, there were significant pathway terms for DEGs, including endocytosis, protein processing in endoplasmic reticulum (ER), and adipocytokine signaling pathway, which are thought to affect adipocyte hypertrophy and hyperplasia.
Eight genes were found to be involved in the endocytosis pathway. Of these, three were up-regulated in the high-NFMP group: beta-arrestin-1 (ARRB1; FC = 2.07, p = 0.019), vacuolar protein sorting-associated protein 37D (VPS37D; FC = 2.54, p = 0.015), and ADP-ribosylation factor 3 (ARF3; FC = 2.35, p = 0.031). The other five genes involved in the endocytosis pathway in the high-NFMP group were downregulated: ADP-ribosylation factor 4 (ARF4; FC = 0.47, p = 0.030), ARF-GAP with coiled-coil, ANK repeat and PH domain-containing protein 2 (ACAP2; FC = 0.31, p = 0.009), heat shock 70 kDa protein 1B (HSPA1A; FC = 0.36, p = 0.005), heat shock 70 kDa protein 6 (HSPA6; FC = 0.24, p = 0.008), and vacuolar protein sorting-associated protein 45 (VPS45; FC = 0.49, p = 0.013). Endocytosis is a vital cellular process in which extracellular ligands, soluble molecules, proteins, and lipids are taken up from the extracellular matrix [18]. Glucose transporter type 4 (GLUT4) is the most abundant glucose transporter regulated by insulin in adipocytes [19]. GLUT4 is translocated from the intracellular matrix to the plasma membrane or vice versa depending on insulin levels; this is orchestrated by exocytosis and endocytosis [20]. Both GLUT4 and CD36, an essential protein for fatty acid endocytosis and uptake in adipocytes, have similar trafficking and recycling systems induced by insulin [21]. Endocytosis is a crucial cellular process through which the uptake of glucose and fatty acid occurs. Through lipogenesis, adipocytes convert glucose into long-chain fatty acids and triglycerides, and finally store lipids as lipid droplets (LDs) [22]. Endocy- plot of RNA-sequencing data between the high and low NFMP groups. Up-regulated genes are colored red, and enes are colored green in the high NFMP group. Two fold-up and down-regulated threshold lines were shown as and green line, respectively, and P-value (< 0.05) threshold line was shown as dotted horizontal line. tosis related to GLUT4 and CD36 trafficking may lead to enlargement of the adipocytes (hypertrophy), which may in turn induce adjoining preadipocytes to proliferate (hyperplasia) by secreting paracrine growth factors [23]. Therefore, endocytosis may directly or indirectly affect marbling fineness. All five DEGs involved in protein processing in the ER pathway were down-regulated in the high-NFMP group: DNAJ homolog subfamily A member 1 (DNAJA1; FC = 0.48, p = 0.003), DNAJ homolog subfamily B member 1 (DNAJB1; FC = 0.48, p = 0.011), endoplasmic reticulum oxidoreductase 1 beta (ERO1B; FC = 0.46, p = 0.042), HSPA1A, and HSPA6. The ER is a critical cellular organelle controlling calcium homeostasis, lipid synthesis, protein folding, and LD formation, and the size and number of LDs can be affected by the ER stress response [24]. The protein-folding or export capacity of the ER can be exceeded when the ER is stressed under some physiological and pathological conditions [25]. Several studies have reported LD accumulation [26], reduced adipocyte insulin sensitivity [27], and adiponectin oligomerization [28] during the ER stress response. Therefore, the involvement of ER protein processing in adipocyte hyperplasia and hypertrophy may contribute to marbling fineness.
All four DEGs involved in the adipocytokine signaling pathway were up-regulated in the high-NFMP group: long chain fatty acid CoA ligase ACSBG2 isoform X1 (ACSBG2; FC = 4.91, p = 0.002), 5′-AMP-activated protein kinase subunit gamma-2 isoform X1 (PRKAG2; FC = 2.15, p = 0.048), leptin receptor (LEPR; FC = 3.70, p = 0.007), and solute carrier family 2, facilitated glucose transporter member 1 (SLC2A1 or GLUT1; FC = 2.72, p = 0.019). Leptin is mainly secreted by mature adipocytes [29], and high adiposity and body fat increased the concentration and expression of leptin in cattle [30]. Leptin also can induce LD formation in several cell types in an mammalian target of rapamycindependent manner [31,32]. Genetically modified mice without the leptin receptor in adipose tissue were lighter than controls [33], and bone marrow stromal cells with impaired leptin receptors did not differentiate into adipocytes [34]. These results suggest that leptin can directly affect the development and function of adipocytes via leptin receptors. SLC2A1 is a class I glucose transporter that is expressed in adipose tissue at moderate levels [35]. PRKAG2 is a subunit of 5′-AMPactivated protein kinase [36]. Collectively, the adipocytokine signaling pathway may affect adipocyte hyperplasia or hypertrophy, possibly contributing to marbling fineness.

Validation of differentially expressed genes using quantitative real-time polymerase chain reaction
To identify candidate genes associated with marbling fineness, we selected 10 DEGs which may be related with adipocyte hyperplasia and hypertrophy, although these were not listed in the top 10 up-and down-regulated genes shown in Table  2. Description of gene function in UniProt website was considered as selection standard (The UniProt Consortium [37]). The qPCR results were consistent with those from the RNA sequencing analysis for all genes tested, although statistical significance was not obtained in the qPCR analysis (Table 4).
We observed higher insulin receptor substrate 2 (IRS2) transcript levels in the LT samples from the high-NFMP group. IRS2 plays an important role in the metabolic activities of insulin and IGF1 in the liver, muscle, and adipose tissue [42]. Cell culture studies suggest that IRS2 also plays an important role in adipocyte differentiation [43]. Therefore, the up-regulated IRS2 expression in the high-NFMP group suggests that it is involved in the formation of MPs in the LT.
Several studies identified DEGs by RNA sequencing analysis and reported transcriptome or genes associated with the marbling score. For example, the DEGs from the contrast between age and marbling score were recently identified by RNA sequencing analysis in longissimus dorsi muscle of Korean cattle steers, and new potential early age markers and many genes for marbling development were found [44].
In this study, we compared DEGs between groups with highand low-NFMP in the LT of Korean cattle, and the NFMP is used for calculation of fineness index. Our study has a novel aspect that we identified DEGs associated with marbling fineness, which is a new concept for marbling score.

CONCLUSION
This study is a first attempt to understand the molecular events associated with the formation of fine MPs in the LT. Using RNA sequencing and KEGG pathway analyses, we found that three pathways (adipocytokine signaling pathway, endocytosis, and protein processing in the ER) involved with hyperplasia, hypertrophy, and adipocyte differentiation were significantly enriched in the DEGs. The SIRT4 and IRS2 genes, which are involved in adipocyte hyperplasia and hypertrophy, were up-regulated in the high-NFMP group. Our study suggests that pathways regulating adipocyte hyperplasia and hypertrophy are involved in marbling fineness. Further study is required to determine whether the balance between adipocyte hyperplasia and hypertrophy is linked to marbling fineness.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.

FUNDING
This work was supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2017R1A2B4003207; 2020R1I1A1A01052628).  The qPCR results were consistent with those from the RNA sequencing analysis for all genes tested, although statistical significance was not obtained in the qPCR analysis.