Comparative Transcriptome Analysis of Adipose Tissues Reveals that ECM-Receptor Interaction Is Involved in the Depot-Specific Adipogenesis in Cattle

Adipocytes mainly function as energy storage and endocrine cells. Adipose tissues showed the biological and genetic difference based on their depots. The difference of adipocytes between depots might be influenced by the inherent genetic programing for adipogenesis. We used RNA-seq technique to investigate the transcriptomes in 3 adipose tissues of omental (O), subcutaneous (S) and intramuscular (I) fats in cattle. Sequence reads were obtained from Illumina HiSeq2000 and mapped to the bovine genome using Tophat2. Differentially expressed genes (DEG) between adipose tissues were detected by EdgeR. We identified 5797, 2156, and 5455 DEGs in the comparison between OI, OS, and IS respectively and also found 5657 DEGs in the comparison between the intramuscular and the combined omental and subcutaneous fats (C) (FDR<0.01). Depot specifically up- and down- regulated DEGs were 853 in S, 48 in I, and 979 in O. The numbers of DEGs and functional annotation studies suggested that I had the different genetic profile compared to other two adipose tissues. In I, DEGs involved in the developmental process (eg. EGR2, FAS, and KLF7) were up-regulated and those in the immune system process were down-regulated. Many DEGs from the adipose tissues were enriched in the various GO terms of developmental process and KEGG pathway analysis showed that the ECM-receptor interaction was one of commonly enriched pathways in all of the 3 adipose tissues and also functioned as a sub-pathway of other enriched pathways. However, genes involved in the ECM-receptor interaction were differentially regulated depending on the depots. Collagens, main ECM constituents, were significantly up-regulated in S and integrins, transmembrane receptors, were up-regulated in I. Different laminins were up-regulated in the different depots. This comparative transcriptome analysis of three adipose tissues suggested that the interactions between ECM components and transmembrane receptors of fat cells depend on the depot specific adipogenesis.


Introduction
Fat is a loose connective tissue composed of adipocytes and functions as not only the main site of energy storage but also the highly active metabolic and endocrine organ. It is restricted to certain areas in the body such as subcutaneous layer between muscle and dermis, around internal organs and between muscles. The amount and distribution of fat are important factors that influence the meat quality in the beef industry [1,2]. For several decades, it has been known that adipose tissues are regionally heterogeneous with respect to metabolic activities and functions [3][4][5]. Evidences showed that the numerous biological and genetic differences among adipose tissues depend on their anatomical locations [6][7][8]. Adipocytes isolated from different depots differ in size, lipoprotein lipase release, lipid synthetic capacity, fatty acid incorporation and other characteristics [9][10][11][12][13]. Regional variations have been also observed in the replicative potential, fatty acid transfer and adipogenic differentiation of preadipocytes originating from various depots of the same individuals under the identical in vitro conditions [14][15][16][17][18]. The inter-depot variations might have the influences of the inherent genetic programming of the adipocyte development and the microenvironment of adipose depots such as their hormonal response, local nutrient availability, innervations and anatomic constraints [14][15][16][17][18]. The fat deposition in muscle, intramuscular fat, might be regulated by the genetic mechanisms different from that in other adipose depots for the occurrence and development. Intramuscular fat content in meat is one of important traits that influence eating quality such as meat tenderness, juiciness, and taste [19]. However, very little is known on the genetic mechanism of intramuscular fat because it is the most difficult depot to study developmentally. Especially, there is the lack of information concerning the transcriptomic differences among omental, subcutaneous and intramuscular adipose tissues.
Gene expression profiling or transcriptome analysis provides new insights to understand the molecular basis of adipogenesis [20]. A few of transcriptome studies in adipose tissues has been reported in cattle using microarray and 39 digital gene expressiontag analysis [21,22] The application of next-generation sequencing (NGS) to this area have started to reveal the complexities of different biological processes at an unprecedented level of sensitivity and accuracy [20]. In the present study, we investigated the difference of the depot specific gene expression from 3 different adipose tissues of omental, subcutaneous and intramuscular tissues in cattle using RNA-seq technology.

Ethics Statement
All animal care and experimental procedures were reviewed and approved by the Institutional Animal Care and Use Committee of the National Institute of Animal Science (No. 2010-042).

Animals
Nine heads of Hanwoo (Korean native cattle, Bos taurus coreanae) cattle were fed and managed in the feeding barn at Hanwoo Experimental Station in National Institute of Animal Science under the high quality beef production program (2007). The steers used for this study were 30 months old at slaughter (Body Weight = 835615.2 kg).

Sample Preparation
Immediately, after stunning and exsanguination, the muscle and fat portions between the 6th to 7th ribs were removed, and the subcutaneous and intramuscular fat depots were sampled from this rib section. The omental adipose tissue was taken within the lesser curvature of the abomasums and stored at 280uC freezer until the analysis.

RNA-seq Library Preparation and Sequencing Analysis
Total RNAs of subcutaneous, intramuscular, and omental fat tissues were isolated using TRIzol (Invitrogen) and a RNeasy RNA purification kit with DNase treatment (Qiagen). mRNA was isolated from the total RNA using oligo-dT beads and were reverse transcribed into double strand cDNA fragments. Constructing and sequencing RNA-seq library for each sample were carried out based on protocols of Illumina HiSeq2000 to generate 90 pair-end reads. Quality of RNA-seq reads from all of the adipose tissues was checked using FastQC ( Figure S1). The reads passed the quality control were mapped to Bovine Taurus genome (bosTau6) from UCSC using Tophat2 (v2.0.2) and were counted using HTseq (v0.5.3p3). The RNA sequencing data from this study have been submitted to the NCBI Gene Expression Omnibus (GEO) under the accession number of GSE39618 (http://www.ncbi.nlm.nih. gov/geo/).

Statistical Analysis and DEG Identification
EdgeR was used to identify differentially expressed genes between adipose tissues [23]. This package is designed for the analysis of replicated count-based expression data and is based on a negative binomial model. RNAseq data may exhibit more variability than expected in a poisson distribution because it is more widely dispersed in the genome. When a negative binomial model is used, the dispersion has to be estimated before the DEG analysis is carried out. EdgeR provided Cox-Reid profile-adjusted likelihood method to estimate dispersion for the pairwise comparisons between Omental and intramuscular (OI), Omental and subcutaneous (OS), intramuscular and subcutaneous (IS), and intramuscular and the combined omental and subcutaneous (IC). After negative binomial models are fitted and dispersion estimate are obtained, a generalized linear model (GLM) likelihood ratio test was used to determine differential expression. The GML likelihood ratio test automatically takes all known sources of variation into account. Significant DEGs each pairwise comparison were selected at FDR,0.01. In order to investigate differentially regulated genes specific to each depot, we selected common DEGs from two-way pairwise comparisons with one of tissues in common. For example, we compared DEGs with the same direction of log FC between OI and OS and selected common DEG to identify the omental-specifically regulated genes. Same methods were applied to detect DEGs specific to the intramuscular and subcutaneous fats.

Functional Annotation of DEGs
The bovine Ensembl gene IDs were converted to official gene symbols by cross-matching to human Ensembl gene IDs and official gene symbols. The official gene symbols of human homologues of bovine genes were used for functional clustering and enrichment analyses using the Database for Annotation, Visualization and Integrated Discovery (DAVID) [24]. The representation of functional groups in omental, intramuscular, and subcutaneous fat tissue relative to the whole genome was investigated using the Expression Analysis Systematic Explorer (EASE) tool [25] within DAVID, which is a modified Fisher's exact test to measure enrichment of gene ontology (GO) terms [26]. To identify enrichemented GO terms functionally clustered genes were filtered by EASE value ,0.01 and selected.

Reverse Transcription and Quantitative Real-time PCR
The 1.5 ug of extracted total RNA was used in a final volume of 20 ml to synthesize cDNA using a High Capacity RNA-to-cDNA Kit (Applied Biosystems, Foster City, CA). With the synthesized cDNA sample of each tissue, quantitative real-time PCR was performed to analyze the quantity of mRNA that coded for the selected genes using the ABI 7300 Real-Time PCR system (Applied Biosystems, Foster City, CA) and a quantitative real-time PCR kit (DyNAmo HS SYBR Green qPCR kit, Finnzymes, Finland). Amplification was carried out for 1 cycle at 95uC for 5 min; 40 cycles at 95uC for 15 s, annealing temperature listed on the Table S1 for 30 s, 72uC for 30 s; 1 cycle of 72uC for 10 min. For amplification, 1 pM of the primers were used for the amplification. The mRNA quantities of the target genes were normalized to that of a reference gene. The qRT-PCR for all target genes was performed twice with two different reference genes, b-actin (ACTB) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH), to make the result more reliable. Using threshold cycle (CT) values for these genes, relative expression levels were calculated with the 2 2DD Ct method [27].

Quality of RNA-sequence Reads in Adipose Tissues
We obtained RNA-seq reads from three different adipose tissues (subcutaneous, intramuscular, and omental) of nine Hanwoo animals [GenBank:GSE39618]. The numbers of total sequence reads and mapping rates for each sample were summarized in Table S2. Only a few amount of sequence reads (,0.05%) did not pass the quality filtering. The average numbers of sequence reads in each tissue were 38, 36, and 35 M reads in subcutaneous (S), intramuscular (I), and omental (O) fat respectively. Among the sequence reads passed the quality control, on average, 96.9% reads in subcutaneous fat, 98.3% in intramuscular, and 97.5% in omental fat were successfully mapped to bovine genome (bosTau6) using TopHat (v2.0.2).

Identification of Differentially Expressed Genes (DEGs)
To investigate the difference in expression profiles of genes among adipose tissues of subcutaneous (S), intramuscular (I), and omental (O), we identified differentially expressed genes (DEGs) in the pairwise comparison of three adipose depots [28] (Figure 1 Table 1, File S1). The number of DEG in OI and IS were two times higher compared with the OS. We also compared between the intramuscular fat (I) and the combined omental and subcutaneous fat (C) and found 5657 DEGs (2992 up-regulated in I and 2665 in C) ( Table 1, File S1), which showed the similar pattern when the intramuscular fat was compared with other adipose fats separately. In order to identify depot specific DEG, we selected common DEGs from the two-way pairwise comparisons with one of tissues in common. We found 853 (302 up-regulated and 551 down-regulated), 4276 (2507 up-regulated and 1769 down-regulated), and 979 (621 up-regulated and 1320 downregulated) DEGs in subcutaneous (S), intramuscular (I), and omental (O) respectively ( Table 1). The number of DEGs in the intramuscular was 4 times higher than that in other tissues.

Validation of DEGs Using qRT-PCR
We performed qRT-PCR to technically validate the DEGs detected in each adipose tissue. A total of 30 genes consisting of ten genes (5 up-regulated and 5 down-regulated) in each tissue were randomly selected with logFC.2 (Table S1). Correlation was calculated to compare the expression levels between RNAseq and qRT-PCR [27]. Fold changes of DEGs between two techniques were significantly correlated in three depots (I: Pearson's r 2 = 0.0,85, O: Pearson's r = 0.96, and S: Pearson's r = 0.95) Table S3, and Table S4). The results confirmed that DEGs identified in this study were very reliable.

DEGs Involved in the Adipogenesis and Lipid Metabolism
Literatures have reported that many genes were involved in the adipogenesis and lipid metabolism. Among the depot-specific DEGs identified, we found 28 genes related to the mechanisms (Table 2). CEBPG, DDIT3, KLF7, and EBF1 that were known to function as anti-adipogenic factor increased in intramuscular fat. UCP1, FOXC2, and PPARGC1A related to brown adipogenesis were detected in the omental and intramuscular depots. AIPOQ, FABP4, and FAS that were involved in the fatty acid synthesis and metabolism were identified in intramuscular fat. FASN, LPL, and DGAT1 associated with the determinants of the deposition of intramuscular fat were down-regulated in the intramuscular depot.

Gene Ontology and Functional Annotation of DEGs
Biological process gene ontology of DEGs from the pairwise comparison was summarized in Table S5. In intramuscular fat,  developmental process was the most significantly enriched term compared to other two tissues (OI: p.value = 1.51e-47 and IS: p.value = 2.30e-49). Also developmental process was significantly enriched in the both omental (p.value = 9.47e-16) and subcutaneous fat (6.23e-12) when they were compared each other. Gene ontology (GO) analysis was also performed with the depot specific DEGs. Biological process GO terms of each adipose tissue were summarized in Figure 3a. Most significantly enriched 3 GO terms of DEGs up-regulated in the intramuscular fat were developmental process (p.value = 2.30e-40), multicellular organismal process (p.value = 1.40e-28), and biological regulation (p.value = 1.40e-10).
Immune system process (p.value = 2.90e-11) and developmental process (p.value = 6.20e-05) were the most significantly enriched GO terms of DEGs up-regulated genes in the omental and subcutaneous fat respectively For the DEGs down-regulated in each depot, developmental process (S: p.value = 1.80e-17, O: p.value = 5.80e-12) and multicellular organismal process (S: p.value = 1.80e-14 O: p.value = 6.90e-11) were the most significantly enriched in both subcutaneous and omental fat. In the intramuscular fat, immune system process (p.value = 1.40e-12) was the most significant GO terms of biological process. When intramuscular fat was compared with the combined omental and subcutaneous fat, developmental process (p.value = 1.28e-42) and immune system process (p.value = 2.33e-18) were most significantly enriched in intramuscular and the combined respectively ( Figure 3b). GO analysis of cellular components and molecular function were summarized in Table S6. Overall cellular component GO terms of depot-specific DEGs were related to the extracellular region and cell part, and molecular function GO terms were involved in the various molecule activity.

KEGG Pathway Analysis of Depot Specific DEGs
KEGG pathway analysis of DEG resulted from the pairwise comparison showed that DEGs were largely clustered into 6 representative terms (higher level of the pathways) including metabolism, cellular processes, environmental information processing, human diseases, organismal system, and genetic information processing (Figure 4, and Table S7). Many DEGs in intramuscular fat were highly enriched in the pathways of environmental information processing. When the omental fat was compared to intramuscular fat (OI), the top significant pathways were valine, leucine and isoleucine degradation (p. value = 3.39E-11), PPAR signaling pathway (p.value = 8.91E-10), lysosome (p.value = 2.59E-09) and fatty acid metabolism (p.value = 7.38E-09) (Table S7 (a)). When the subcutaneous fat was compared to intramuscular fat (SI), similar pathways were observed with less significant p-value (Table S7 (b)). The most significant pathway in both comparisons is valine, leucin and isoleucine degradation which is an amino acid metabolism. KEGG pathways enriched in the intramuscular fat showed the similar pattern of the pathways when it was compared to omental and subcutaneous fats (Table S7 ( (Table S7, (e) and (f)). When intramuscular fat (I) was compared to the combined omental and subcutaneous fats (C), same KEGG pathways were enriched in the intramuscular fat and in the combined fats as the intramuscular fat was compared to other fats separately (Table S7, (g) and (h)). KEGG pathway analysis of depot specific DEG identified 4 significantly enriched pathways of up-regulated genes in omental fat, 16 pathways in intramuscular fat, and 2 in subcutaneous fat tissue (p,0.01) (Table S8). In omental fat, DEGs related to lipid metabolisms were highly up-regulated, while those related to signaling modules and interactions were down-regulated. In intramuscular fat, DEGs involved in the human cardiac muscle disease such as cardiomyopathy (HCM, ARVC, and DCM), cancer, and Type II diabetes mellitus were most significantly upregulated but those related to the lipid metabolism and immune response were down regulated. In subcutaneous fat, genes in the pathways of host defense were up-regulated and those involved in the disease and cancer were down -regulated.

Developmental Process and ECM-receptor Interaction
Significant numbers of DEGs from the adipose tissues were enriched in various GO terms of the developmental process and their KEGG pathway analysis showed that several pathways including, axon guidance, pathways in cancer, and ECM-receptor interaction were common in all three fat tissues (Table 3). Among those pathways, interestingly, the ECM-receptor interaction also involved in other enriched pathways such as pathways in cancer, dilated cardiomyopathy, hypertrophic cardiomyopathy, and Arrhythmogenic right ventricular cardiomyopathy, which were the most significant in intramuscular fat. Although ECM-receptor interaction was common pathway in the three adipose tissues, genes involved in the pathways were differentially regulated ( Figure 5). Among genes in ECM-receptor interaction, 9, 4, and 18 genes were significantly up-regulated in subcutaneous, omental and intramuscular fat respectively ( Figure 5). Collagen genes (COL1A1, COL1A2, COL3A1, COL5A2, and COL6A3) were significantly up-regulated in subcutaneous fat and integrins (ITGA1, 2, 7, 8, 11, V, and ITGB1) were significantly up-regulated

Discussion
The amount and distribution of fat are known to be important factors to influence the meat quality in the beef industry. Especially, intramuscular fat (marbling) plays important roles in eating quality and composition of meat compared to other fats [1,2]. Some literatures reported that intramuscular fat is different from other fats in three respects: adipocyte size, metabolic activities, and developmental timing [29][30][31]. Intramuscular adipose tissue is a later developing tissue and behaves different from other adipose tissues in development of cellularity and metabolic capacity [31]. Intramuscular adipocytes are smaller than non-muscular adipocytes in cattle and pigs [29,32]. Metabolic activities including lipogenesis and lipolysis in intramuscular adipocytes were lower than other adipocytes by proteomic analysis, gene expression and/or enzyme activity [29,30]. Genetic mechanisms in intramuscular fat might be different from those in other fat depots. However little has been known about the genetic and functional difference of the intramuscular fat compared to others. In order to investigate the genetic profiles of fat tissues based on depots and understand the difference of the genetic mechanisms, we performed the transcriptome analysis from 3 different adipose tissues of omental, subcutaneous, and intramuscular fat using Illumina HiSeq 2000.

Sequencing Quality and Coverage
We obtained RNA-seq transcriptome data from omental, subcutaneous and intramuscular fats of nine cattle. The transcriptome data generated from the NGS provided very reliable sequences reads because over 99% of the sequence reads passed the QC control and the mapping rate of the sequence reads to bovine genome was very high. Our data also showed that the sequencing depth is deep enough for this study. We obtained more than 30 M reads from most of samples. Although it might not be meaningful to mention the sequencing coverage of transcriptome, we would like to have an idea how well the depth of the sequence could cover the adipose transcriptome. There are several studies showing that about 30 M reads is sufficient to detect .90% of annotated genes in yeast, chicken, and cattle [33][34][35], which suggested that our data had enough coverage, even though the tissues used for the study were different.

Different Genetic Profile in Intramuscular Fat
Identification of DEGs among the fat tissues from the 3 different adipose depots supported that the intramuscular fat had more different genetic profiles than both the omental and subcutaneous fats did. The numbers of DEG showed two times higher when the intramuscular fat was compared to other adipose tissues than when omental fat and subcutaneous fat were compared each other (Table 1). When the intramuscular fat was compared with the combined omental and subcutaneous fat, similar numbers of DEG were identified. This also confirmed that omental and subcutaneous fat shared many adipogenetic genes and pathways compared to that in intramuscular depot. The functional annotation analysis also presented the similar GO biological process terms and KEGG pathways in the omental and subcutaneous fats, but not in intramuscular fat (Table S5, Table S7, Table S8, Figure 3 and

Lower Metabolic Activities of Intramuscular Fat
One of the significantly enriched pathways of up-regulated genes in the omental fat was PPAR signaling pathway. Recent studies proved that PPARs functions as a key regulator of adipocyte differentiation, accumulation and phenotypes [36]. Genes enriched in the pathway are involved in lipid metabolisms including fatty acid transport (LPL, ACSL1, and OLR1), fatty acid oxidation (EHHADH, CPT1B, ACADL, and ACADM), adaptive thermogenesis (UCP1), and gluconeogenesis (PCK1) [37][38][39][40][41][42][43][44][45][46][47][48][49][50]. Most of them were also involved in the fatty acid metabolism pathway. PPARs (peroxisome proliferators-activated receptors) are known to play a key role in the metabolic pathways involving fatty acid oxidation and lipid metabolism [51]. We found 28 genes known to be involved in the adipogenesis and lipid metabolism in the previous studies (Table 2). Based on the expression patterns, many of them were significantly down-regulated in the intramuscular fat compared to other tissues. PPARG, a master regulator of adipogenesis, was down-regulated and other anti-adipogenic genes (CEBPG, DDIT3, GATA2, and KLF7) were up-regulated in intramuscular. Our results in part supported that intramuscular adipocytes had the low metabolic activities related to lipid metabolism based on transcriptome analysis with more various range of metabolic pathways and genes. Some genes involved in intramuscular deposition (FASN, FABP4, LPL, THRSP, and DGAT1) were also down-regulated in intramuscular fat. They were highly expressed in the early growth of cattle with high starch feeding and showed significant positive correlation with intramuscular fat content [52][53][54]. It seems that the expression results of those genes in this study were different from the previous. However it is difficult to simply compare those results because our study doesn't have information of intramuscular fat deposition. Further study might need to clear out if the difference of gene expression was caused by the age difference because animals used in this study are over 30 month old or if the expression of the genes in intramuscular fat was relatively lower because our analysis was based on the direct comparisons between different fat tissues. In addition, it might be more efficient to identify the functional effect of intramuscular fat deposition on meat quality with different tissue sampling of intramuscular fats which was grouped based on the meat quality indices such as intramuscular fat contents and water contents.

Pathways of Interested in each Adipose Depot
There were depot specifically enriched pathways in each adipose depot (Table S8). In omental fat, several members of CLDN family were differentially expressed in the pathways of leukocyte transendotheilial migration (CLDN1, 3, and 15), and cell adhesion molecules (CLDN5, and 10). The claudin (CLDN) family plays a major role in tight junction-specific obliteration of the intercellular space through calcium-independent dell-adhesion activity. CLDN6 is known to be an important regulator in adipogenesis and fat deposition but is not detected in this study [55]. Genetic manipulation of CLDN1 expression induced changes in cellular phenotype [56]. This family might affect the depotspecific morphology of adipocytes by regulating tight junctions between cells. DEGs down-regulated in the subcutaneous fat were significantly enriched in pathway of melanogenesis. A recent study demonstrated that melanic biosynthesis pathway was functional in adipose tissues, and melanin and other proteins in the melanogenic pathway of adipocytes were expressed in higher level in the obese [57]. WNT signaling in the melanogenesis pathway is an important regulator for adipogenesis or insulin secretion. Expression of WNT5B gene was increased at an early phase of adipocyte differentiation in mouse and involved in the pathogenesis of the type 2 diabetes through the regulation of adipocyte function [58]. It was interesting that the cardiovascular disease related pathways (dilated cardiomyopathy, hypertrophic cardiomyopathy, and Arrhythmogenic right ventricular cardiomyopathy) were most significantly enriched in the intramuscular fat. There is no direct evidence showing that fat deposition in skeletal muscle caused the cardiac muscle diseases. However, excess body weight caused by chronic-over nutrition and high-fat feeding was associated with metabolic and/or diabetic cardiomyopathy in human [59].

ECM-receptor Interaction and Differentially Regulated Genes Based on Adipose Depot
Based on the functional annotation analysis of DEGs from different adipose tissue, we identified the ECM-receptor interaction pathway was involved in the various pathways of the three adipose tissues. The extracellular matrix (ECM) is essential for tissue architecture and has an important role in adipogenesis [60]. However it has not received sufficient attention in adipose tissue because of the difficulties associated with analysis of ECM components. ECM consists of a complex mixture of structural and functional macromolecules including glycosaminoglycans (GAGs) and fibrous proteins (collagen, elastin, fibronection, and lammin) [60]. The main constituents of ECM in adipose tissue are collagen (type I, IV, and VI), laminin (LN1, 8), fibronectin (FN), hyaluronan, and proteoglycan [61]. Specific interactions between cells and the ECM are mediated by transmembrane molecules such as integrins. In this study, different members of ECM fibrous proteins and integrins were up-regulated in the different adipose tissues. Collagens and integrins were significantly enriched in subcutaneous and intramuscular fat respectively. It seems that the interactions between ECM components and integrins were differentially regulated based on the adipose depots. Several studies reported that members of integrins showed the expression change associated with adipocyte differentiation [62]. Other genes involved in the ECM-receptor interaction were also showed that they were differentially regulated depending on the adipose depots. This comparative transcriptome analysis of 3 adipose tissues showed that the interaction between ECM components and transmembrane receptors of fat cells might influence the depot specific adipogenesis.

Supporting Information
Figure S1 Data quality control using FastQC.