RNA sequencing of bovine subcutaneous adipose tissues

To obtain a global view of the bovine adipose tissue transcriptome, total RNA from the subcutaneous adipose tissues of fetal bovines, adult bulls, adult heifers and adult steers (ATFB, ATAB, ATAH and ATAS) were used to construct RNA libraries for Illumina sequencing. In total, we acquired an average of 27 billion clean reads from each library. Approximately 80% of the total reads uniquely mapped to the reference genome (http://hgdownload.soe.ucsc.edu/goldenPath/bosTau7/bigZips/bosTau7.fa.masked.gz). A total of 3% of sequences matched multiple positions in the reference genome, and 17% did not map to the reference genome (Table 1). Only the uniquely mapped reads were considered in this analysis. Gene expression level was calculated using the reads per kilobase transcriptome per million mapped reads (RPKM) method. To validate the RNA-Seq results, 15 genes that were uniquely or highly expressed in one of the four samples were selected and normalized to β-actin via the qRT-PCR method (Figure S1). We calculated the correlation efficiency of qRT-PCR and RNA-seq data using SPSS software (version 1.70), and the result (Figure 1) supported our sequence data (P<0.01, R = 0.83).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Verification of gene expression analysis by quantitative realtime PCR (qRT-PCR).

Note: Individual gene expression ratios were calculated using RPKM data generated by RNA-seq and plotted against calculations done for the same gene using qRT-PCR.

https://doi.org/10.1371/journal.pone.0101261.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of sequence read alignments to the reference genome.

https://doi.org/10.1371/journal.pone.0101261.t001

Analysis of Gene Expression in Adipose tissue

In our study, we detected a total of 12,233 genes expressed in the four samples. To better categorize these genes, which presented differential expression levels, gene expression RPKM values were categorized into three groups: high (≥500 RPKM), medium (10 to 500 RPKM), and low (<10 RPKM) (Table 2). There were 184 (1.60%), 181 (1.58%), 169 (1.45%) and 167 (1.44%) genes highly expressed in ATFB, ATAB, ATAH and ATAS, respectively. Genes that are highly expressed in a specific tissue may be responsible for the basic metabolism and function of that tissue. In the four expression patterns that were detected, 109 genes were expressed commonly, among which 63 genes were related to the protein components of ribosomes, and all of these genes were more highly expressed in ADFB than in the adipose tissue of adult cattle (Table S1). To ascertain the other potential functions of adipose tissue, we performed GO analysis for the other 46 genes. These results showed that those genes were involved in various cell component, molecular function and biological processes, and no GO term was significantly enriched (P<0.05) (Figure S2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. RNA-Seq gene expression results for the four patterns of adipose tissues.

https://doi.org/10.1371/journal.pone.0101261.t002

Age plays an important role in the development of adipose tissue. We compared the expressional differences in adipose tissue between fetal bovine at 200∼250 days, at which time the white adipose tissue just emerges and its development starts, and adult cattle at 20 months, whose adipose tissue increases mainly due to the accumulation of triglycerides. In this study, the expression profiles of adipose tissue changed considerably from fetuses to adults. There were 3,188, 4,167 and 3,903 genes differently expressed in ATAB, ATAS and ATAH compared to ATFB, respectively. We selected a total of 2,703 genes that were differentially expressed between the two stages using the edgeR method [19]. During development, sex is not only an important turning point in metabolism, but is also a key factor resulting in individual differences. To detect the effect of sex on adipose tissue development, we compared the expression profiles of adipose tissues in adult bulls, adult cows and adult steers at 20 months of age. A total of 847 genes were found to differ significantly in expressional levels between ATAB and ATAH. The number of DEGs between ATAS and ATAB (1,464 DEGs) was approximately three times higher than between ATAS and ATAH (495 DEGs). When mixed together, there were 57 genes differentially expressed among the three adipose tissues of adult cattle and 38 genes differentially expressed among all the four adipose tissues. The distribution of differential expression genes among the four adipose tissues was shown in Figure 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Distribution of differential expression genes among the four adipose tissues.

Note: ATFB: adipose tissues of fetal bovines; ATAB: adipose tissues of adult bulls; ATAS: adipose tissues of adult steers; ATAH: adipose tissues of adult heifers; AT: adipose tissues

https://doi.org/10.1371/journal.pone.0101261.g002

Gene Ontology Analysis of the DEGs

To further investigate the biological processes associated with the 2,703 differentially expressed genes, we performed Gene Ontology (GO) analysis by running queries for each differentially expressed gene against the GO database, which provides information related to three ontologies: molecular function, cellular component, and biological process. First, we arranged GO analysis for DEGs between ATFB and adult bovine adipose tissue, and the results are shown in Table S2. In total, we found that 2979 GO terms were assigned to the DEGs, and among these genes, there were 251 GO terms corresponding to molecular function, 519 GO terms corresponding to cellular component and 2209 GO terms corresponding to biological process. Structural molecule activity and receptor activity were the only two significantly enriched molecular functions. In the cellular component category, 9 terms were significantly enriched (P<0.05). The four most significantly enriched GO terms included the structural constituents of ribosomes, followed by genes expressed in the extracellular region, ribonucleoprotein complex, extracellular matrix, and genes that are integral to membrane structures. In the case of biological processes, 22 GO terms were significantly enriched and were related to various processes, such as organ development, cell proliferation, and the immune system. Lipid metabolic processes were not found to be enriched significantly, with a comparatively low corrected P-value (0.053).

GO analysis was performed for DEGs between the adipose tissues of male and female cattle (Table S3–S5). When comparing ATAB and ATAS, we detected 2311 GO terms, including 193 GO terms in the molecular function category, 400 GO terms in the cellular component category and 1718 GO terms in the biological process category. In the molecular function category, GO terms for carbohydrate derivative binding and receptor activity were significantly enriched. There were 9 GO terms in the molecular function category that were significantly enriched, which were mostly related to the cell membrane, especially to the mitochondrial membrane. In the biological process category, we found 33 GO terms that were significantly enriched, including the GO terms that we expected, such as steroid hormone stimulus and response to lipids. Three GO terms related to cytokines were also significantly enriched.

When comparing ATAB and ATAH, we detected 1762 GO terms, including 148 GO terms in the molecular function category, 283 GO terms in the cellular component category and 1330 GO terms in the biological process category. In the molecular function category, there were 7 significantly enriched GO terms. Apart from the GO terms that were significantly enriched between ATAB and ATAS, cytokine activity, signal transducer activity, molecular transducer activity, motor activity and transmembrane signaling receptor activity were also significantly enriched. In the cellular component category, the DEGs that were most enriched included those categorized to the extracellular region. In the biological process category, we found 18 GO terms that were significantly enriched. Similar to the comparison between ATAB and ATAS, GO terms related to steroid hormone stimulus, response to lipids and cytokines were significantly enriched.

Comparing ATAH and ATAS, we detected 1310 GO terms, including 125 GO terms in the molecular function category, 215 GO terms in the cellular component category and 970 GO terms in the biological process category. In the molecular function category, the most significantly enriched were carbohydrate derivative binding, followed by oxidoreductase activity, transferase activity and cell surface binding. The DEGs were significantly enriched in cellular components of the extracellular matrix, intrinsic to membrane, extracellular region part and in the extracellular region. In the biological process category, 21 GO terms were significantly enriched. Compared to the other two comparisons between adipose tissues of adult bovine, GO terms related to steroid hormones were also included in the 21 GO terms, and more processes related to lipid metabolism were significantly enriched, including the following: fatty-acyl-CoA metabolic process, lipid metabolic process, neutral lipid metabolic process, acylglycerol metabolic process, acyl-CoA metabolic process, triglyceride metabolic process and regulation of the lipid biosynthetic process. The estrogen stimulus response process was also found to be significantly enriched, while no GO term related to cytokines was found to be significantly enriched.

KEGG Pathway Analysis of the DEGs

Different genes usually cooperate with each other to exercise their biological functions. Pathway-based analysis helps to further understand the biological functions of genes. Overall, 228 pathways (17 pathways significantly enriched, P<0.05) were assigned to the DEGs in the comparison of ATFB and adipose tissue from adult bovines (Table S6). The most significantly enriched pathway was ribosomal. The PPAR signaling pathway, one the most important pathways related to the development of adipose tissue, was also found to be significantly enriched. Comparing ATAB and ATAS, 201 pathways (16 pathways significantly enriched, P<0.05) were assigned to the DEGs (Table S7). The pathway term showing the highest level of significance was type I diabetes mellitus, followed by oxidative phosphorylation, pathways in cancer, ECM-receptor interaction, proximal tubule bicarbonate reclamation, p53 signaling pathway, sphingolipid metabolism and cytokine-cytokine receptor interaction. Comparing ATAB and ATAH, we detected 183 pathways (13 pathways significantly enriched, P<0.05) assigned to the DEGs (Table S8). Type I diabetes mellitus and cytokine-cytokine receptor interaction were also among the significantly enriched pathways. In the comparison between ATAH and ATAS, there were 170 pathways (12 pathways significantly enriched, P<0.05) assigned to the DEGs (Table S9). Steroid biosynthesis, instead of type I diabetes mellitus, was the most significantly enriched pathway. Furthermore, biosynthesis of unsaturated fatty acids, fatty acid elongation and the PPAR signaling pathway also showed significance in our analysis.

Differential alternative splicing in adipose tissue

The study of alternative splicing (AS) has long been a valuable subfield of molecular biology, but this subfield was not given enough attention until recent years. AS plays a major role in the generation of proteomic and functional complexity in higher organisms [20]. In our data, we detected 12,233 genes, of which up to 4,753 genes underwent AS events. Therefore, in total, approximately 38.85% of the genes in adipose tissues were found to be alternatively spliced in our study. However, single-sample analysis showed that more than 50% of the AS events were specifically expressed in the adipose tissues we studied, thereby indicating a strong association of these AS events with the regulation of adipose tissues in specific physiological environments. Taking this phenomenon into consideration, it must be noted that the number of genes that underwent AS events in adipose tissues may have been far higher than 4,753 because a single tissue can encounter various physiological environments during development. These data suggested a major role for AS in contributing to the complexity of gene expression in adipose tissue.

To identify the types of AS mechanisms being used by genes expressed in adipose tissues, we conducted an analysis of AS events in the four types of adipose tissue (Figure 3). All seven known types of AS models were found in ATAB. These seven types included exon skipping (ES), intron retention (IR), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), alternative first exons (AFE), alternative last exons (ALE) and mutually exclusive exons (MXE). In the other three adipose tissues, there were six known types of AS models. Here, no AS events were categorized as the mutually exclusive exon model, and this type of AS was only found twice in ATAB. Almost all AS events in the four adipose tissues belonged to the four models: ES, IR, A5SS, A3SS, and only a few AS events occurred in the other three models. However, analysis of the distributions of the main four AS models highlighted a few differences between ADFB and the adipose tissue of adult cattle: ES and A3SS had a comparative advantage in ADFB, while IR and A3SS were more abundant than ER and A5SS in ADAC.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Statistic charts of alternative splicing events and genes in the four adipose tissues.

Note: ATFB: adipose tissues of fetal bovines; ATAB: adipose tissues of adult bulls; ATAS: adipose tissues of adult steers; ATAH: adipose tissues of adult heifers.

https://doi.org/10.1371/journal.pone.0101261.g003

Identification of novel transcribed units

The current gene annotation is almost exclusively based on gene prediction programs that detect protein coding regions and do not include gene coding for small proteins (<100 codons) [21]. Non-coding RNAs (ncRNAs) appear to comprise a hidden layer of internal signals that control various levels of gene expression in physiology and development [22]. However, very few ncRNAs, with the exception of microRNAs, have been found and analyzed in adipose tissue to date. To identify novel transcription units, we compared the distribution of assembled reads with the annotation of reference genes using stringent criteria. In total, 4,020, 4,539, 4,713, and 4,183 novel transcription units (NTUs) were found in ATAB, ATAS, ATAH, and ATAF, respectively. The Coding Potential Calculator was used to predict the coding ability of the NTUs. The majority of the NTUs were not expected to encode proteins, and only approximately 30% of the NTUs were considered to have coding ability. Of the NTUs with coding ability, 45, 43, 46, and 31 NTUs had fewer than 100 codons in ATAB, ATAS, ATAH, and ATAF, respectively. Next, we analyzed the distribution of all of the NTUs. First, the increase in the number of NTUs per 1,500 bp was calculated and suggested two basic clusters of tendency according to the coding ability of the NTUs. In table 3, the numbers of NTUs under 1,500 bp both with and without coding ability occupied over 50% in the four samples and decreased as an index per 1,500 bp. This illustrated that both NTUs with and without coding ability were enriched under1,500 bp and that the main difference in the numbers of the two types of NTUs was also concentrated in this range.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Distributions of the novel transcript units with and without coding ability.

https://doi.org/10.1371/journal.pone.0101261.t003

Detecting SNP Variants in Pooled Transcriptome Samples by RNA-Seq

RNA-seq has been shown to be an efficient way to detect SNP variants on a large scale at the mRNA level. In bovine milk, Angela Cánovas and his colleagues discovered 100,734 SNPs in Holstein samples using RNA-seq and validated these results by Sanger sequencing technologies, which confirmed that RNA-Seq technology is an efficient and cost-effective method to identify SNPs in transcribed regions [23]. In our study, to detect the SNPs at the transcriptome level, we applied SOAPsnp software and detected 56,564, 65,154, 78,061 and 86,965 putative SNPs in ATFB, ATAB, ATAH and ATAS, respectively. These results are shown in Table_S10 to Table_S13.

Ethics statement

All animals in this study were approved by Institutional Animal Care and Use Committee (IACUC) of Northwest A&F University and Qinbao Animal Husbandry Co., Ltd, respectively. All surgery was performed under sodium pentobarbital anesthesia, and all efforts were made to minimize suffering. Bovine embryos of slaughtered cows were collected from the Tumen abattoir, a local slaughterhouse of Xi'an, P.R. China. Heifers, bulls, and steers were obtained from Shanxi Kingbull Livestock Co., Ltd.

Sample collection and RNA sequencing

Bovine embryos at approximately 250 days (gestation period 280 days) were collected from the reproductive tracts of slaughtered cows. The adult animals were kept in free-stall housing, fed total mixed ration, and slaughtered at 20 months. In this study, five bovine embryos, five adult cows, five adult bulls, and five adult steers of Qinchuan cattle were selected. The subcutaneous adipose tissues were collected immediately after the cattle were slaughtered, snap-frozen in liquid nitrogen, and stored at −80°C until use.

Using the above samples, we constructed four pooled RNA libraries. Total RNA was extracted from adipose tissues using the Trizol method according to the manufacturer's instructions (TaKaRa, Dalian, China). The RNA was treated with DNase I for 30 min at 37°C to remove residual DNA. Then, poly(A) mRNA was isolated using beads with oligo(dT). The purified mRNA was first fragmented to approximately 200 nt fragments using the RNA fragmentation kit. First-strand cDNA synthesis was performed using random hexamer primers and reverse transcriptase. The first strand cDNA was synthesized using the mRNA fragments as templates. After the first strand was synthesized, a custom second strand primer and strand synthesis buffer (Illumina) were both added, followed by dNTPs, RNase H and DNA polymerase I to nick translate the second-strand. The DNA product was then purified using a QiaQuick PCR column and eluted in EB buffer; the samples were then subjected to end repair, addition of a single A base, adaptor ligation, agarose gel isolation of approximately 200 nt cDNA, and PCR amplification. Then, the cDNA libraries were prepared according to Illumina's protocols and sequenced on the Illumina platform in BGI, Shenzhen, China.

Read Mapping on the Bovine Reference Genome and Data Analysis

After the 3′ adaptor sequence was cleared away and redundant and low-quality reads were removed, the clean reads were aligned to the reference genome (http://hgdownload.soe.ucsc.edu/goldenPath/bosTau7/bigZips/bosTau7.fa.masked.gz) for assembly. Clean reads were 100% matched to the reference genome using the SOAPaligner program by not allowing any mismatches in reads. For the reads that could not be aligned with the reference sequences, we performed the SOAP alignment again while allowing less than 3 bp mismatches in the reads. A similar strategy was used to filter the unique sequence reads and multi-position reads. The gene expression level was normalized by calculating the number of reads per kilobase transcriptome per million mapped reads (RPKM) method. We calculated the fold change and P-value [decided by controlling the false discovery rate (FDR)] of each gene in two samples. Finally, the following criteria were used to judge the significance of gene expression difference: (1) Log2Ratio≥1 and (2) FDR≤0.001 [51]. All the basic data series were submitted to National Center for Biotechnology Information Gene Expression Omnibus with accession number GSE47653.

Gene ontology and pathway enrichment analysis of DEGs

Gene ontology (GO) is an international standard gene functional classification system [52]. We first mapped the differentially expressed genes (DEGs) to GO terms in the database (http://www.geneontology.org/) to calculate gene numbers for every term and then used the corrected P-value≤0.05 as a threshold to find significantly enriched GO terms in the input list of DEGs. The P-value formula was as follows:where N is the number of all genes with GO annotation; n is the number of DEGs in N; M is the number of all genes annotated to certain GO terms; and m is the number of DEGs in M. The calculated P-value was subjected to Bonferroni correction [50].

Next, KEGG (http://www.genome.jp/kegg/) was used to perform a pathway enrichment analysis of DEGs [53]. The calculation formula was the same as that in the GO analysis. Here, N is the number of all genes with KEGG annotation, n is the number of DEGs in N, M is the number of genes annotated to a specific pathway, and m is the number of DEGs in M. Biological pathways with an FDR≤0.05 were considered significant.

Identification of novel transcribed units

To discover novel transcribed regions, we compared our assembled transcripts and annotated genomic transcripts from a reference annotated sequence. The assembled transcripts that did not align with the reference transcripts were then aligned with the reference genome. To be defined as a novel transcript, an assembled transcript must meet three requirements: (a) the transcript must be at least 200 bp away from an annotated gene, (b) the transcript has a length of over 180 bp, and (c) the sequencing depth is no less than 2. After discovering novel transcripts, we used the Coding Potential Calculator (CPC, http://cpc.cbi.pku.edu.cn/) to assess the protein-coding potential of the novel units.

Analysis of alternative splicing

We used the SOAPsplice software to detect alternative splicing (AS) events. First, all of the reads were mapped to the reference genome for an intact alignment. Then, the unmapped reads were prepared for the next spliced alignment using the spliced alignment algorithm and the splice junctions were reported. Finally, we used both the mapping results and the known splice junctions to detect seven AS events according to Wang et al [42].

Detection of Putative SNPs

We used SOAPsnp software to detect SNPs in the four samples based on the massively parallel Illumina GA technology. SOAPsnp is a member of the SOAP (Short Oligonucleotide Analysis Package) software [54]. The program calculates the likelihood of each genotype at each site based on the alignment of short reads to a reference sequence together with the corresponding sequencing quality scores. It then infers the genotype with highest posterior probability at each site based on Bayes' theorem (the reverse probability model). Thus, we have taken into account the intrinsic bias or errors that are common in Illumina GA sequencing data and recalibrated the quality values for use in inferring consensus sequence. The SNP detection stringency conditions were as follows: (1) at least two unique mapping reads that support the polymorphic nucleotide and (2) quality score ≥20.

Validation of differentially expressed genes

To validate the sequence results, 15 genes that were specially or highly expressed in one of the four samples were selected and normalized to β-actin via the quantitative real-time PCR (qRT-PCR) method. Total RNA (1.0 µg from tested tissues) was reverse-transcribed to cDNA using the PrimeScript RT reagent Kit with gDNA Eraser according to the manufacturer's instructions. Primers were designed and are shown in Table_S13. qRT-PCR was performed using a Bio-Rad CFX 96 Real Time Detection System and SYBR Green PCR Master Mix in a 20 µl reaction. The relative gene expression values were calculated using the 2^−ΔΔCt method.