This study collected 65 RNA sequencing data sets from 14 tissues (visceral adipose tissue, subcutaneous adipose tissue, intramuscular adipose tissue, muscle, heart, liver, spleen, lung, kidney, skin, brain, cerebellum, pituitary, and small intestine) of cattle (Bos Taurus). All these data were downloaded from NCBI GEO database. In order to identify important variations in adipose tissue-specific genes, this study collected and isolated DNA samples from 516 Shandong black cattle (with corresponding data of carcass traits, including the weight of cervical bone, bull collar, tendon, round small intestine, topside, hollow bone, rump steak, chuck tender, rib eye, tender loin, beef brisket, flank, chuck roll, and striploin), 59 Luxi cattle, 60 Bohai black cattle, and 37 Mengshan cattle in Shandong province, China.

The quality of raw sequencing data was assessed using FastQC v0.11.8 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) (15). In order to obtain clean reads, the sequencing data were filtered. The connectors at both ends of reads were removed, and the sequences containing low-quality bases were cut off. The clean reads were compared to the bovine reference genome (GCF_002263795.1_ars-UCD1.2_genomic.fna) using HISAT2 version 2.1.0 (16). The expression level (fragments per kilobase of transcript per million, FPKM) of all genes was calculated based on the bovine genome file (GCF_002263795.1_ars-UCD1.2_genomic.gFF) (15). A principal components analysis (PCA) was then performed, and some outliers were removed. Finally, 52 data sets were reserved for subsequent analysis, and each tissue containing 2–11 SRA data sets (Supplementary Table 1).

Since many genes were not expressed or had very low expression in the sequencing data collected in this study, genes with FPKM <1 in all tissues were removed, and a total of 20,446 genes were left after filtering. Then the tissue-specific index τ of each gene was calculated (17, 18). The τ value varies from 0 to 1. A τ value close to 1 indicates high tissue-specificity of gene expression, while τ value close to 0 indicates that the gene is widely expressed in all tissues (17, 18). For tissues with multiple biological replicates, the average value of FPKM values was used. Next, the top 20% genes were selected as candidate tissue-specific genes according to τ values. If the FPKM value of a candidate gene in a tissue is more than 1, and the expression level in that tissue is among the top three of all detected tissues, the candidate gene will be defined as a tissue-specific gene (19).

HTSeq Version 0.11.2 software was used to calculate the number of reads mapped to each gene according to the bam files generated after comparison (20). Subsequently, differential expression of genes was analyzed by DEseq2 package in R (P adjust < 0.05).

The KOBAS (http://kobas.cbi.pku.edu.cn/kobas3) was utilized for Gene Ontology (GO) clustering analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of candidate genes.

The bovine samples used for gene expression profile analysis were obtained from Kingbull Livestock Co., Ltd., (Yangling, Shaanxi, China). Tissues (spleen, lung, kidney, muscle, visceral fat, brain, testis, ovary) from four calves (two male and two female) were collected. RNA was isolated from tissue samples using Trizol reagent (Takara, Dalian, China). Consequently, RNA was used as template to prepare cDNA using Prime Script™ RT Reagent kit (Takara, Dalian, China). The prepared cDNA was preserved at −20°C.

The PCR and agarose gel electrophoresis were used to determine in which tissues the gene is expressed. Then, genes expressed in more than four tissues were selected for quantitative reverse transcription-polymerase chain reaction (qRT-PCR) assay to analyze their expression levels in different tissues. The primers were designed using NCBI Primer-BLAST (Table 1). The 25 μL reaction system contained 12.5 μL 2 × ChamQ SYBR qPCR Master Mix (Vazyme biotech co., ltd. Nanjing, China), 10 ng cDNA, 5 pM of each primer. The reaction procedure is as follows: pre-denaturation at 95°C for 45 s; 40 cycles of denaturation at 95°C for 15 s, annealing and elongation at 60°C for 40 s. The 2^−ΔΔCt method was applied to compute the relative gene expression levels in tissues (21). The glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was employed as an internal control.

The candidate important genetic variations of adipose tissue-specific genes were selected from the Ensembl database. Two pairs of primers were designed using NCBI Primer-BLAST to detect the authenticity and identify frequencies of the candidate variation (rs720343880) in different cattle breeds (Table 1). The touchdown-PCR method and 3.5% agarose gel electrophoresis were utilized to identify the genotypes of individuals. Sanger sequencing was used to confirm the accuracy of genotyping. Genetic parameters were computed on the MSRcall website (http://www.msrcall.com) according to the methods described as Nei (22), including observed heterozygosity (Ho), expected heterozygosity (Ne), polymorphism information content (PIC), and Hardy-Weinberg equilibrium (HWE) (22). The association analysis was performed based on the model similar to that in a previously published study (23). Significance between different genotypes and carcass traits was determined using an independent-sample t-test (the number of individuals in a group of less than three was not calculated).

Principal component analysis was carried out on the collected transcriptome data to eliminate outlier data sets. The remaining 52 data sets were analyzed, and the results showed that the rest of the data sets replicated well (Figure 1A). After removing genes with low expression in all tissues (FPKM < 1), a total of 4,090 genes with tissue-specific expression were detected. At the same time, cluster analysis based on the level of tissue-specific gene expression showed that different data sets from the same tissue aggregated well, and the three adipose tissues (subcutaneous adipose, visceral adipose, and intramuscular adipose tissues) clustered together (Figure 1B). The above results indicated that the data sets could represent the gene expression levels in different tissues, which can be used for subsequent analysis.

To date, a total of 35,158 genes have been annotated in the bovine genome. Analysis of bovine transcriptome annotation data showed that 58.87% of the genes are protein-coding genes, 14.73% are lncRNA genes, 2.27% are miRNA genes, and other RNA genes (including tRNA genes, pseudogenes, piRNA genes, etc.) account for 23.14% (Figure 2A). Among the 4,090 tissue-specific genes, 69.29, 19.61, 1.64, and 9.46% were annotated as protein-coding genes, lncRNA genes, miRNA genes, and other RNA genes, respectively (Figure 2B). These results indicated that lncRNA have higher tissue expression specificity (3).

Results of statistical analysis showed the number of genes specifically expressed in tissues, including brain (n = 987), pituitary gland (n = 883), small intestine (n = 878), cerebellum (n = 859), skin (n = 654), liver (n = 542), kidney (n = 518), spleen (n = 455), lung (n = 425), muscle (n = 337), heart (n = 279), intramuscular adipose tissue (n = 223), subcutaneous adipose tissue (n = 159), and visceral adipose tissue (n = 117) (Figure 2C; Supplementary Table 2). In all tissues, protein-coding genes were the main gene type (Figure 2D). Tissue-specific expressed genes may be related to tissue multifunctionality, tissue-specific function and cell type. These tissue-specific expressed genes are crucial for revealing the characteristics of different tissues and the functions of different genes (24, 25).

Due to the essential role of fat in cattle development and beef quality, we focused on genes specifically expressed in adipose tissue. Among the ATSGs, we found 23 genes specifically expressed in all three adipose tissues, including 16 protein-coding genes and seven lncRNA genes (Figure 3). The 16 protein-coding genes are as follows: ADAMTS16 (ADAM metallopeptidase with thrombospondin Type 1 Motif 16), ADIG (Adipogenin), COMP (cartilage oligomeric matrix protein), EN1 (Engrailed 1), GIPR (Gastric inhibitory polypeptide receptor, or glucose-dependent insulinotropic peptide), GSC (Goosecoid homeobox), HCAR1 (Hydroxycarboxylic acid receptor 1), KERA (Corneal protein), LEP (Leptin), LHCGR (luteinizing hormone/choriogonadotropin receptor), LOC101907335 (antigen WC1.1-like), LOC107131843 (insulin receptor substrate 1-like), LOC616957 (apolipoprotein L6, APOL6), MDFIC2 (MyoD Family Inhibitor Domain containing 2), S100A5 (S100 Calcium binding protein A5), and TRARG1 (also known as SLC2A4 (Solute carrier family 2 member 4) and GLUT4 (Glucose transporter type 4). The 7 lncRNA genes are LOC100847835, LOC104969981, LOC104970976, LOC112446369, LOC112448034, LOC112448366, and LOC112449245. There were 49 genes belonged to subcutaneous adipose tissue-specific genes (SATSGs) and visceral adipose tissue-specific genes (VATSGs), including 29 protein-coding genes, 18 lncRNA genes, one pseudogenes, and one miRNA gene (MIRN214) (Figure 3). Moreover, SATSGs and intramuscular adipose tissue-specific genes (IATSGs) included 12 genes: 10 protein-coding genes and two lncRNA genes (Figure 3). In addition, there were 14 genes belonged to IATSGs and VATSGs, including nine protein-coding genes, four lncRNA genes, and one pseudogene (Figure 3). Further statistical analysis showed that a total of 378 genes and 257 protein-coding genes were specifically expressed in adipose tissue (Figure 3).

Analysis of the expression levels of genes specifically expressed in visceral adipose tissue (VAT) in all tissues showed that the expression level of these genes was highest in VAT, followed by two other types of adipose tissue, subcutaneous adipose tissue (SAT) and intramuscular adipose tissue (IAT) (Figure 4A). Analysis of the expression levels of genes specifically expressed in SAT in all tissues exhibited that the expression level of these genes was greatest in SAT, followed by the other two types of adipose tissues (Figure 4B). The expression level of genes specifically expressed in IAT was the highest in IAT, but was lower in the other two adipose tissues, indicating that the gene expression patterns of VAT and SAT were more similar (Figure 4C). Besides, many IATSGs were highly expressed in muscle, suggesting that gene expression patterns between IAT and muscle may be more similar (Figure 4C).

Correlation analysis based on gene expression levels (log₂(FPKM+1) of 378 adipose tissue-specific genes revealed a significant correlation between the VAT and SAT (R = 0.819, P = 7.64E-93), a smaller correlation coefficient between VAT and IAT (R = 0.137, P = 0.008), and surprisingly, there was no significant correlation between SAT and IAT (R = 0.050, P = 0.328) (Figure 5). Next, as expected, significant positive correlations were found between IAT and muscle tissue (R = 0.360, P = 5.37E-13), while significant negative correlations were found between muscle and SAT (R = −0.261, P = 2.556E-07), and muscle and VAT (R = −0.197, P = 0.0001) (Figure 5). This might be due to the fact that the genes we used for the correlation analysis are ATSGs, and these genes are more likely to have a positive regulatory effect on adipogenesis, but a negative effect on other tissue development.

In order to investigate whether the above correlation is due to the genes we used were ATSGs, this study used all expressed genes to conduct correlation analysis between samples. The results displayed that the correlation coefficient R between the three types of adipose tissues and muscle was >0.50 (P < 0.001, Table 2). Further analysis demonstrated that using the expression levels of all genes has a problem: some frequently expressed genes and some genes with very low expression levels have a large influence on correlation and significance analysis. Therefore, tissue-specific genes were used for correlation analysis. As a result, the correlation coefficient R between VAT and SAT, VAT and IAT, SAT and IAT, IAT and muscle were 0.975 (P < 0.001), 0.244 (P = 2.81E-56), 0.233 (P = 1.68E-51), and 0.116 (P = 9.68E-14), respectively. In addition, there was no significant relationship between the muscle and the other two adipose tissues. Thus, these results indicated that although the gene expression patterns in IAT and muscle were similar, IAT was still more similar to the other two adipose tissues. Furthermore, regardless of which gene group were used, the gene expression patterns of SAT and VAT were the most identical, followed by IAT.

Among tissue-specific genes, protein-coding genes account for the largest proportion, and the current annotation information on protein-coding genes is the most comprehensive. Therefore, the function of the adipose tissue-specific protein-coding genes was analyzed. Based on Go analysis, three categories were identified: molecular function, cellular component and biological process. The first 20 terms with a corrected P < 0.05 in biological process and all corrected P < 0.05 terms in the other two categories were listed (Figure 6). The KEGG analysis results were given in Figure 7, and P-value was the corrected P-value. Some terms were annotated in all three adipose tissues, such as extracellular space and calcium ion binding. Additionally, some pathways were tissue-specific, such as metalloendopeptidase activity, G protein-coupled peptide receptor activity, and actin filament binding, which may be related to tissue-specific functions.

The above analysis identified 117, 159, and 223 genes specifically expressed in VAT, SAT, and IAT, respectively. However, the expression level of some of these genes was low. This may be because the top three genes of FPKM in all tested tissues were selected as tissue-specific genes, but some of these genes rank second or third in adipose tissue are poorly expressed. Therefore, in order to screen for genes that regulate adipose tissue function, the current study analyzed genes significantly higher expressed in adipose tissue (ATHEGs) from 52 cattle transcriptome data sets and combined them with ATSGs.

A total of 27,984 genes were detected in this study, and the differentially expressed genes in the various tissues were examined based on the number of reads of these genes. Based on the results of the correlation analysis of gene expression levels in tissues, the differential expression analysis in this part considers these three adipose tissues as an adipose tissue group and the other tissues as a non-adipose tissue group, and calculates the differentially expressed genes between the two groups. Consequently, a total of 1,575 genes were ATHEGs, while 5,840 genes had significantly lower expression in adipose tissues.

Furthermore, differential expression analysis of the three adipose tissues and other tissues showed that 798, 2, and 127 genes were significantly highly expressed in subcutaneous adipose (SATHEGs), visceral tissues (VATHEGs), and intramuscular adipose tissues (IATHEGs), respectively. The genes highly expressed in the three adipose were further analyzed, and no gene were found to be common among the three sets of genes. Only three genes were both SATHEGs and IATHEGs (Figure 8A). This result indicated that the gene expression patterns in the three adipose tissue types were identical, and the use of the three adipose tissue types as an adipose tissue group to analyze differentially expressed gene was reasonable.

A total of 23 genes (16 protein-coding genes and seven lncRNA genes) were found to belong to both ATSGs and ATHEGs (Figure 8B). However, no genes belonging to both VATSGs and VATHEGs were detected; 56 genes belonged to SATSGs and SATHEGs; and 19 genes were a memberr of IATSGs and IATHEGs (Figures 8C–E). These tissue-specific genes with high expression in adipose tissue may be related to their specific functions and locations.

In accordance with the transcriptome data of 23 genes, we made an expression profile of these genes, which displayed that these genes are highly and specifically expressed in bovine adipose tissue (Figure 9). Then, we used RT-PCR and qRT-PCR to determine the expression profiles of the 23 candidates in 8 bovine tissues, including spleen, lung, kidney, muscle, fat, brain, testis, and ovary. The results of RT-PCR showed that the mRNA expression of TRARG1 was detected only in fat in the analyzed tissues (Figure 10A). Some other genes, such as ADIG, HCAR1, APOL6, LOC107131843, LOC100847835, LOC104969981 were not only expressed in fat, but expressed more strongly in fat (Figure 10A). Besides, for genes expressed in more than four tissues, we detected its expression profiles using qRT-PCR. The results exhibited that GSCs had the highest expression in fat (Figure 10B). Meanwhile, several genes, involving LHCGR, MDFIC2, LOC104970976, LOC112448366 were highly expressed in fat, but also had high expression levels in some other tissues (Figure 10B). However, a few genes, consisting of S100A5, LOC112446369, and LOC112449245 are weakly expressed in fat (Figure 10B). Experimental verification results show that some candidate genes, such as TRARG1, LOC107131843, LOC100847835, and LOC104969981 screened from sequencing data can be preferentially selected for the study of bovine adipose tissue development.

The promoter is very important for gene expression. In this study, we tested a 12 bp insertion/deletion (indel) variation (rs720343880) in the promoter region of the bovine adipose tissue-specific lncRNA gene, LOC100847835. From the results of the sequence alignment based on the Ensembl database, it can be found that the lncRNA gene LOC100847835 and its promoter are conserved among species, especially among ruminants (Figure 11A). Based on sequences alignment of the candidate variation (rs720343880) region, it can be seen that the 12 bp insertion sequence is fully consistent with that of the domestic yak genome (Figure 11B). The sequencing result and agarose gel electrophoresis pattern revealed that the 12 bp indel (rs720343880) does exist in the tested population of interest (Shandong black cattle, Luxi cattle, Bohai black cattle, and Mengshang cattle) (Figures 11C,D). Furthermore, genotyping results showed that there were three genotypes (deletion/deletion, DD), insertion/deletion (ID), and insertion/insertion (II) in Shandong black cattle and Luxi cattle, but only two DD and ID genotypes in Bohai black cattle and Mengshang cattle (Table 3). The results of the genetic parameters calculation demonstrated that rs720343880 was in Hardy Weinberg equilibrium in the four breeds (Table 3). The association analysis results revealed that different genotypes of rs720343880 were associated with Shandong black cattle cervical bone, bull collar, tendon, hollow bone, round small intestine, topside, hollow bone, and rump steak (Table 4). Collectively, these results mean that this variation has an important impact on cattle production.