Bioinformatics Analyses of Bovine Adipose Tissue Transcriptome from Lilu Beef Cattle at Different Stages of Growth

Fat is an unpopular part of meat for consumers and is considered to be unhealthy to their bodies in many countries. However, the quantity and distribution of fat in the adipose tissue are important contributors to meat quality, such as taste and nutritional quality of the meat (Lozeman et al., 2001; Wood et al., 2008; Bong et al., 2010). The beef carcass can be categorized into quality grades based on subcutaneous fat thickness, so breeding for optical carcass fat is one of the major goals toward better profitability in beef industry (Taniguchi et al., 2008). The growth and development of adipose is a valuable trait that affects marbling and an important factor determining the price of beef in the Chinese beef market.

Fat accumulation is affected by age, genetics, and nutrition. Metabolic processes affecting fat accumulation occur in adipose tissue, there are some certain relationships exist between adipose tissue metabolism and fat accumulation (Rule et al., 1992). Gene expression in fat depots provides further proof that the mechanisms for fat accumulation differ significantly among animal species (Hishikawa et al., 2005). In fact, it is a complicated biochemical process that involves many genes, developmental stages, and pathways. Some studies have reported regional differences in the expression of individual genes among different adipose depots (Kim et al., 2000; Wang et al., 2005). Adipocyte biology plays pivotal roles in the regulation of fat metabolism (Marzolla et al., 2012). Uncovering novel genes may provide deeper appreciation of the genetic basis for the growth and development of adipose.

Lilu cattle are the result of crossing fixed combined with directional excellent-choosing that includes the pedigree of 62.5% Limousin and 37.5% Chinese Luxi cattle. Lilu cattle have bigger and well-balanced body type. Adult bull and cow weight is above 800kg and 500kg, respectively. Hybrid off-springs have good production performance and adaptability. Currently, this breed will be under-review as a new beef cattle variety in China.

Microarray technology has been widely used to discover new genes and functions, to understand biochemical pathways (Bourzac et al., 2011; Wang et al., 2018). We used a novel bioinformatics approach, the dynamic impact approach (DIA) (Bionaz et al., 2012) to analyze microarray data at four different age stage in Lilu cattle. In addition, we used ingenuity pathway analysis (IPA) to study downstream regulators of transcriptomics differences. The primary purpose of this study was to uncover the gene expression and signaling or metabolic pathways of subcutaneous adipose tissues during development of Lilu beef cattle. As a result, the expression pattern of the adipose transcriptome of Lilu beef cattle could be discerned.

A total of 12 Li-Lu beef cattle, four groups (12, 18, 24, and 30 months of age, n=3) were used in this experiment. The whole procedure for experimental animals was performed in strict accordance with guideline (IACC20060101, Jan, 2006) of the Institutional Animal Care and Use Committee of Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences, and all efforts were made to minimize suffering. The subcutaneous adipose tissues near the last thirteenth or fourteenth rib were collected from animal carcasses after slaughter, immediately frozen in liquid nitrogen, and kept at -80°C until analyzed.

Total RNA extraction was performed by homogenizing the fat tissue samples with Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified using the RNAeasy mini kit (QIAGEN; Cat# 74106). RNA quality and purity were determined by NanoDrop ND-1000 spectrophotometer at 260/280 nm (NanoDrop Technologies, USA). The integrity of total RNA was assessed by an Agilent 2100 Bioanalyzer profile (Agilent Technologies Inc., USA); only samples with RNA Integrity Numbers (RIN) values more than 7 were used for microarray and qRT-PCR analysis.

Purified total RNA was transcribed into cDNA using the Low RNA Input Linear Amplication kit (Agilent Technologies Inc., USA). Then cDNA was transcribed into cRNA which was labeled with Cyanine-3 (Cy3) NHS ester (GE healthcare, USA) using T7 RNA polymerase (Agilent Technologies Inc., USA).

The Cy3 labeled RNA samples were hybridized with Agilent GF Bovine 4×4 gene expression microarray slides (Covering the entire bovine genome, including 43803 probes) at 65°C for 17 h. The hybridized microarray slides were washed with a Gene Expression Washing Buffer Kit (Agilent Technologies Inc.; Cat # 5188-5327) and were scanned with an Agilent DNA Microarray Scanner (#G2565BA, Agilent Technologies) at a resolution of 5μm.

Microarrays were adjusted for dye and array effect (Quantile normalization and array centering), duplicated spot intensities were not averaged and were subsequently used for statistic alanalysis. A mixed model with repeated measures was then fitted to the normalized log2-transformed adjusted ratios (sample/reference standard) using MIXED Procedure of SAS (version 9.1.3; SAS Institute, Cary, NC, USA). The model included the fixed effects of time (12, 18, 24, 30). Beef cattle were considered as a random effect. Differences in relative gene expression were considered significant at P<0.05 and fold change greater than 6. Statistical significance was accepted at the P<0.05 level.

Candidate genes were selected based on microarray data for qRT-PCR validation. qRT-PCR primers were designed and synthesized (Biosune, Shanghai) to assay 8 differentially expressed genes (Table I). qRT-PCR were performed using the ABI Prism 7700 Sequence Detector

Gene abbreviations: FOX1, forkhead box O1; FASN, fatty acid synthase; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; LEAP2, liver-expressed antimicrobial peptide-2; LRL, lipoproteinlipase; MHC, major histocompatibility complex; MYH, MutY human homologue; SIRT1, sirtuin 1.

(Applied Biosystems, Foster City, CA, USA) in a final volume of 20 μL containing SYBR Green I (Invitrogen). The cycling parameters of qRT-PCR amplification were as follows: 94°C for 2 min, 35 cycles of 94°C for 30 s, appropriate annealing temperature for 30 s, 72°C for 30 s.

The qRT-PCR data were genometric mean of the internal control genes and log2 transformed before statistical analysis. All means were compared using the PDIFF statement of SAS, and difference were considered statistically significant at P<0.05.

Quantile normalization: The goal of the quantile method is to make the distribution of probe intensities for each array in a set of arrays the same. The method is motivated by the idea that a quantile-quantile plot shows that the distribution of two data vectors is the same if the plot is a straight diagonal line and not the same if it is other than a diagonal line.

Microarrays were adjusted for dye and array effect (Quantile normalization and array centering), duplicated spot intensities were not averaged and were subsequently used for statistical analysis. A mixed model with repeated measures was then fitted to the normalized log2-transformed adjusted ratios (sample/reference standard) using MIXED Procedure of SAS (version 9.1.3; SAS Institute, Cary, NC, USA). The model included the fixed effects of time (12, 18, 24, 30). Beef cattle were considered as a random effect. Differences in relative gene expression were considered significant at P≤0.05 and no fold change for 6 comparisons. Statistical significance was accepted when P<0.05.

The qPCR data were genometric mean of the internal control genes and log2 transformed before statistical analysis. All means were compared using the PDIFF statement of SAS, and difference were considered statistically significant at P≤0.05.

The dynamic impact approach (DIA) was a method to study the biological impact of experimental conditions through the transcriptome level and it can interpret the biology of the impact by providing the direction of the impact (Bionaz et al., 2012). In this study, the DIA were used to uncover impact and direction (flux) of the biology pathway in four different age stages of Lilu cattle. Briefly, the data of gene chip was categorized as follows: Entrez gene ID, FDR, fold-change (FC), and P-value, and then the data was uploaded to DIA and calculated. The results are showed in figures and tables.

Ingenuity pathway analysis (IPA) software was used to analyze the upstream transcription regulators and their connections with other downstream genes that were differentially expressed. A list of DEG (P-values≤0.05) was uploaded to the IPA using core analysis. The results of upstream transcription regulators were downloaded and saved for further analysis.

All differently expressed genes of different ages were listed in Table II, the number of differently expressed genes were 188, 707, 1016, 907, 1026 and 165 between 12 vs 18, 12 vs 24, 12 vs 30, 18 vs 24, 18 vs 30 and 24 vs 30, respectively. The differential expression genes mainly focus on 12 vs 24, 12 vs 30, 18 vs 24, and 18 vs 30 and the numbers of genes are more than the other two groups. The down-regulated genes of these groups are more than up-regulated genes. There are no big differences between 12 vs 18 and 24 vs 30.

The differential expressed genes were selected according to P-value ≤ 0.05. No fold change in this study.

To validate the microarray hybridization results, 8 genes were selected from the differentially expression gene list for quantitative qRT-PCR assays. Because they are either important components of the lipid metabolic process and are related to meat quality of cattle (FOX1, FASN, HMGCR, LPL, and SIRT1) or because they are important and believable genes in the process of chip analyzing (LEAP2, MHC, MYH). The adipose tissues used in qRT-PCR assays were total RNA of all 4 age stages from 12 animals. For all 8 genes, the expression trends are similar both in microarray analyses and qRT-PCR experiment. However, the expression abundance was higher in the microarray analyses than in qRT-PCR of FASN, LEAP2, LPL genes. Altogether, validation of microarray data by qRT-PCR revealed a high correspondence between both analyses (Table III).

A functional analysis of DEG (P<0.05) was made with DIA. The DIA provide a summary of the KEGG pathways in the form of categories and sub-categories (Fig. 1), in accordance with the number of DEG (Table II) all the differentially expressed genes were more focus on 12 vs 30, 18 vs 24 and 18 vs 30. The summary of the KEGG pathways showed that KEGG pathway categories were more impacted 18 vs 24 and 18 vs 30 in parallel with the other time comparisons. The categories “metabolism” was the most impacted in the adipose tissue between 18 vs 24 and 18 vs 30, and the 18 month is very important age for beef cattle. The categories “metabolism” was the most impacted in the adipose tissue between 18 vs 24 and 18 vs 30, and the 18 month is very important age for beef cattle. The subcategories of pathway “energy metabolism” were activated from 12 month to 24 month, and then this pathway was inhibited on 30 month. From Figure 1, the most metabolic pathways were inhibited after 18 month of beef cattle.

*, means P<0.05; **, means P<0.01; #, means P>0.05; +, means up-regulation; - means down-regulation.

In metabolic pathway, sulfur metabolism; thiamine metabolism; glycosaminoglycan degradation; alpha-linolenic acid metabolism; linoleic acid metabolism; glycine, serine and threonine metabolism are the top 6 most affected metabolism pathways in the adipose tissue between 18 vs 24 (Fig. 2). The thiamine metabolism; glycine, serine and threonine metabolism; sulfur metabolism; cysteine and methionine metabolism; histidine metabolism and ascorbate and aldarate metabolism are top 6 most affected metabolism pathways in the adipose tissue between 18 vs 30 (Fig. 3).

In nonmetabolic pathway, fat digestion and absorption; vasopressin-regulated water reabsorption; renin-angiotensin system; leukocyte transendothelial migration; PPAR signaling pathway; dorso-ventral axis formation are the top 6 most affected non metabolism pathways in the adipose tissue between 18 vs 24 (Fig. 4). The thiamine metabolism; glycine, serine and threonine metabolism; sulfur metabolism; non-homologous end-joining; cysteine and methionine metabolism are the top 6 most affected non metabolism pathways in the adipose tissue between 18 vs 30 (Fig. 5).

We used IPA to find the relationships between transcription factors and all the DEGs between 18 vs 24 and 18 vs 30. Eighty-eight transcription regulators were found in the adipose tissue between 18 vs 24, KLF5, one of activated transcription regulators, which is relevant to cancer, cellular development, cellular growth and proliferation, and this regulator is unregulated FAM110A, MYC, DUSP1, ACTA2, NOTCH1 and PREB, down regulated RUNX1 (Fig. 6). Another eighty-eight transcription regulators were tested in the adipose tissue by IPA between 18 vs 30, PPARGC1A, one of inhibited transcription, which is related to lipid metabolism, down regulated 20 target molecules and up regulated 6 target molecules (Fig. 7).

LuXi yellow cattle is a good local beef cattle in Shandong province, which is famous for its excellent meat quality, however, its growth rate is slow. Therefore, we have used introduced breed to improve the local breed. Lilu beef cattle is breeded through the method of crossing fixed and combined with directional excellent-choosing individual which include 62.5% adventitious Limousin blood and 37.5% local Luxi Yellow cattle blood. We have analyzed the growth performance, carcass characteristics, and meat quality in four age stages. Most of growth traits are superior to Luxi yellow cattle in the same stage and there are not many dominant variations between them in meat quality. These results have been reported previously (Song et al., 2011, 2012).

In beef cattle, postnatal adipose tissue development was thought to occur as follows: internal, intermuscular, subcutaneous, and, last intramuscular. The accumulation of fat is associated with the genetic background, development, and nutrition of an animal. The adipose are very important because they can store fats as energy sources (Lee et al., 2010), and fat related carcass traits are important to the beef industry with value of the meat (Taniguchi et al., 2008). However, the growth and differentiation of adipose tissue were controlled by many factors, including heredity, nutrition and others. A number of quantitative trait loci (QTL) studies have shown candidate genes and their location on chromosomes that may be associated with backfat thickness in cattle (Moore et al., 2003; Wu et al., 2005). Some results indicated that adipocyte differentiation were controlled by a number of adipogenic transcription factors (Rangwala and Lazar, 2000; Rosen and Spiegelman, 2000), and adipocyte differentiation were also regulated by differentiation-inducing agents (Girard et al., 1994; Reusch et al., 2000).

From overall differently expression gene pattern of this study, we found that the differently expressed genes were focusing on 18 vs 24, 18 vs 30 and 24 vs 30. However, there are fewer genes existed between 12 vs 18, 24 vs 30. As the reported, during the development of beef cattle, the growth rates of different parts and tissues of the body are not the same. Head, limbs and skeleton are the fastest at the early age, body length and muscle mainly grow after the early age, while after maturity; the growth of body weight and fat is the main thing (Wang et al., 2005). Form these differently expression genes numbers, we can concluded that 18 month age is very important development stage in beef cattle.

Fat deposition and metabolism, is affected by several lipid metabolism pathway. In metabolic pathways, the most of pathways are inhibited, especially between 18 vs 24 and 18 vs 30. The lipid metabolism, metabolism of other amino acids and metabolism of cofactors and vitamins were most inhibited between 18 vs 24, also in 18 vs 30. In beef cattle, people always to improve lipid level by castration method. In fact, the steers have higher lipid than bulls, because liver fat deposition is affected by several lipid metabolism pathways, steers had lower hepatic mRNA levels of LPL and FATP5 than bulls. Steers also had higher adipogenic SREBPs mRNA levels (Baik et al., 2015). The differently expressed genes were analyzed by Gene chip. Gene chip has been widely used in bioscience research as a technical tool. Gene chip has been used in many species, cattle (Okumura et al., 2007; Jiang et al., 2009), pig (Sun et al., 2013), sheep, and others (Yu et al., 2010). Agilent GeneChip oligonucleotide arrays are a popular platform for the high-throughput analysis of gene expression or discovering new genes (Lee et al., 2010). The differently expressed genes have three main trends in four age stages of Lilu beef cattle and these genes are critical in PPAR signaling pathway. In our group, we have researched that three-hydroxy-3-methylglutaryl-CoA reductase (HMGCR) catalyzes the conversion of HMG-CoA to mevalonate and is considered the rate-limiting enzyme in the overall pathway of cholesterol biosynthesis (Istvan and Deisenhofer, 2000) and have an important role in PPAR signaling pathway. These genes involved signaling pathway which would contribute to our further research in Lilu beef cattle.

The induction of adipocyte differentiation and the resulting appearance of new small adipocytes after induction of PPAR have been reported by Yamauchi and Kadowaki (2001) in murine models of insulin resistance. PPAR, a key regulator in adipogenesis and fat storage, controlling the expression of many adipocyte-specific genes (Picard et al., 2004). However, the development of adipocyte size in cattle from 10 to 26 months of age does not correspond with a continuous increase in s.c. fat, which is consistent with steady-state levels of PPAR gene expression (Albrecht et al., 2011). Up regulation of PPARγwas observed in the backfat tissue of Lilu cattle with the increasing age (Liu et al., 2014). In present research, PPAR signal pathway is be active between 12 vs 30, and the most impact exists between 18 vs 30. As for beef cattle, after maturity, the fat will deposit very quickly in different body part.

We acknowledge the National Natural Science Foundation of China under Grant No. 31402098), Young Talents Training Program of Shandong Academy of Agricultural Science, MATS-Beef Cattle System, the Sustentative Research Project of China Ministry of Science and Technology under Grant No. 2015BAD03B04), Breeding New Varieties Projects of Transgenic Organisms under Grant No. 2016ZX08007-002), Agricultural Science and Technology Innovation Project of Shandong Academy of Agricultural Sciences under Grant No. CXGC2016A04).

Authors have declared that there is no conflict of interests regarding the publication of this article.