Comparative Transcriptome Analysis of Key Genes and Pathways Activated in Response to Fat Deposition in Two Sheep Breeds With Distinct Tail Phenotype

Fat tail in sheep presents a valuable energy reserve that has historically facilitated adaptation to harsh environments. However, in modern intensive and semi-intensive sheep industry systems, breeds with leaner tails are more desirable. In the present study, RNA sequencing (RNA-Seq) was applied to determine the transcriptome profiles of tail fat tissues in two Chinese sheep breeds, fat-rumped Altay sheep and thin-tailed Xinjiang fine wool (XFW) sheep, with extreme fat tail phenotype difference. Then the differentially expressed genes (DEGs) and their sequence variations were further analyzed. In total, 21,527 genes were detected, among which 3,965 displayed significant expression variations in tail fat tissues of the two sheep breeds (P < 0.05), including 707 upregulated and 3,258 downregulated genes. Gene Ontology (GO) analysis disclosed that 198 DEGs were related to fat metabolism. In Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, the majority of DEGs were significantly enriched in “adipocytokine signaling,” “PPAR signaling,” and “metabolic pathways” (P < 0.05); moreover, some genes were involved in multiple pathways. Among the 198 DEGs, 22 genes were markedly up- or downregulated in tail fat tissue of Altay sheep, indicating that these genes might be closely related to the fat tail trait of this breed. A total of 41,724 and 42,193 SNPs were detected in the transcriptomic data of tail fat tissues obtained from Altay and XFW sheep, respectively. The distribution of seven SNPs in the coding regions of the 22 candidate genes was further investigated in populations of three sheep breeds with distinct tail phenotypes. In particular, the g.18167532T/C (Oar_v3.1) mutation of the ATP-binding cassette transporter A1 (ABCA1) gene and g.57036072G/T (Oar_v3.1) mutation of the solute carrier family 27 member 2 (SLC27A2) gene showed significantly different distributions and were closely associated with tail phenotype (P < 0.05). The present study provides transcriptomic evidence explaining the differences in fat- and thin-tailed sheep breeds and reveals numerous DEGs and SNPs associated with tail phenotype. Our data provide a valuable theoretical basis for selection of lean-tailed sheep breeds.


INTRODUCTION
Fat tail is a valuable trait that helps sheep (Ovis aries) adapt to harsh conditions, such as extremely cold winters, food shortages, and drought seasons. Fat-tailed sheep characteristically deposit a mass of fatty tissue in the tail region during summer and autumn seasons when nutritious pastures are available. In winter when temperatures are extremely low and grasslands are covered with heavy snow for long periods of time, these animals obtain energy by decomposing the fat deposits in their tails to sustain life (Moradi et al., 2012;Gan et al., 2013a). According to archaeological and genomic findings, modern sheep breeds were domesticated at the Fertile Crescent region of Iraq about 9,000 years ago (Zeder, 1999;Handley et al., 2007;Chessa et al., 2009). Similar to their Asian mouflon ancestor, modern sheep breeds were thin tailed at the early stages. Fat-tailed sheep appeared ∼5,000 years ago through long-term artificial and natural selection during the long evolution process to adapt to harsh local climate conditions (Cooper, 1986;Wade et al., 2010;Muigai and Hanotte, 2013;Beynon et al., 2015). However, in modern society, dietary habits and health concepts have undergone a profound revolution, and mutton with lower fat content is preferred for consumption. Furthermore, the fat deposition requires more energy than the growth of lean tissues. Thus, the efficiency of meat production is higher than that of fat in intensive and semi-intensive sheep industry systems, and the tail fat, which constitute about 20% of the total carcass weight, markedly reduce their economic value (Nejati-Javaremi et al., 2007). For the above reasons, the fat-tailed sheep breeds are gradually becoming less preferable by producers and consumers, and lean-tailed sheep breeds are more desirable (Momani-Shaker et al., 2002;Gokdal et al., 2004;Kashan et al., 2005;Unal et al., 2006;Khaldari et al., 2007;Marai et al., 2009). Therefore, elucidation of the key genes that regulate the deposition and decomposition of tail fat in sheep and the molecular mechanisms controlling fat metabolism would greatly accelerate the breeding course of lean-tailed sheep and production of more healthy mutton to the market.
Fat is not only used for energy storage but is also an important endocrine tissue involved in regulating crucial physiological and biochemical reactions in organisms (Scherer, 2006;Galic et al., 2010;Poulos et al., 2010;Ohashi et al., 2011). Traditionally, research on the molecular regulatory mechanisms of sheep tail adipose tissue development and deposition has only involved the investigation of the functions of single candidate genes. Earlier, Kumar et al. (1998) reported higher expression of the leptin (LEP) gene in perirenal, backside, and omental fat tissues of fat genotype than lean genotype Coopworth sheep. Both activity and expression of lipoprotein lipase (LPL) gene were increased in fat tissue of the fat genotype sheep group. Moreover, its expression was tissue specific but not affected by the nutrition level (Bonnet et al., 2000). In Guangling large-tailed sheep, the uncoupling protein 1 (UCP1) gene was expressed at significantly higher levels in perirenal fat than other tissues but showed very low expression in subcutaneous fat (Yuan et al., 2012). The cell death-inducing DFFA-like effector c (CIDEC) gene was highly expressed in rump fat tissue of Altay sheep, which decreased significantly after a 4week fasting period . Expression of fatty acid binding protein 4 (FABP4) gene was significantly elevated in tail fat of Lori-Bakhtiari, a fat-tailed sheep breed, compared with Zel, a thin-tailed sheep breed (Bakhtiarizadeh et al., 2013). In Altay sheep, FABP4 gene was abundantly expressed in intestinal and rump fat tissues and showed no significant changes when the nutritional status of sheep altered, suggesting its fundamental role in adipose metabolism . The collective findings highlight the involvement of several critical genes in the regulation of tail fat deposition and mobilization in sheep. In general, traits of animals, including tail size, are modulated by multiple interconnected genes that form a refined regulatory network to manage complex internal and external environments. Therefore, research specifically focusing on one or several genes cannot completely elucidate such a network.
In recent years, due to the rapid development of nextgeneration sequencing and application of RNA sequencing (RNA-Seq) technology, numerous genes expressed in tail fat tissue of sheep have been identified via analysis of transcriptome data in attempts to clarify the molecular network regulating adipose deposition and metabolism in fat tail of sheep. RNA-Seq has been used successfully to investigate the genes expressed in fat tissue of sheep. Wang et al. (2014) applied RNA-Seq technology to compare the transcriptome profiles of two sheep breeds, Kazak (fat-tailed) and Tibetan (short-tailed). Their study led to the identification of 646 differentially expressed genes (DEGs) between the two sheep breeds, including 280 upregulated and 366 downregulated genes. Moreover, the genes NEL-like 1 (NELL1) and flavin containing monooxygenase 3 (FMO3), which displayed the most significant fold changes, were highly correlated with adipose deposition in the tail (Wang et al., 2014). Guangling large-tailed sheep and small-tailed Han sheep are two typical fat-tailed sheep breeds in China. Using RNA-Seq, a total of 4,131 DEGs were determined in tail fat tissues of these two breeds, with FABP4, fatty acid binding protein 5 (FABP5), adiponectin (ADIPOQ), and cluster of differentiation 36 (CD36) identified as the four most highly transcribed genes (Li et al., 2018). A research on small-tailed Han and Dorset sheep revealed 602 DEGs, and Gene Ontology (GO) analysis showed that several of these genes were enriched in the triglyceride biosynthetic process (Miao et al., 2015). By applying RNA-Seq, Ma et al. (2018) investigated the tail fat transcriptome of Lanzhou fattailed (long fat-tailed), small-tailed Han (short fat-tailed), and Tibetan (short thin-tailed) sheep and identified several DEGs and long non-coding RNAs (lncRNAs). GO and pathway analysis of DEGs and target genes of differentially expressed lncRNAs revealed that the majority were enriched in fatty acid metabolism and fatty acid elongation-related pathways that contribute to fat deposition (Ma et al., 2018). Previous studies clearly indicate that the tail fat deposition ability of sheep with different tail phenotypes is a complex quantitative trait regulated by multiple genes. The molecular mechanisms underlying tail fat deposition remain to be elucidated.
Due to the multiple functional roles of fat tissues (Scherer, 2006;Galic et al., 2010;Poulos et al., 2010;Ohashi et al., 2011), the fat deposition in the animal's body is most likely to be affected by various factors, such as environment temperature, food, and breeds. Thus, to systematically investigate the relationship between the fat deposition ability and genetic difference of sheep, such influencing factors should be taken into consideration. Altay and Xinjiang fine wool (XFW) sheep are both distributed in the Xinjiang province of China for a long time in history, and their tail fat deposition ability is extremely distinct. Altay sheep is one of the most popular fat-tailed breeds. Their tail and rump are fused together, and they are characterized by their ability to deposit rump fat. Our research showed that the rump fat weight of adult male Altay sheep accounted for about ∼25% of carcass weight on average during autumn (the data have not yet been published). In contrast, XFW sheep, a typical long thin-tailed sheep breed, stores almost no fat tissue in its tail. Thus, these two sheep breeds with distinct tail fat characteristics represent good models for investigating DEGs involved in the regulation of tail fat deposition (Figure 1). Accordingly, we focused on rump and tail fat tissues of Altay and XFW sheep as model animals. RNA-Seq technology was applied to identify DEGs and the associated signaling pathways, with the aim of highlighting candidate genes and mechanisms that play critical roles in regulating adipose deposition in tail of sheep. Genetic variations of these DEGs were also further investigated. Our collective data would provide fundamental information and theoretical guidelines for efficient breeding of lean-tailed sheep.

Ethics Approval and Consent to Participate
All animal experimental procedures were approved by the Biological Studies Animal Care and Use Committee, Xinjiang Production and Construction Corps, Peoples Republic of China, and the ethics committee of Xinjiang Academy of Agricultural and Reclamation Sciences, People's Republic of China approved this study (Permit Number: XAARS-AE-2018012).

Tail Fat Collection and RNA Extraction
Three healthy male Altay and XFW sheep (3-4 years of age) were randomly selected, respectively, from the sheep farm in the Animal Husbandry Institute of Xinjiang Academy of Agricultural and Reclamation Sciences, Shihezi, Xinjiang, China. The body weight of Altay and XFW sheep were 105 ± 5 and 100 ± 5 kg, respectively, and the individuals of each sheep breed were randomly selected to ensure that they did not have a genetic relationship. The sheep were reared according to the feeding standard of meat-producing sheep and goats of China (NY/T 816-2004). At the morning, noon, and afternoon of each day, sheep were provided sufficient forage and clear water. Each evening, an additional 250 g corn per animal was supplied. Six months later, the rump region of Altay sheep were full and round owing to massive fat tissue deposition while no obvious fat tissue were present in tail and rump of XFW sheep. Before sampling, the sheep were anesthetized by injecting 30 mg/kg body weight of pentobarbital sodium (Ningbo, Zhejiang, China) through ear vein (Yan et al., 2013). The pentobarbital sodium is a kind of medium-efficiency barbital hypnotics that could inhibit the uplink activation system of the brain stem reticular structure. Five minutes later, the sheep were slaughtered when they were under anesthesia status, and 100 g of tail fat tissue from each sheep were rapidly collected. The fat tissues were collected from the left rump of Altay sheep, and from the whole tail of XFW sheep due to their very little tail fat. Samples were immersed in liquid nitrogen for transportation and maintained at -80 • C in the laboratory.
Total RNA was extracted from tail adipose tissue using TRIzol reagent (Invitrogen, California, United States) according to the manufacturer's protocol. In short, the fat tissue was homogenized using homogenizer in TRIzol reagent and centrifuged for 5 min at 12,000 × g at 4 • C. The clear supernatant was transferred to a new tube and incubated for 5 min to allow complete dissociation of the nucleoproteins complex. Then the chloroform was added, and the mixture was mixed thoroughly by shaking the tube. Subsequently, the tube was incubated for 2-3 min and centrifuged for 15 min at 12,000 × g at 4 • C. The aqueous phase containing the total RNA was transferred to a new tube, then isopropanol was added and the mixture was centrifuged for 10 min at 12,000 × g at 4 • C; the total RNA precipitate forms a white gel-like pellet at the bottom of the tube. The concentration and integrity of total RNA were evaluated using the 2100 Bioanalyzer instrument (Agilent Technologies, Waldronn, Germany) by measuring absorbance at 260 and 280 nm.

Blood Sample Collection and Genomic DNA Extraction
The fat-rumped Altay sheep, thin-tailed XFW sheep, and short fat-tailed Hu sheep were selected to investigate the distribution of single nucleotide polymorphisms (SNPs) in different sheep populations with distinct tail phenotypes. A total of 104 individuals of each sheep breed were randomly selected, and 5 ml venous anticoagulation blood was collected from each sheep through the anterior vena cava. Genomic DNA was extracted from blood samples using the Blood DNA Extraction kit (Tiangen, Beijin, China). DNA concentration and purity were evaluated using the 2100 Bioanalyzer instrument (Agilent Technologies, Waldronn, Germany). The OD 260 /OD 280 value of DNA samples was ∼1.8. DNA was dissolved in TE buffer (pH = 8.0), and the concentration was adjusted to 200 ng/µl, then stored at -20 • C before being used to detect SNPs in the sequences of candidate genes.

The cDNA Library Construction and Sequencing
Total RNA from three Altay sheep was mixed to generate a cDNA library and similarly from three XFW sheep to generate another cDNA library. After treatment of total RNA with DNase I, poly(A) mRNA was isolated using oligo (dT) magnetic beads (Invitrogen, California, United States). Then the isolated mRNA was cut into short fragments in fragmentation buffer, and the first-strand cDNA was synthesized using random hexamer primers and reverse transcriptase. Subsequently, second-strand cDNA was synthesized, purified using the QiaQuick PCR extraction Kit (QIAGEN, Hilden, Germany), and the poly(A) fragment was added to both ends. The short fragments were connected with sequencing adapters and separated on gels via electrophoresis. Suitable fragments were selected as templates for amplification to construct cDNA libraries. Finally, the two libraries were paired-end sequenced using Illumina HiSeq 2000 at the Beijing Genomics Institute (Shenzhen, China).

Analysis of Sequencing Data
Applying trimmomatic (Bolger et al., 2014), the adapter sequences, reads in which the percentage of low-quality bases (quality value Q ≤ 10) were > 50%, and reads in which unknown bases were >2% were filtered to avoid disturbing the subsequent assembling and analysis (Li et al., 2018). Clean reads were imported into FastQC using FastQ format for further quality control and sequencing quality evaluated (Brown et al., 2017). After that, the clean reads were aligned against the sheep reference genome (Ovis aries, v3.1) using Spliced Transcripts Alignments to a Reference (STAR) (Dobin et al., 2013), and then assembled using StringTie (Kovaka et al., 2019) by comparing with the reference genome of Bos taurus (v3.1) and Ovis aries (v3.1). Meanwhile, GffCompare was used to test the sensitivity and precision of comparison result (Pertea and Pertea, 2020). Finally, the assembled reads were annotated based on protein databases. To get comprehensive annotation information, NR (RefSeq non-redundant proteins) (Pruitt et al., 2007), Swiss-Prot (Boutet et al., 2016), and Cluster of Orthologous Groups of proteins (COG) (Tatusov et al., 2001) databases were used simultaneously.
Expression levels of genes were determined based on gene coverage and fragments per kilobase million (FPKM) values. The gene coverage value represents the ratio of the number of bases on the unigene covered by reads and unigene sequence length and FPKM represents fragments per kilobase of transcript per million reads mapped (Mortazavi et al., 2008;Li et al., 2018). The FPKM method could be used to eliminate the influence of different gene lengths and sequencing levels on calculation of gene expression, allowing direct comparison of gene expression differences between two samples. Expression levels of genes were subsequently normalized and the ratios calculated using Cufflinks 1 . The cuffdiff module of cufflinks was applied to identify DEGs between two samples (Trapnell et al., 2013). The t-test was used to get the P-value of unigenes. The relative expression of unigenes in tail fat tissues of Altay and XFW sheep was determined according to log 2 fat-rump-ratio/log 2 thin-tail-ratio, and DEGs screened taking log 2 ratio ≥ 1 or ≤ -1 and P < 0.05 as standard criteria. Finally, multiple hypothesis testing was applied to revise the P-value of each DEG to guarantee the low false discovery rate in whole. We filtered out DEGs for subsequent analyses by setting the false discovery rate (FDR) as ≤ 0.001 and absolute value of log 2 ratio ≥ 2 as the threshold.
SNPs of different transcriptomes were further detected using the Genome Analysis Toolkit (GATK) software package (v4.0.10) (Franke and Crowgey, 2020). Based on comparison of short sequences, SNP sites were filtered using the Short Oligonucleotide Alignment Program (SOAPsnp) package to obtain high-quality variants and the conditions set as variation detection quality ≥ 30, site depth 10-100 × , and adjacent SNPs > 10 bp (Li et al., 2008).

GO and Kyoto Encyclopedia of Genes and Genomes Analyses of DEGs
GO analysis was employed to annotate functions and further enriching the GO terms of DEGs using David online software and DEGs were mapped to GO terms in the database 2 . Gene numbers for each term were ultimately calculated. Significantly enriched GO terms for DEGs were determined using a hypergeometric test (Cao and Zhang, 2014). Next, the P-value was subjected to Bonferroni correction and threshold P < 0.05 used to define significantly enriched GO terms for DEGs. Meanwhile, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 3 was applied for pathway enrichment analysis of DEGs. Pathways with P < 0.05 were considered significantly enriched for DEGs.

Quantitative Real-Time PCR Validation of RNA-Seq Results
Based on the results from GO, KEGG, and SNP detection of DEGs, 22 candidate genes related to tail fat deposition were selected, including 11 upregulated and 11 downregulated genes.
Using non-pooled RNA samples (n = 3 for each breed), the expression patterns of these genes in tail fat tissues of Altay and XFW sheep were examined via quantitative real-time PCR (qRT-PCR) to validate the reliability of RNA-Seq data. Betaactin (β-actin) gene served as the internal control gene. qRT-PCR primers were designed using Oligo6.0 and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). The 22 genes and their primers are listed in Supplementary Table 1.
Using the PrimeScript II 1st-Strand cDNA Synthesis Kit (Takara, Dalian, China), cDNAs for detection of DEGs involved in tail fat deposition were generated via qRT-PCR conducted on a Roche 480 instrument (Roche, Mannheim, Germany) using SYBR Green PCR Master Mix kit (QIAGEN, Germany) according to the manufacturer's instructions. The reaction mixture comprised 10 µl 2 × Quanti Fast SYBR Green PCR Master Mix, 1 µl cDNA (<100 ng), 0.5 µl forward and reverse primers (10 µM), respectively, and ddH 2 O to a total volume of 20 µl. The following conditions were used for amplification: 95 • C for 5 min, followed by 45 cycles of 95 • C for 10 s, 65 • C for 30 s, and 72 • C for 7 min. qRT-PCR analysis was performed in triplicate for each sample and relative expression for each gene estimated with the 2 − Ct method. Data were analyzed using the Statistical Analysis System (version 6.12; SAS Institute, Inc., Cary, NC, United States) and results expressed as means ± SD. Significance of differences was analyzed using one-way analysis of variance (ANOVA). Differences were considered significant at P < 0.05 and highly significant at P < 0.01.

Analysis of Variations in Candidate Functional Genes Involved in Fat Deposition
Based on SNP analysis of candidate genes, 10 SNPs existing in coding regions of the 22 candidate functional genes involved in tail fat deposition were selected for analysis. Sequences containing the selected SNP sites were exported from transcriptome data and mapped onto the sheep genome (ISGC Oar_v3.1/oviAri3, August 2012) by applying Blat of the UCSC Genome Browser 4 , and subsequently 1,000 bp sequences around the SNP sites were cut out to design primers using Oligo 6.0 (Supplementary Table 2). The PCR fragment sizes were 200-300 bp, and SNP sites were located near the middle of the amplified fragments. Primers were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China).
Genomic DNA of Altay (n = 104), Hu (n = 104), and XFW (n = 104) sheep were amplified using the above primers (Supplementary Table 2), and the distribution of the 10 SNPs in these three breeds were determined via single-strand conformation polymorphism analysis of PCRamplified fragments (PCR-SSCP) or restriction fragment length polymorphism analysis of PCR-amplified fragments (PCR-RFLP) technology. The PCR reaction mixture comprised 2.5 µl 10 × PCR buffer, 0.5 µl DNA template (100 ng), 2 µl dNTPs (2.5 mM each), 0.5 µl forward and reverse primers (10 µM), and 0.5 µl Taq DNA polymerase (5 U/µl) (TaKaRa, Dalian, 4 http://genome.ucsc.edu/ China), with ddH 2 O to a total volume of 25 µl. The following conditions were used for amplification: 95 • C for 5 min; 45 cycles of 95 • C for 30 s, 55 • C for 30 s, and 72 • C for 30 s; and 72 • C for 10 min. Amplified products were detected via 1.5% agarose gel electrophoresis and used for SSCP or RFLP analysis. Then the distributions of SNPs in populations of three sheep breeds were calculated, and the relationship between the SNP and tail fat deposition ability was determined by its genotype frequency. If the genotype frequency of a SNP was significantly different in three sheep breeds (P < 0.05), we speculated it might relate to tail fat deposition ability. The Hardy-Weinberg equilibrium of SNPs in the three groups was verified by calculating the expected frequencies and numbers, and tested using the goodness-of-fit χ 2 test.

Interaction Analysis of Candidate Proteins Involved in Adipose Metabolism
Search Tool for the Retrieval of Interacting Genes (STRING/Proteins) 11.0 (Szklarczyk et al., 2019) was applied to analyze the protein-protein interactions of candidate genes involved in adipose metabolism with the threshold score > 0.4 (medium confidence), and then the interaction network was illustrated.

Summary of Transcriptome Sequencing Data
Two cDNA libraries were constructed using mRNAs extracted from fat-rumped Altay sheep and thin-tailed XFW sheep, sequenced, and two sets of raw reads were obtained containing 51,943,518 and 51,770,440 raw reads, respectively. Lowquality raw reads and adapter sequences were then filtered, ultimately resulting in 46,614,192 and 46,646,110 clean reads. Approximately 84 and 81% clean reads could be mapped to the sheep reference genome (Ovis aries, v3.1). The clean reads were finally assembled into Unigenes, which were categorized to two classes, specifically, clusters and singletons. Clusters were labeled by the prefix "CL, " followed by the cluster id. A single cluster included several Unigenes with > 70% sequence similarity. Singletons were indicated by the prefix "Unigene" (Supplementary Table 3). In total, 153,914 and 117,254 clusters and 78,065 and 56,293 singletons were obtained from the two sample sets, respectively. The mean lengths of clusters were 335 and 317 nt, while mean lengths of singletons were 696 and 629 nt for Altay and XFW sheep groups, respectively (Supplementary Figure 1). Clusters and singletons were further analyzed and filtered, resulting in a final total of 48,894 Unigenes. Transcriptome sequencing data are summarized in Table 1.

Annotation and Expression Analysis of Unigenes
Comparison of the Unigenes obtained with known gene sequences of Bos taurus (v3.1) and Ovis aries (v3.1) revealed a total of 21,527 genes (E-value < 0.00001), which were subsequently matched to NR, Swiss-Prot, KEGG, and COG, leading to 57, 53, 61, and 45% annotation, respectively (Evalue < 0.00001). GO analysis was applied to clarify the biological functions of the above genes (Supplementary Table 3).
Calculation of gene coverage revealed that 65% (13,993/21,527) genes of Altay sheep and 68% (14,638/21,527) genes of XFW sheep had 90-100% coverage (Figures 2A,B). In total, 19,878 annotated genes with FPKM > 0 were detected in the two samples. The FPKM trends of the two samples were comparable, indicating similar expression patterns of the majority of genes in tail fat tissues of Altay and XFW sheep ( Figure 2C). The largest proportions of genes were expressed at low (1 < FPKM < 10) and moderate (10 < FPKM < 100) levels and only a small fraction expressed at high levels (FPKM > 100).
The results indicate that high-throughput sequencing technology has an obvious advantage in detection of low-abundance genes. Further analysis revealed 94.08% (18,701/19,878) of the total genes, including 384 uniquely expressed genes, in the fat rump of Altay sheep and 93.25% (18,536/19,878) of total genes, including 219 uniquely expressed genes, in the thin tail of XFW sheep. Overall, we detected 18,317 common genes in the tail fat tissues of the Altay and XFW sheep ( Figure 2D).

Analysis of DEGs Between Tail Fat Tissues of the Two Sheep Breeds
In total, 8,042 DEGs were identified between the two sheep breeds using FDR ≤ 0.001 and | log 2 ratio| ≥ 1 as filter criteria ( Figure 3A and Supplementary Table 4). Within this gene set, differences in levels of 3,965 DEGs in tail fat tissues of the two sheep breeds were highly significant (FDR ≤ 0.001 and | log 2 ratio| ≥ 2), including 707 highly upregulated and 3,258 highly downregulated genes in fat-rumped Altay sheep, compared to thin-tailed XFW sheep ( Figure 3B and Supplementary Table 4).
Based on significant differences in expression of these genes between tail fat tissues of Altay and XFW sheep breeds and their participation in regulation of fat metabolism, we speculate that the DEGs identified play potentially important roles in influencing tail phenotypes of different sheep breeds. We further focused on 22 DEGs showing highly significant up-or downregulation in tail fat tissue of Altay sheep as candidate genes.

qRT-PCR Validation of RNA-Seq Data
The qRT-PCR expression patterns of the investigated genes were consistent with RNA-Seq data (Figure 4 and Supplementary  Figure 2), supporting the reliability of the expression profile generated with RNA-Seq.

GO and KEGG Analyses of DEGs Between Altay and XFW Sheep
GO was applied for functional analysis of the 8,042 DEGs (Supplementary Table 4). In total, 847 terms were enriched in cellular component, of which 24 terms were significantly enriched, such as "membrane, " "membrane part, " and "cell periphery" (Figure 5A). Overall, 7,014 biological process terms were enriched, 71 of them to a significant extent (P < 0.05), including "cell communication, " "response to stimulus, " and "multicellular organismal process" (Figure 5B). Among the 1,903 terms enriched in molecular function, 32 terms were significantly enriched (P < 0.05), including "substrate-specific transporter activity, " "insulin receptor binding, " and "protein kinase A regulatory subunit binding" (Figure 5C).  The KEGG pathways with P < 0.05 were considered significantly enriched (Supplementary Table 4). In total, 256 pathways were enriched, among which 134 were significantly enriched (P < 0.05). The top 20 pathways, including "MAPK signaling, " "insulin signaling, " "Jak-STAT signaling, " and "phospatidylinositol signaling, " are listed in Figure 5D. We identified 245 and 187 DEGs enriched in the MAPK and insulin signaling pathways, respectively, indicating potential roles of these mechanisms in tail/rump fat metabolism of Altay and XFW sheep.
Other pathways, including "metabolic pathway, " "insulin signaling pathway, " and "glycerolipid metabolism, " additionally play important roles in fat metabolism. The top 30 enriched pathways are presented in Figures 7A,C. Heatmaps were generated and clearly depicted significantly enriched pathways and DEGs (Figures 7B,D) (P < 0.05).
Since the majority of the 198 DEGs were enriched in key pathways significantly related to fat metabolism, we proposed  that these genes are potentially critical for sheep tail phenotype regulation and required to be further investigated.

Detection of Variants in Candidate Genes Related to Tail Fat Metabolism
A total of 41,724 and 42,193 SNPs were detected in tail fat tissue transcriptomes of Altay and XFW sheep, respectively (Supplementary Tables 5,6). We specifically focused on the 22 candidate genes related to tail fat metabolism, which led to the identification of 13 SNPs in the coding regions of nine genes, among which 12 induced amino acid alterations ( Table 3). The distribution of seven SNPs that induced amino acid substitutions were further investigated in Altay, XFW, and Hu sheep populations with different tail phenotypes ( Table 4).
Based on data obtained from 104 individuals of each sheep breed, the distribution of g.18167532T/C (Oar_v3.1) mutation of ABCA1 and g.57036072G/T (Oar_v3.1) mutation of SLC27A2 in these three populations showed significant differences. For the g.18167532T/C mutation of ABCA1, 89.4% individuals in the fat-rumped Altay sheep population were CC genotype, 96.2% of fat-tailed Hu sheep (a breed with a short fat-tailed phenotype between Altay and XFW sheep) were TC genotype, and 86.5% individuals in the long thin-tailed XFW group were TT genotype. The results of the Chi-square test showed that this SNP was not in Hardy-Weinberg equilibrium in XFW and Hu sheep populations (P < 0.01) while the Altay population was in Hardy-Weinberg equilibrium at this site (P > 0.05).
For the g.57036072G/T mutation of SLC27A2, G allele was main genotype in fat-rumped Altay sheep and short fat-tailed Hu sheep populations (63.0 and 68.3% of individuals had the G allele, respectively) while in the thin-tailed XFW sheep population, 89.9% of individuals had the T allele. The three sheep populations were not in Hardy-Weinberg equilibrium at this SNP.
Due to the ABCA1 and SLC27A2 were extremely expressed in tail fat tissues of Altay and XFW sheep, and the distribution of mutations in their gene sequences were significantly different in sheep breeds with different tail phenotype, we hypothesized that they might play potential roles in influencing the tail phenotype in sheep breeds and required to be further investigated.

Transcriptome Studies of Tail/Rump Fat From Sheep
Tail/rump fat deposition in sheep breeds, which distribute in high-latitude regions, facilitate their adaptation to harsh environments. However, these fat-tailed/rumped sheep are undesirable for many modern sheep industry and consumers. Uncovering the mechanisms underlying tail fat deposition would aid in not only improving the meat quality of existing breeds but also selecting the lean-tailed sheep breeds. Meanwhile, the developments in RNA-Seq have facilitated elucidation of the regulatory mechanisms underlying a single trait.
Frontiers in Genetics | www.frontiersin.org (Wang et al., 2014;Miao et al., 2015;Li et al., 2018;Ma et al., 2018), leading to the identification of several important genes and pathways related to fat deposition. However, since the sheep breeds selected for these studies have not always been suitable, the results may not comprehensively reflect differences in the molecular mechanisms regulating tail fat deposition. Guangling large-tailed sheep, Lanzhou fat-tailed sheep, and small-tailed Han sheep studied by Li et al. (2018) and Ma et al. (2018) are fat-tailed breeds that display some differences in tail size. Their transcriptomes are highly similar, and thus limited DEGs have been identified among these breeds. Furthermore, fattailed sheep breeds investigated in these earlier studies are distributed in temperate rural regions of China where the winter is not extremely cold and feed is sufficient over the whole year. Consequently, the speed of deposition and mobilization of tail fat in these breeds may be slower, compared to the breeds in high-latitude areas.
Here, we used two highly suitable sheep breeds, fat-rumped Altay and thin-tailed XFW sheep, for investigating transcriptome differences in tail fat tissue. These two breeds were selected for several reasons. Firstly, both groups of sheep are traditionally localized in the Xinjiang province of China, which avoids the potential effects of distinct environment factors on transcriptome data. Secondly, the tail types are extremely different. Thus, we speculated that there may be significant differences between their tail fat transcriptomes. Thirdly, the living environment of Altay sheep is characterized by long and cold winter and relatively short summer and autumn seasons. The average temperature across the whole year is about 0.7 • C to 4.9 • C and can drop as low as -51.5 • C in winter. To adapt the extreme climate, Altay sheep get the ability to rapidly deposit massive amounts of fat in the rump by consuming sufficient nutritious grass during the warmer seasons and rapidly mobilize in extreme winter conditions to obtain energy and sustain life. Owing to this ability, Altay sheep can live outdoors, even at extreme temperatures of -30 to -40 • C. Thus, the two sheep breeds represent good models to compare transcriptome and DEGs profiles that regulate tail fat deposition.
Using RNA-Seq, a total of 19,878 genes were identified in tail fat tissues, among which 8,042 were differentially expressed between the two sheep breeds (FDR ≤ 0.001 and | log 2 ratio| ≥ 1). Li et al. (2018) reported 5,395 DEGs between tail fat tissues of Guangling large-tailed sheep and small-tailed Han sheep. In addition, 646 DEGs between Kazak and Tibetan sheep were reported by Wang et al. (2014), 390 DEGs between Lanzhou fattailed and Tibetan sheep by Ma et al. (2018), and 602 DEGs between small-tailed Han and Dorset sheep by Miao et al. (2015). The numbers of total genes and DEGs identified in the current study were significantly higher than previously reported figures, which may be attributable to the higher suitability of our animal models.
Among the 3,965 DEGs with significantly different expression (FDR ≤ 0.001 and | log 2 ratio| ≥ 2), 707 were upregulated and 3,258 downregulated in rump fat tissue of Altay sheep. The total number of upregulated genes was significantly lower than that of downregulated genes, similar to the trend observed by the groups of Wang et al. (2014) and Miao et al. (2015). We speculated that the majority of downregulated genes may play important roles in tail fat metabolism in thin-tailed sheep breeds from as early as the time of domestication. Under both human-induced artificial and natural selection, sheep breeds with the powerful ability of efficient tail fat deposition evolved via upregulation of genes promoting fat deposition. Since these types of sheep were more adaptable to harsh environments, their numbers gradually increased and new breeds developed, such as Altay, Guangling large-tailed, Lanzhou fat-tailed, and other fat-tailed sheep breeds worldwide.

DEGs in Sheep Tail/Rump Fat Tissue
DEGs in tail fat tissues of Altay and XFW sheep may play key roles in determining tail phenotype differences. To examine this hypothesis, we focused on the 198 DEGs showing significant differences in expression (FDR ≤ 0.001 and | log 2 ratio| ≥ 2) and closely related to fat metabolism. Among the KEGG pathways of these genes, 134 pathways, which related to adipose metabolism, were significantly enriched, including adipocytokine signaling, PPAR signaling, fat digestion and absorption, and glycerolipid metabolism.
The PPAR signaling pathway regulates cellular differentiation, energy balance, and lipid metabolism (Gross et al., 2017). PPAR exists as α, β, and γ isoforms (Wu et al., 1999), and the expression of PPARγ is necessary for adipocyte differentiation (Farmer, 2006;Floyd and Stephens, 2012). Activation of PPARγ is reported to be essential for deposition of intramuscular fat (Hausman et al., 2018). Previous research indicates that the majority of these upregulated genes are critical for fat deposition while downregulated genes are related to fat mobilization. In bovine mammary glands, mRNA abundance at 60 days postpartum of FABP3 and ACSL1 were 80-and 7-fold greater relative to 15 days antenatal, with peak expression of SLC27A2 and SLC27A6 at 240 and 15 days relative to parturition, respectively, which are significantly associated with milk fat synthesis (Bionaz and Loor, 2008). ANGPTL4 promotes LPL protein intracellular degradation and triglyceride levels in adipocytes (Dijk et al., 2016). PLIN1, an Fsp27 activator, interacts with the CIDE-N domain of Fsp27 and markedly enhances lipid droplet growth by promoting lipid exchange and transfer (Sun et al., 2013). ADIPOQ, an important adipocytokine, modulates glucose and fatty acid oxidation (Abbasi et al., 2004), and its gene sequence polymorphisms are associated with adipose deposition in pig and cattle (Dall'Olio et al., 2009;Choi et al., 2015). ACAA1 and ACADL play critical roles in beta-oxidation of fatty acids (Houten and Wanders, 2010).
AKT, an upregulated gene, has three isoforms in mammals (designated AKT1, AKT2, and AKT3) that are implicated in the regulation of widely divergent cellular processes, such as metabolism, differentiation, proliferation, and apoptosis (Lawlor and Alessi, 2001). In mouse adipocytes, upon rapid activation of AKT2, glucose transporter 4 (GLUT4) translocates to the cell surface and glucose transportation is accelerated (Ng et al., 2008). Most of the downregulated genes were negatively correlated with fat deposition. In adipose tissue of high fat diet-induced obese rats, MAPK8 is significantly downregulated and apoptosis of adipocytes inhibited, which may be the main contributory factor to obesity (Qiu et al., 2010). Mutations of PRKAB1 and PRKAB2 are significantly associated with meat quality traits in pigs (Fontanesi et al., 2003).
Upon further analysis of these DEGs and their pathways, we observed involvement of a number of DEGs in multiple pathways. For instance, ACSL1, ACSL3, ACSL4, and ACSL6 contribute to regulation of adipocytokine signaling, PPAR signaling, fatty acid degradation, and metabolic pathways. MAPK8, AKT2, and AKT3 are involved in MAPK signaling, adipocytokine signaling, and insulin signaling pathway, while FASN and ACACA participate in regulation of fatty acid biosynthesis, insulin signaling, and metabolic pathways ( Table 2). The mechanisms underlying fat metabolism are complex. Fat is not only a tissue used to store energy in animals but also an important endocrine tissue involved in regulating crucial physiological and biochemical reactions (Galic et al., 2010;Ohashi et al., 2011). Accordingly, fat metabolism is regulated by an elaborate network composed of numerous signaling pathways. We speculated that these genes involved in multiple pathways may play bridging roles to connect these signaling mechanisms.
Among the 22 candidate genes, ABCA1, ACACA, and CIDEC were significantly upregulated in rump fat tissue of Altay sheep (P < 0.01) with 5.37, 6.75, and 5.86 times higher expression, compared to tail fat tissue of XFW sheep. C1QTNF1 and HSL were significantly downregulated in rump fat tissue of Altay sheep with 0.15 and 0.13 times expression relative to tail fat tissue of XFW sheep (P < 0.01). ABCA1 is a membrane transporter protein that plays an essential role in the efflux of cholesterol from peripheral tissues back to the liver for participating in lipid metabolism (Zhao et al., 2011). In sheep reared under intensive conditions and offered sufficient feed, ACACA in muscle was significantly upregulated and the fat deposition accelerated (Dervishi et al., 2011). CIDEC (FSP27) located on the surfaces of lipid droplets of adipocytes could promote enlargement or fusion of lipid droplets via clustering and lipid transfer (Li et al., 2017;Slayton et al., 2019;Price et al., 2020). Meanwhile, CIDEC (FSP27) suppressed HSL located on lipid droplet surfaces and inhibited lipolysis (Li et al., 2014). In Altay sheep fasted for 4 weeks, the CIDEC level in rump fat tissue was significantly downregulated . In human liver samples of individuals with obesity and diabetes mellitus, CIDEC was significantly upregulated (Keller et al., 2008;Karki, 2019). The collective studies confirmed essential roles of these candidate genes in fat metabolism.
In view of the significant differences in expression levels of these genes between tail fat tissues of Altay and XFW sheep, and their signaling pathways being closely related to fat metabolism, we speculated that the genes may be potential regulators of tail phenotype differences of sheep and need to be further investigated.

Relevance of Gene Variants and Tail Fat Deposition of Sheep
Previous research has confirmed that a number of gene variants are closely related to tail phenotypes of sheep. Using Ovine SNP50k BeadChip, Moradi et al. (2012) investigated gene variants in two Iranian sheep breeds with fat tail and thin tail phenotypes, respectively. The group identified several mutations that were significantly different between the fat-tailed and thin-tailed sheep breeds in three regions located on chromosomes 5, 7, and X (Moradi et al., 2012). Our group additionally showed that polymorphisms of g.59571364, g.59912586, g.60149273, and g.59383635 loci on Chromosome X are markedly related to tail fat deposition ability of sheep breeds (Gan et al., 2013a,b;Zhang et al., 2013).
RNA-Seq offers novel opportunities for the efficient detection of transcriptome variants (SNPs and short insertion/deletions) in different tissues and species (Cánovas et al., 2010;Cox et al., 2012). Besides whole-genome sequencing, RNA-Seq offers a cost-effective approach for identifying variations and potentially causal mutations underlying the analyzed phenotypes (Hudson et al., 2012;Suárez-Vega et al., 2015). Using RNA-Seq, Suárez-Vega et al. (2017) detected 57,795 variants in the regions harboring quantitative trait loci (QTL) for mild yield, protein, and fat percentages in sheep, among which 21.44% were novel. In the current study, we detected 41,724 and 42,193 SNPs in tail fat tissue transcriptomes of Altay and XFW sheep, respectively, using RNA-Seq. We further focused on SNPs in 22 candidate genes related to tail fat metabolism, 12 of which altered the encoded amino acid ( Table 3).
The distributions of g.18167532T/C mutation of ABCA1 and g.57036072G/T mutation of SLC27A2 were significantly different in the three sheep breed populations with distinct tail phenotypes ( Table 4). ABCA1 encodes a key protein regulating apolipoprotein-mediated efflux of cholesterol and phospholipid from peripheral cells to high-density lipoprotein-cholesterol (HDL-C) (Srivastava, 2002;Fredenrich and Bayer, 2003). Most previous studies have focused on the association of ABCA1 gene polymorphisms with human disease. The SNPs rs4149267, rs1800977, rs1800978, and rs2230806 of ABCA1 are associated with HDL-C concentrations in Caucasian, Sacramento, and French populations (Prochay et al., 2012;Clifford et al., 2013). A significant association was observed between the SNP of ABCA1 and type 2 diabetes in patients of Han Chinese ancestry and Japanese population (Daimon et al., 2005;Kong et al., 2015). A rs2230806 genetic variation of ABCA1 was significantly related to the development and severity of coronary artery disease (CAD) in an Iranian population. Moreover, the K allele of ABCA1 R219K polymorphism has been shown to exert a protective effect against CAD risk and is correlated with decreased severity of CAD, independently of plasma lipid levels (Kolovou et al., 2011;Ghaznavi et al., 2018). Here, we detected a g.18167532T/C mutation of ABCA1. Its distribution in different sheep breeds was significantly related to tail fat deposition ability, which required to be further investigated.
The significance of mutations in SLC27A isoforms has been established in previous reports. Wang et al. (2007) identified a SNP in exon 7 leading to an amino acid alteration in Large White and Meishan pig breeds, which was significantly correlated with growth and carcass traits. SNPs at SLC27A1 and SLC27A4 were associated with saturated fatty acid and stearic acid contents in longissimus dorsi muscle of pig (Melo et al., 2013). Two SNPs in exon 3 and 3'UTR of bovine SLC27A1 exerted effects on milk production traits, such as milk protein and milk fat percentages, in Chinese Holstein cattle (Lu et al., 2010). In the present investigation, the genotype of the g.57036072G/T mutation of SLC27A2 was distinct in three sheep breeds with different tail phenotypes and the G allele was significantly related to fat rump phenotype in Altay sheep ( Table 4).
The two mutations of ABCA1 and SLC27A2 identified in this study led to alterations in the encoded amino acids. In view of their significant association with tail phenotype of sheep, the issue of whether these mutations affect rump fat deposition in Altay sheep by influencing the functions of the corresponding proteins requires further investigation. Here, we confirmed relationships between limited SNPs and tail fat deposition traits of sheep. Further research is warranted to ascertain the associations of several other detected genes with tail fat metabolism.

CONCLUSION
In this research, we examined the differences in transcriptome profiles and sequences of tail fat tissue from Altay and XFW sheep breeds that are distributed within the same province in China but display distinct tail fat deposition traits by applying RNA-Seq followed by qRT-PCR validation to confirm the reliability of our findings. DEGs were identified and their functions evaluated via GO and KEGG analyses. The genes associated with fat metabolism were filtered out for further analysis. Based on the data, 22 candidate genes and two SNPs were identified that potentially contribute to differences in tail fat deposition abilities of sheep. Determination of the specific roles of these DEGs and candidate genes in tail fat deposition ability may aid in the selection of lean-tailed sheep breeds.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material. The raw reads data was submitted to the Short Read Archive (SRA) under the accession number SRR11586694-11586695 and BioProject accession number PRJNA627341 (https://dataview.ncbi.nlm.nih.gov/?search=SUB7313538).

ETHICS STATEMENT
The animal study was reviewed and approved by the Biological Studies Animal Care and Use Committee, Xinjiang Production and Construction Corps, China, and the Ethics Committee of Xinjiang Academy of Agricultural and Reclamation Sciences.

AUTHOR CONTRIBUTIONS
WZ and SG conceived the study, conducted the data analysis, and prepared figures and tables. JW, MX, SW, and JY performed the sample collection and total RNA preparation. WZ, JW, and MX performed the qRT-PCR validation and SNP detection. WZ wrote the manuscript with the assistance of JW and SG. All authors read and approved the final manuscript.