Gene Expression Analysis Provides New Insights into the Mechanism of Intramuscular Fat Formation in Japanese Black Cattle

Japanese Black cattle (Japanese Wagyu) have a unique phenotype in which ectopic intramuscular fat accumulates in skeletal muscle, producing finely marbled beef. However, the mechanism of intramuscular fat formation in Japanese Black cattle remains unclear. To investigate the key genes involved in intramuscular fat accumulation, we comprehensively analyzed mRNA levels in subcutaneous and intramuscular fat tissues using RNA sequence (RNA-seq) analysis, which detected 27,606 genes. We identified eight key genes, namely carboxypeptidase E, tenascin C, transgelin, collagen type IV alpha 5 (COL4A5), cysteine and glycine-rich protein 2, PDZ, and LIM domain 3, phosphatase 1 regulatory inhibitor subunit 14A, and regulator of calcineurin 2. These genes were highly and specifically expressed in intramuscular fat tissue. Immunohistochemical analysis revealed a collagen network, including COL4A5, in the basement membrane around the intramuscular fat tissue. Moreover, pathway analysis revealed that, in intramuscular fat tissue, differentially expressed genes are related to cell adhesion, proliferation, and cancer pathways. Furthermore, pathway analysis showed that the transforming growth factor-β (TGF-β) and small GTPases regulators RASGRP3, ARHGEF26, ARHGAP10, ARHGAP24, and DLC were upregulated in intramuscular fat. Our study suggests that these genes are involved in intramuscular fat formation in Japanese Black cattle.


Introduction
Japanese Black cattle (i.e., Japanese Wagyu) are among the most expensive meats, characterized by excellent marbling, rich and sweet aroma (the so-called Wagyu beef aroma) [1,2]. Japanese Black cattle are pure Wagyu species that originate from pedigrees, such as Tajima, Kedaka, and Itozakura in Japan [3], whereas hybrid breeds with other cattle species are widespread all over the world. Marbling is a phenotype in which ectopic intramuscular fat accumulates in muscle tissue. In humans, excessive dietary lipid intake and age-related senescence lead to ectopic intramuscular fat in skeletal muscle. Fatty acid accumulation correlates with lipotoxicity-related cell dysfunction and is associated with sarcopenic obesity [4]. However, unlike cattle, humans accumulate intramuscular fat due to senescence and illness. Therefore, the marbling traits of Japanese Black cattle are considered a valuable research model for elucidating the molecular mechanism of ectopic intramuscular fat formation in various livestock, as well as in humans [5].
Intramuscular fat is thought to be formed by the myogenic transdifferentiation of myogenic stem cells [6] or a multi-step differentiation process from fibroblast-like preadipocytes to mature adipocytes [7]. The preadipocytes also regulate proliferation, differentiation, and lipogenesis with surrounding connective tissues through intercellular communication via paracrine factors, extracellular matrix (ECM), and cell-cell adhesion [8]. A genome-wide analysis has been conducted to study the genetic traits of intramuscular fat in livestock [9,10]; however, the molecular mechanism of intramuscular fat formation remains unknown. Intramuscular fat formation is a quantitative trait affected by environmental factors and genomic sequences. The difficulty in identifying key genes related to intramuscular fat accumulation is the requirement to consider environmental factors. Therefore, we performed a transcriptome analysis to attempt to elucidate this molecular mechanism. Recent improvements in the comprehensive analysis of large-scale RNA sequencing (RNAseq) data have been generated, and several research institutes attempted to study mRNA expression profiles in cattle [11][12][13][14]. However, RNA-seq analysis has not been performed to report the transcriptome associated with pure Japanese Black cattle.
To search for the key genes involved in intramuscular fat accumulation, we performed comparative RNA-seq analyses of subcutaneous and intramuscular fat tissues from Japanese Black cattle, which exhibit low genetic and environmental variation. Additionally, we performed quantitative PCR (qPCR) analysis and tissue immunostaining of the identified characteristic genes to explore the molecular mechanism associated with intramuscular fat formation. Furthermore, we examined the molecular switches related to the signal transduction pathways involved in intramuscular fat formation using pathway analysis.

Japanese Black Cattle
In this study, we analyzed the intramuscular and subcutaneous fat of four Japanese Black cattle. The four cattle are of a similar breed. Their pedigree is a hybrid of the Tajima, Kedaka, and Itozakura strains, exhibiting excellent bodyweight growth [1,15,16]. The sampled cattle were selected from representative breeds (>27 months old; meat grade: ≥4; steers) produced in Japan. The cattle were fattened for 20 months on a common beef cattle diet consisting of grass hay plus compound feed (wheat bran, barley, feed corn, soybean meal, rice bran, minerals, vitamins, etc.). Immediately after slaughter, the sternocleidomastoid muscle of the cattle was collected by the technical staff of the Commercial Meat Processing Center (Kobe, Japan). This experiment does not include animal experimentation because the sampling was designed in the process of commercial distribution of beef with the cooperation of general livestock farmers.

RNA Preparation from Fat Tissue
Subcutaneous fat and intramuscular fat were collected and immersed in Nucleic acid preservation (NAP) buffer (19 mM EDTA, 18 mM trisodium citrate, and 3.8 M ammonium sulfate pH 5.2). After complete removal of fine connective tissue under a stereomicroscope, each tissue was cryopreserved at −80 • C.
Total RNA was purified using~100 mg of fat tissue by the Maxwell RSC simply RNA tissue kit using a Maxwell RSC instrument (Promega K. K., Tokyo, Japan) according to manufacturer instructions. The concentration of the obtained RNA was measured using a NanoDrop spectrophotometer (Thermo Fisher Scientific K.K., Tokyo, Japan). RNA samples were quality checked by the Agilent 2200 TapeStation (Agilent Technologies Japan, Tokyo, Japan), and high-purity samples with a degradation index RNA integrity number equivalent ≥ 7.6 were subjected to an RNA-seq analysis.

Preparation of Sequence Library
Purified RNA (10 ng) was amplified by PCR amplification (7 cycles) using the Clontech SMART-Seq v4 Ultra Low Input RNA kit (Takara Bio, Kusatsu, Japan) for construction and sequencing of the total RNA-seq library. The PCR product was purified by the magnetic bead method using AMPure XP (Beckman Coulter, Tokyo, Japan). Double-stranded cDNA (0.2 ng) was generated and barcoded by PCR amplification (11 cycles) using the Nextera XT DNA Library Prep kit (Illumina K.K., Tokyo, Japan). The sequence library was validated using a fragment analyzer (Agilent Technologies, Tokyo, Japan).
The Genedata Profiler Genome (v.13.0.11; Genedata K.K., Tokyo, Japan) was used for sequence analysis. The read sequence was mapped onto the genome sequence using the mapping software STAR (v.2.6.0c) annotated from the obtained genome position, and the expression level was calculated for each gene and transcript. Data were analyzed using the bovine genome sequence ARS-UCD1.

Gene Expression Analysis
Differentially expressed genes (DEGs) were determined by comparing the difference between the log 2 transcripts per million (TPM) values of intramuscular and subcutaneous fat tissue, which was statistically evaluated with a t-test (significance: p < 0.05). Clustering analysis was performed using a 1625 DEGs MultiExperiment Viewer (Mev v.4.9, https: //sourceforge.net/projects/mev-tm4/, accessed on 12 April 2021). A Gene Ontology (GO) analysis and a Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed on 1625 DEGs using the Database for Annotation, Visualization, and Integrated Discovery (DAVID v.6.8, http://www.genome.jp/kegg/, accessed on 12 April 2021).

qPCR Analysis
cDNA libraries were constructed using ReverTra Ace qPCR RT master mix with gDNA Remover (Toyobo, Osaka, Japan). qPCR was performed using the One-step TB Green Primescript RT-PCR kit II (Takara Bio) and a Dice Real-Time System Lite thermocycler (Takara Bio) according to manufacturer instructions, with bovine β-tubulin used as an internal standard. The thermocycling conditions were as follows: heating at 95 • C for 30 s, followed by 40 to 47 cycles of denaturation at 95 • C for 15 s, annealing for 60 s at the indicated temperatures for each primer ( Figure S1), and extension at 60 • C for 30 s.
For tissue immunostaining, deparaffinized sections were activated with an antigenretrieval solution (Histo VT one solution; Nacalai Tesque, Kyoto, Japan) at 90 • C for 20 min. After permeabilization with 0.3% Triton-X in PBS (−) for 30 min and inactivation of endogenous peroxidase with 0.3% H 2 O 2 solution for 30 min, the sections were treated overnight with an antibody diluted with Can Get Signal immunostaining solution A (Toyobo). The bound antibody was visualized using Histofine Max-Po (Nichirei Biosciences, Tokyo, Japan) and ImmPACT DAB (Vector Laboratories, Burlingame, CA, USA). Counterstaining was performed with Mayer's hematoxylin solution (Fujifilm Wako Chemicals), and image analysis was performed using an all-in-one microscope system (BZ8000; Keyence, Osaka, Japan).

Statistical Analysis
The density of protein bands was measured using an image analysis software (Image J, National Institutes of Health, Bethesda, MD, USA). A multivariate data analysis of the Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) was performed using the SIMCA14 software (Inforcom, Tokyo, Japan) [15,18]. Statistical significance was determined at a p < 0.05 according to Student's t-test with JMP12 (SAS Institute Japan, Tokyo, Japan) and Excel 2019 (Microsoft Japan, Tokyo, Japan). The q-values for the GO and KEGG analyses were performed using DAVID v.6.8 according to the manual.

RNA-seq Analysis of Subcutaneous Fat and Intramuscular Fat
During the analysis of a small amount of intramuscular fat samples (like marble in muscle) collected from Japanese Black cattle, we also collected subcutaneous fat from the same cattle to compare with intramuscular fat (Figure 1a). We found no difference in the appearance of adipocytes between subcutaneous and intramuscular fat tissues, although the size of adipocytes was significantly smaller (p < 0.05) in intramuscular fat (Figure 1b After extracting expressed genes from the primary analysis data, a total of 27,60 genes were detected. Gene expression levels were normalized as TPM, and false-positiv genes were removed. To comprehensively view gene expression in each fat tissue, w performed a multivariate analysis of the data and displayed the genes as a visual plo After extracting expressed genes from the primary analysis data, a total of 27,606 genes were detected. Gene expression levels were normalized as TPM, and false-positive genes were removed. To comprehensively view gene expression in each fat tissue, we performed a multivariate analysis of the data and displayed the genes as a visual plot ( Figure 2). Actin alpha 2 (ACTA2), tropomyosin 2 (TPM2), and myosin heavy chain 11 (MYH11) are typical genes of muscle fibers, predicted as contaminants, and removed from the analytical data. ard deviation, * p < 0.05). (d) Summary of RNA-seq analysis. We analyzed the intramuscular and subcur Japanese Black cattle by RNA-seq analysis.
After extracting expressed genes from the primary analysis data, a total of 27,606 genes were detected. Gene expression levels were normalized as TPM, and false-positive genes were removed. To comprehensively view gene expression in each fat tissue, we performed a multivariate analysis of the data and displayed the genes as a visual plot ( Figure 2). Actin alpha 2 (ACTA2), tropomyosin 2 (TPM2), and myosin heavy chain 11 (MYH11) are typical genes of muscle fibers, predicted as contaminants, and removed from the analytical data. We then compared gene expression between intramuscular fat and subcutaneous fat. First, the TPM values were log-transformed and then compared between the two groups of subcutaneous fat and intramuscular fat, and the genes that obtained positive or negative significant differences by t-test were extracted as DEGs. The extraction of DEGs revealed 1625 genes expressed at significantly higher levels in intramuscular fat than in subcutaneous fat and 130 genes showing higher expression in subcutaneous fat than intramuscular fat (p < 0.05). We first investigated genes showing significantly higher expression in intramuscular fat than in subcutaneous fat. Of the 1625 DEGs in intramuscular fat, 65 were identified as highly expressed genes (count value ≥ 1000, the expression ratio of int/sub ≥ 2; p < 0.05; Figure 3a) and 25 as specifically expressed genes (count value ≥ 16, expression ratio of int/sub ≥ 4; p < 0.05; Figure 3b). Of these, we identified eight new genes, namely carboxypeptidase E (CPE), tenascin C (TNC), transgelin (TAGLN), collagen type IV alpha 5 (COL4A5), cysteine and glycine-rich protein 2 (CSRP2), PDZ and LIM domain 3 (PDLIM3), phosphatase 1 regulatory inhibitor subunit 14A (PPP1R14A), and regulator of calcineurin 2 (RCAN2), as highly and specifically expressed genes (Figure 3c). The 65 highly expressed genes and 25 specifically expressed genes are listed in Figure S2, and the top 15 highly expressed genes in intramuscular fat are shown in Figure S3.

Verification of Characteristic Genes Expressed in Intramuscular Fat
We then verified the reliability of the expression data using qPCR. We verified the eight genes identified as highly expressed in intramuscular fat, and significant differences were confirmed for six genes (Figure 4). The reproducibility of the gene expression was confirmed for each gene. 65 were identified as highly expressed genes (count value ≥ 1000, the expression ratio of int/sub ≥ 2; p < 0.05; Figure 3a) and 25 as specifically expressed genes (count value ≥ 16, expression ratio of int/sub ≥ 4; p < 0.05; Figure 3b). Of these, we identified eight new genes, namely carboxypeptidase E (CPE), tenascin C (TNC), transgelin (TAGLN), collagen type IV alpha 5 (COL4A5), cysteine and glycine-rich protein 2 (CSRP2), PDZ and LIM domain 3 (PDLIM3), phosphatase 1 regulatory inhibitor subunit 14A (PPP1R14A), and regulator of calcineurin 2 (RCAN2), as highly and specifically expressed genes (Figure 3c). The 65 highly expressed genes and 25 specifically expressed genes are listed in Figure S2, and the top 15 highly expressed genes in intramuscular fat are shown in Figure S3.

Expression Analysis of the COL4 Isoform
Intramuscular fat comprises adipocytes that ectopically proliferate and differentiate around the epimysium, perimysium, and endomysium [19] (Figure 5a). In the present study, we focused on COL4A5. COL4A has six different isoforms, the chains of which form heterotrimers to construct a complex collagen network in the basement membrane [20]. Comparison of the expression of the six isoforms revealed COL4A5 and COL4A6 as showing the highest expression in intramuscular fat relative to subcutaneous fat (Figure 5b), with this confirmed by qPCR ( Figure S4). We then confirmed COL4A isoform expression in the basement membrane around intramuscular fat by tissue immunostaining. Unlike other collagens, COL4A is not expressed in muscle tissue. We observed that COL4A isoforms were unevenly distributed in the perimysium and endomysium around the muscle fibers (Figure 5c). COL4A5 and COL4A6 were only distributed in the thickened basement membrane around adipocytes in the endomysium. By contrast, COL4A1 and COL4A2 abundantly localized in the perimysium, whereas COL4A3 and COL4A4 (data not shown) were detected at low levels in the endomysium.

Verification of Characteristic Genes Expressed in Intramuscular Fat
We then verified the reliability of the expression data using qPCR. We verified the eight genes identified as highly expressed in intramuscular fat, and significant differences were confirmed for six genes (Figure 4). The reproducibility of the gene expression was confirmed for each gene.

Expression Analysis of the COL4 Isoform
Intramuscular fat comprises adipocytes that ectopically proliferate and differentiate around the epimysium, perimysium, and endomysium [19] (Figure 5a). In the present study, we focused on COL4A5. COL4A has six different isoforms, the chains of which form heterotrimers to construct a complex collagen network in the basement membrane [20]. Comparison of the expression of the six isoforms revealed COL4A5 and COL4A6 as showing the highest expression in intramuscular fat relative to subcutaneous fat ( Figure  5b), with this confirmed by qPCR ( Figure S4). We then confirmed COL4A isoform expression in the basement membrane around intramuscular fat by tissue immunostaining. Unlike other collagens, COL4A is not expressed in muscle tissue. We observed that COL4A isoforms were unevenly distributed in the perimysium and endomysium around the muscle fibers (Figure 5c). COL4A5 and COL4A6 were only distributed in the thickened basement membrane around adipocytes in the endomysium. By contrast, COL4A1 and COL4A2 abundantly localized in the perimysium, whereas COL4A3 and COL4A4 (data not shown) were detected at low levels in the endomysium.  Data represent the mean ± standard error (n = 3). ** p < 0.01, (c) Tissue immunostaining for COL4A family proteins in paraffin sections of muscle tissue from Japanese Black cattle. DAB staining (brown) was performed using primary antibodies targeting specific isoforms, with hematoxylin used for counterstaining. White arrowheads indicate perimysium, and gray arrowheads indicate endomysium. Size bar, 100 µm. Int, intramuscular adipocytes; M, muscle tissue.

Pathway Analysis of Differentially Expressed Genes
We then evaluated the intracellular functions of DEGs using a KEGG enrichment analysis and a GO analysis. KEGG results showed the characteristic pathways for intramuscular and subcutaneous fat ( Figure 6). DEGs in intramuscular fat were significantly associated with 20 KEGG pathways, including focal adhesions (34 genes; p = 0.0004), pathways in cancer (50 genes; p = 0.0065), PI3K-Akt signaling (42 genes; p = 0.0231), MAPK signaling (32 genes; p = 0.0316), Ras signaling (29 genes; p = 0.0493), and transforming growth factor (TGF)-β signaling (13 genes; p = 0.0518). By contrast, DEGs in subcutaneous fat were associated with cellular metabolism pathways, including metabolic pathways (38 genes; p = 0.0001), lysosomes (10 genes; p = 0.0002), and the regulation of adipocyte lipolysis (7 genes; p = 0.0002). To further narrow down the genes involved in intramuscular fat formation, we analyzed the cellular functions of these genes using GO (Figure 7). A total of 1625 DEGs were significantly associated with ECM binding (3.8%; q-value, 2.60) and zinc ion binding (10.8%; q-value, 2.05) as a molecular function. Intramuscular fat was also associated with circulatory system development (24.6%; q-value, 5.12), actin cytoskeleton organization (16.9%; q-value, 3.44), cell surface receptor signaling (30.7%; q-value, 3.19), and cell migration (20.0%; q-value, 2.78) among biological processes. GO analysis revealed that the DEGs of intramuscular fat were associated with ECM binding, cell-cell adhesion, angioplasty, and receptor signal. These findings suggest that the surrounding environment influences intramuscular fat accumulation. The intracellular function of intramuscular adipocytes resembles that of infiltrating cells and tumors in response to the microenvironment [21]. This unique cellular environment is seemingly involved in intramuscular fat accumulation.  To further narrow down the genes involved in intramuscular fat formation, we analyzed the cellular functions of these genes using GO (Figure 7). A total of 1625 DEGs were significantly associated with ECM binding (3.8%; q-value, 2.60) and zinc ion binding (10.8%; q-value, 2.05) as a molecular function. Intramuscular fat was also associated with circulatory system development (24.6%; q-value, 5.12), actin cytoskeleton organization (16.9%; q-value, 3.44), cell surface receptor signaling (30.7%; q-value, 3.19), and cell migration (20.0%; q-value, 2.78) among biological processes. GO analysis revealed that the DEGs of intramuscular fat were associated with ECM binding, cell-cell adhesion, angioplasty, and receptor signal. These findings suggest that the surrounding environment influences

Analysis of COL4-Related Cellular Signals in Intramuscular Fat
We hypothesized that the ECM of the basement membrane contributes to intramuscular fat formation. TGF-β, integrins, and cell migration factors represent ECM-related signals (Table 1). TGF-β is a growth factor that induces COL4A secretion in connective tissue [22]. TGFB2 and TGFB3 were highly expressed in intramuscular fat tissue compared to subcutaneous fat tissue, whereas TGFB1 was expressed at a lower level. SMAD, a major transcription factor that functions in the canonical pathway of TGF-β, showed similar expression levels between subcutaneous and intramuscular fat tissues. However, the expression levels of SMAD7, CHRD, and TGIF1, which are regulators of TGF-β signaling, are significantly differed between tissues. Furthermore, heterodimers, namely ITGA2, ITGB1, and ITGB6, were significantly highly expressed in intramuscular fat compared to subcutaneous fat tissue. The expression of heterodimers, integrin α (ITGA) and integrin β (ITGB), ITGA2, ITGB1, and ITGB6, was significantly higher in intramuscular compared to subcutaneous fat tissue. By contrast, ITGA6, which encodes the major ITGA in adipocytes [23], was highly expressed in subcutaneous fat tissue. Concerning cell migration factors, CXCL5, FGF9, and FGF13 showed significantly higher expression in intramuscular fat relative to subcutaneous fat, whereas CXCR4, which is a receptor for CXCL5, showed similar expression levels between subcutaneous and intramuscular fat tissues.

Analysis of COL4-Related Cellular Signals in Intramuscular Fat
We hypothesized that the ECM of the basement membrane contributes to intramuscular fat formation. TGF-β, integrins, and cell migration factors represent ECM-related signals (Table 1). TGF-β is a growth factor that induces COL4A secretion in connective tissue [22]. TGFB2 and TGFB3 were highly expressed in intramuscular fat tissue compared to subcutaneous fat tissue, whereas TGFB1 was expressed at a lower level. SMAD, a major transcription factor that functions in the canonical pathway of TGF-β, showed similar expression levels between subcutaneous and intramuscular fat tissues. However, the expression levels of SMAD7, CHRD, and TGIF1, which are regulators of TGF-β signaling, are significantly differed between tissues. Furthermore, heterodimers, namely ITGA2, ITGB1, and ITGB6, were significantly highly expressed in intramuscular fat compared to subcutaneous fat tissue. The expression of heterodimers, integrin α (ITGA) and integrin β (ITGB), ITGA2, ITGB1, and ITGB6, was significantly higher in intramuscular compared to subcutaneous fat tissue. By contrast, ITGA6, which encodes the major ITGA in adipocytes [23], was highly expressed in subcutaneous fat tissue. Concerning cell migration factors, CXCL5, FGF9, and FGF13 showed significantly higher expression in intramuscular fat relative to subcutaneous fat, whereas CXCR4, which is a receptor for CXCL5, showed similar expression levels between subcutaneous and intramuscular fat tissues. We then evaluated the Ras signal pathway involving Ras protein, which at acts as a molecular switch to control various signal transduction pathways. GDP/GTP-exchange factor (GEF), GTPase activating factor (GAP), and effectors are listed as regulators of small GTPases [24]. Among GEFs, RASGRP3 and ARHGEF26 were highly expressed in intramuscular fat relative to subcutaneous fat. For GAP, ARHGAP10 and ARHGEF24 were highly expressed in intramuscular fat, as were RASSF9 and RASSF10, which are effector molecules for Ras GTPase.

Discussion
Japanese Black cattle are not crossed with other cattle breeds, and their pedigree is strictly controlled; therefore, the use of Japanese Black cattle in research introduces low genetic variation. We previously analyzed the metabolite [15] and lipid composition [18] of intramuscular fat in Japanese Black cattle. In this study, we performed RNA-seq analysis to explore the key genes involved in the intramuscular fat accumulation of Japanese Black cattle. Previous studies indicated that comparisons between different cattle breeds tend to detect DEGs related to popular adipocyte functions, such as lipid metabolism and adipogenesis [12]. In the present study, we focused on intramuscular fat formation-related genes. Therefore, to avoid bias, we collected intramuscular and subcutaneous fat tissues from the same cattle.
Microscopic observation showed little contamination of other tissues (Figure 1b). Intramuscular fat is known to comprise a smaller portion of adipocytes than subcutaneous fat [25]. We confirmed that the intramuscular fat used in this analysis comprised a smaller portion of adipocytes than subcutaneous fat (Figure 1c). RNA-seq analysis was performed on a small sample (100 mg) to avoid contamination with other connective tissues (Figure 1d). A higher number of DEGs was detected in intramuscular than in subcutaneous fat tissues (Figure 2). The muscle-derived myosin and actin genes were detected by RNA-seq. These genes might be derived from a fragment of muscle tissue attached to the sample or intramuscular fat-specific cytoskeletal proteins [16]. To render our search more comprehensive, we explored characteristic genes of intramuscular fat using two search criteria for genes, high expression and specificity. Of the identified genes, we have identified eight novel genes highly and specifically expressed in intramuscular fat tissues from Japanese Black cattle ( Figure 3). These eight genes were validated as highly expressed in intramuscular fat tissue by quantitative PCR assays ( Figure 4).
The role of these eight genes in intramuscular fat formation was inferred from their respective functions. TNC is a high molecular weight glycoprotein that plays an ECM to stabilize the basement membrane with collagen. Interestingly, TNC is reportedly induced by obesity-related inflammatory responses in the adipose tissue [26] and tumor microenvironment [27]. Tenascin has four isoforms, TNC, TNR, TNXB, and TNW, in cattle. In our RNA-seq data, only TNC was highly detected in intramuscular fat. TAGLN is a TGF-β inducible gene that functions as an actin cross-linking protein [28]. TAGLN also promotes cell migration and adipocytic differentiation of human bone marrow-derived stromal stem cells [29]. CPE acts as an exopeptidase and is involved in the biosynthesis of neuropeptides and peptide hormones in endocrine tissues and the nervous system. Deficiency in CPE activity has been linked to diabetes, obesity, and reduced learning ability [30]. Genetic analyses highlighted that single nucleotide polymorphisms in CPE are associated with beef quality [31]. Of these eight genes, CPE, TNC, and TAGLN have been detected in the intramuscular fat of other cattle [32]. These three genes are presumed to be essentially expressed in intramuscular fat regardless of cattle pedigree.
CSRP2 contains three isoforms (CSRP1, CSRP2, CSRP3). These CSRPs act as transcription factors in the nucleus and promote cell differentiation in various cell types. CSRP2 expression is reportedly induced by TGF-β [33], an essential gene for epithelialmesenchymal conversion [28]. Although CSRP1, CSRP2, and CSRP3 were highly expressed in intramuscular fat tissue ( Figure S2), their detailed function is unknown. RCAN2 is involved in muscle differentiation and cancer progression [34]. PPP1R14A is a typically upregulated gene in obesity-related hypertension [35]. PDLIM3 are sarcomere-related proteins. Of note, PDLIM3 regulates cell-cell adhesion and migration and promotes malignant metastasis [36]. Unlike other collagens, COL4A is a unique collagen isoform in the basement membrane [37]. COL4A5 forms supramolecular networks with COL4A [20].
Tissue immunostaining showed that COL4A5 colocalized with COL4A6 in the thickened basement membrane of the endomysium (Figure 4). Based on their expression patterns (Figure 5b), COL4A1 and COL4A2 seemed to contribute to both subcutaneous and intramuscular fat accumulation, whereas COL4A5 may be involved in the "microintramuscular fat" (Figure 5a). Recently, the specific expression of COL4A5 in the tumor microenvironment has been reported [38]. COL4A5 expression in the basement membrane may involve the infiltration and proliferation of ectopic cells, such as intramuscular fat and tumors.
Next, we investigated the intracellular functions of DEGs in each fat tissue. KEGG analyses revealed that intramuscular fat is associated with gene pathways related to cell migration and proliferation ( Figure 6). Notably, PI3K-Akt, MAPK, and Ras signaling pathway upregulation suggest the activation of downstream membrane receptors and ECM binding in intramuscular fat tissues. The cellular function associated with intramuscular fat was consistent with the GO analysis results (Figure 7).
From these results, we focused on the ECM and compared the TPM values of genes related to TGF-β, cell adhesion, and cell migration factors (Table 1). TGF-β regulates gene expression by SMAD-mediated signal transduction [39]. TGF-β reportedly induces CPE, TNC, TAGLN, CSRP2, and COL4A5 expression [28,33,40]. Therefore, it is suggested that TGF-β/SMAD pathway-mediated promotion of gene transcription might underly intramuscular fat formation in Japanese Black cattle. Furthermore, ITGA6 and ITGB1, the subunits of integrin α6β1 that specifically bind to the COL4A5 heterotrimer [20], were increased in intramuscular fat. The expression of cell migration factors CXCL5, FGF9, and FGF13, was also increased in intramuscular fat tissue (Table 1). Therefore, in the intramuscular fat tissue, increased cell adhesion to the basement membrane may promote preadipocyte cell migration by increasing ITGB1, CXCL5, FGF9, and FGF13 expression.
Surprisingly, cancer-related pathways were significantly associated with intramuscular fat tissue ( Figure 6), including TGF-β, Wnt-1, Ras, and Rho signaling-related genes ( Figure S5). The TGF-β is known as a non-canonical pathway. This pathway activates Ras/MAPK signaling and the Rho family signaling of small GTPase in a SMAD-independent manner [40]. These small GTPase-mediated TGF-β signaling reportedly promote cell migration and differentiation in various inflammatory tissues and tumors [39]. These findings suggest that the cell function specific to intramuscular fat accumulation might be mediated by TGF-β. Therefore, TGF-β is considered a pivotal gene in intramuscular fat formation in Japanese Black cattle. Comparison between two groups by DEGseq2 is shown in ( Figure S6).
The TPM levels did not differ between intramuscular and subcutaneous fat tissues in the Ras and Rho family of small GTPase. Many previous studies have shown that small GTPase activity is specifically regulated during signal transduction [24,41]. Therefore, we compared the expression of key regulators between intramuscular and subcutaneous fat tissues. RNA-seq analyses detected 66 genes as regulators of the Ras or Rho family of small GTPases. Of these, RASGRP3 and ARHGEF26 were significantly more highly expressed in intramuscular fat tissue than subcutaneous fat tissue (Table 1). RASGRP3 activates the Ras signals and promotes invasion and cell migration [42]. ARHGEF26 activates the Rho family of RhoG and also regulates cancer cell migration [43]. Besides, ARHGAP10, ARHGAP24 [44], and DLC3, which act as GAPs for the Rho family [45], were highly expressed in intramuscular fat. Moreover, ARHGEF33 expression was significantly lower in intramuscular fat tissue than subcutaneous fat tissue. Differential expression of GEFs and GAPs has been widely reported in tumors. These regulators are reportedly sufficient to regulate small GTPase activities by altering their expression levels [46]. These regulators may also regulate the small GTPase activities required for the proliferation, differentiation and maintenance of adipocytes around the basement membrane of the endomysium that compose the intramuscular fat.

Conclusions
In this study, our RNA-seq analysis led to identifying eight genes highly and specifically expressed in the intramuscular fat of Japanese Black cattle. Among them, COL4A5 was found to be highly expressed in the basement membrane around intramuscular fat. Pathway analysis also revealed that cell adhesion, cell proliferation, and cancer pathways are involved in intramuscular fat formation. In addition, the relationship between TGF-β signaling and the regulators of Ras and Rho families are seemingly involved in intramuscular fat accumulation. These genes might be used as molecular markers to better understand intramuscular fat formation [15].

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/genes12081107/s1. Figure S1: List of primer sequences used for qPCR. Figure S2: List of 65 highly expressed genes and 25 genes specifically expressed in intramuscular fat. Figure S3 List of top 15 genes showing elevated expression in intramuscular fat and subcutaneous fat. Figure S4: Expression of COL4A isoforms by qPCR. Figure S5: List of pathway analyses. Figure S6: Comparison between two groups by DEGseq2 package.

Data Availability Statement:
The data presented in this study are contained in the article.