Population genomics reveals that natural variation in PRDM16 contributes to cold tolerance in domestic cattle

Chun-Long Yan, Jun Lin, Yuan-Yuan Huang, Qing-Shan Gao, Zheng-Yu Piao, Shou-Li Yuan, Li Chen,Xue Ren, Rong-Cai Ye, Meng Dong, Han-Lin Zhang, Hui-Qiao Zhou, Xiao-Xiao Jiang, Wan-Zhu Jin,＊,Xu-Ming Zhou,＊, Chang-Guo Yan

1 College of Agriculture, Yanbian University, Yanji, Jilin 133000, China

2 Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences, Beijing 100101, China

3 University of Chinese Academy of Sciences, Beijing 100049, China

4 North-East Cold Region Beef Cattle Science & Technology Innovation Ministry of Education Engineering Research Center, Yanbian University, Yanji, Jilin 133000, China

5 Annoroad Gene Technology Co.Ltd, Beijing 100176, China

ABSTRACT

Keywords: Population genomics; Cattle; Cold tolerance; PRDM16; Brown adipose tissue

lNTRODUCTlON

Temperature is one of the most important environmental factors driving evolutionary change in organisms (Parsons,2005).Mammals require a constant body temperature toensure optimal biological activity (Haim & Levi, 1990; Hayes &Garland, 1995).This leads to strong selection pressure on the heat production system, including shivering and non-shivering thermogenesis (Cannon & Nedergaard, 2004).Shivering thermogenesis produces heat in the short term (Heldmaier et al, 1989), whereas non-shivering thermogenesis is a noncontractile process that can compensate for the defects of shivering thermogenesis and effectively maintain body temperature (Cannon & Nedergaard, 2004).Although white adipose tissue (WAT) stores excessive energy as triglycerides, brown adipose tissue (BAT), which is activated by cold exposure, is recognized as a major source of adaptive non-shivering thermogenesis (Hughes et al, 2009; Nicholls &Locke, 1984; Rowlatt et al, 1971; Saito et al, 2008).For example, uncoupling protein-1 (UCP1) in BAT dissipates energy into heat through uncoupled respiration, resulting in increased fatty acid oxidation and heat production(Klingenberg, 1999).The thermogenic capacity of BAT is particularly effective for maintaining core body temperature in small mammals and infants (Cannon & Nedergaard, 2004).Nevertheless, the thermogenic program in adipose tissue is a complex transcriptional regulation process that has not been fully dissected.The widely reported transcriptional regulators of adipocytes include peroxisome proliferator-activated receptor-gamma (PPARγ), peroxisome proliferator-activated receptor-gamma coactivator 1α (PGC1-α), Forkhead box C2(FoxC2) and PRD1-BF-1-RIZ1 homologous domain-containing protein-16 (PRDM16) (Kajimura et al, 2010).Among these proteins, PPARγ plays a leading role in the differentiation of all adipocytes (Barak et al, 1999; Nedergaard et al, 2005;Tontonoz et al, 1994).PGC1-α acts together with PPARγ or the thyroid hormone receptor for adaptive thermogenesis(Handschin & Spiegelman, 2006; Puigserver et al, 1998).FoxC2 can increase BAT levels to enhance insulin sensitivity,and PRDM16 can induce the browning of WAT and fibroblasts by driving brown adipogenesis while suppressing white fat adipogenesis (Seale et al, 2007).

Cattle are intimately associated with human civilization and culture.At present, there are about 53 cattle breeds in China,and two recognized species: i.e.,B.taurusandB.indicus(Lai et al, 2006; Lei et al, 2006).Archaeological studies support the claim thatB.tauruswas imported into northern China and northeast Asia from north Eurasia between 5 000–4 000 BP(Cai et al, 2014), and thatB.indicusmigrated from the Indian subcontinent to East Asia around 3 000 BP (Payne & Hodges,1997).Interestingly, the habitats of these cattle and the average annual temperature in which they were domesticated vary widely.Several recent studies have investigated cold adaptation mechanisms in cattle at the genomic level,providing valuable resources for future research (Buggiotti et al, 2021; Ghoreishifar et al, 2020; Hu et al, 2021; Igoshin et al,2021); however, most reported candidate genes/variations lack validation.Here, to detect the molecular footprints underlying cold adaptations in domestic cattle, we sequenced the genomes of 28 cattle, including 14 cold-tolerant cattle lineages (annual average temperature of habitat: 2–6 ℃) and 14 cold-intolerant cattle lineages (annual average temperature of habitat: 20–25 ℃).Through characterization of population history and selective sweeps, we identifiedPRDM16as a candidate gene under selection, which is responsible for the modification of BAT function and underpins cold-tolerance in northern cattle.

MATERlALS AND METHODS

Genome sequencing

We sampled a total of 28 cattle from four different regions in China (i.e., Mongolia, Yanbian, Hainan, and Yunnan).DNA was extracted from the blood of each individual, and degradation was monitored based on its concentration by spectrometry, fluorometry, and 1% agarose gel electrophoresis.Paired-end libraries with an insert size of 150 bp were constructed for each individual and sequenced using the HiSeq X Ten Sequencing System (Illumina, USA).Other cattle genomes were obtained from the NCBI database(Supplementary Table S1).We mapped clean reads after filtering sequencing data to theB.taurusgenome assembly(version ARS-UCD1.2) using BWA v0.7.17 (Li & Durbin,2009).Duplicate reads were removed using Picard tools MarkDuplicates (http://broadinstitute.github.io/picard/).All potential single nucleotide polymorphism (SNP) sites were extracted and filtered using GATK (Mckenna et al, 2010) with HaplotypeCaller.Filtering was performed under the following settings: QD＜2.0, ReadPosRankSum＜-8.0, FS＞60.0,QUAL＜30.0, DP＜4.0, MQ＜40.0, MappingQualityRankSum＜-12.5.ANNOVAR (Wang et al, 2010) and an existing genome annotation file (GFF/GTF) were used to make corresponding annotations on the detected SNPs.All experimental procedures were performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals approved by the State Council of the People’s Republic of China (Document No: 1 118 091 400 014).

Phylogenetic and population structure

Principal component analysis (PCA) was carried out using EIGENSOFT (Price et al, 2006).A phylogenetic tree was constructed from the SNP data using the neighbor-joining method in PHYLIP (Plotree & Plotgram, 1989), and graphical demonstration was performed using Newick Utilities (Junier &Zdobnov, 2010).Population structure was further inferred using ADMIXTURE (Alexander et al, 2009) with component(K) set from 2 to 10 and the bestKdetermined using crossvalidation (CV) analysis.

Linkage disequilibrium (LD) and pairwise sequentially Markovian coalescent (PSMC) analysis

The LD patterns for different breeds were calculated using the squared correlation coefficient (r2) between pairwise SNPs with PopLDdecay script (https://github.com/BGI-shenzhen/PopLDdecay).The PSMC model (https://github.com/lh3/psmc)parameters were set to: -N25 -t15 -r5 -p "4+25＊2+4+6", and mutation rate and generation time were set to: μ=1.1×10-8and g=5, respectively.The mutation rate was estimated usingbasemlin the PAML package.

Selective sweep analysis

The population-differentiation statistic (FST) (using VCFtools)(Danecek et al, 2011) and nucleotide diversity (Pi) and Tajima’sD(using VariScan v2.0) were estimated using 50 kbsliding windows with a 25 kb step size along each chromosome.Windows in the top 5% ofFSTvalues were selected as candidate windows to obtain corresponding candidate genes.Fisher’s exact test was performed on synonymous and non-synonymous SNPs in the exon region using PLINK v1.9 (Purcell et al, 2007) to determine the final candidate genes.Before this step, PLINK v1.9 was used to remove sites with strong LD correlation (--indep-pairwise 50 5 0.5), and non-synonymous sites were used for Fisher’s exact test (--fisher).Finally, theQ-value was calculated using the R package fdrtool, and the site with q＜0.01 was selected as the candidate locus to obtain corresponding candidate genes.Enrichment analysis was conducted using gprofiler2 (Kolberg et al, 2020).

Cell culture

Lentiviruses with PRDM16 variants were produced by transfecting HEK293T cells with core plasmids and two helper plasmids (psPAX2 and pMD2G).The transfections were implemented using the polyethylenimine (PEI) method at a PEI:core plasmid:psPAX2:pMD2G ratio of 27:4:3:2.The medium was changed 4–6 h after transfection.After 48 h, the virus-containing medium was harvested and filtered.The 3T3-L1 cells were then incubated overnight (37 ℃, 5% CO2) with the viral supernatant and 8 μg/mL polybrene.For browning differentiation, confluent 3T3-L1 cells were incubated for 2 days in a brown adipogenic induction cocktail (Dulbecco’s Modified Eagle Medium (DMEM) containing 10% fetal bovine serum (FBS), 20 nmol/L insulin, 1 nmol/L 3,3,5-triiodo-L-thyronine (T3), 0.5 mmol/L isobutylmethylxanthine, 0.125 μmol/L indomethacin, and 1 mmol/L dexamethasone).The cells were then maintained in differentiation medium (DMEM containing 10% FBS, 20 nmol/L insulin, and 1 nmol/L T3) for 6 days (37 ℃, 5% CO2).The induction medium was changed every 2 days.At day 8, fully differentiated brown adipocytes were applied for all experiments in this study.

RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA from tissues and cells was extracted with Trizol reagent (Thermo Fisher Scientific, USA).Reverse transcription of 2 μg of total RNA was performed with a highcapacity cDNA reverse transcription kit (Promega, USA).qRTPCR was performed with a SYBR Green Master Mix(Promega, USA) and detected using a Prism VIIA7 Real-Time PCR System (Applied Biosystems, USA).Primers were designed using Primer Quest (Integrated DNA Technologies,USA).Primer sequences are provided in Supplementary Table 2.

Western blot analysis

Cells were lysed in RIPA buffer containing 150 mmol/L sodium chloride, 1.0% TritonX-100, 0.5% sodium deoxycholate, 0.1%sodium dodecyl sulfate (SDS), and 50 mmol/L Tris with freshly added protease and phosphatase inhibitor cocktail (Roche Diagnostics Corp, USA).Equal amounts of protein were distributed in 10% SDS-polyacrylamide gel After electrophoresis, the proteins were transferred to a polyvinylidene fluoride (PVDF) membranes, incubated with blocking buffer (5% fat-free milk) for 1 h at room temperature,and blotted with the following antibodies overnight (4 ℃): anti-PRDM16 (Cat# AF6295, RRID:AB_10 717 965; R&D Systems,USA), anti-UCP1 (Cat# ab209483, RRID: AB_2 722 676;Abcam, UK), anti-PPARγ (Cat# 2 430; RRID: AB_823 599;CST, USA), anti-HSP90 (Cat# 4 874; RRID: AB_2 121 214;CST, USA) and anti-β-actin (Cat# A5441, RRID:AB_476 744,Sigma, USA).The dilution ratio of anti-PRDM16, anti-UCP1,anti-PPARγ and anti-HSP90 was 1:1 000 and the dilution ratio of anti-β-actin was 1:10000.The membranes were then incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 h at room temperature.Signals were visualized using a Mini ChemiTM580 (Sage Creation Science, China) with Super Signal West Pico Chemiluminescent Substrate (Pierce, USA).

Statistical analysis

Data are expressed as mean±standard error (SE).Comparisons between groups were performed with one-way analysis of variance (ANOVA) or Student’st-test.Statistical significance was set toP＜0.05.

RESULTS

Genome sequencing and population history

Whole-genome sequencing of 28 cattle with an average depth of 33.66× obtained 17.3 billion clean reads (Figure 1A;Supplementary Figure S1 and Table S3).In total, 45.2 million single nucleotide polymorphisms (SNPs) were identified, most of which were located in the intergenic (61.51%) and intron(35.75%) regions (Supplementary Table S4).Neighbor-joining trees and PCA based on total SNPs clustered the cattle into two main groups: i.e., northern and southern groups(Figure 1B, C).The first principal component (PC1),representing 32.41% of total variation, separated the samples into northern and southern cattle (Figure 1C).We further analyzed the genomes and found that the rates of LD decay were greater in the southern cattle than in the northern cattle.Half distances (half ofr2) were 18.3 kb (r2=0.37), 12.9 kb(r2=0.26), and 6.3 kb (r2=0.27) for the northern (Mongolia: MG and Yanbian: YB) cattle, Hainan (HN) cattle, and Yunnan (YN)cattle, respectively (Figure 1D).ADMIXTURE analyses with different component (K) values, includingK=2, clearly indicated that the cattle samples could be classified into northern and southern groups (Figure 1E).

The demographic history of cattle was determined using the PSMC model (Li & Durbin, 2011).Results showed two expansions and two bottlenecks, with population peaks at ～50 and ～700 kilo years ago (kya) and population bottlenecks at～30 and 400 kya, respectively (Figure 1F).There were two sharp declines in population, which both occurred during the glacial period (Naynayxungla Glaciation and Last Glacial Maximum), consistent with the idea that environmental temperature has a determinable impact on population size.Similar historical patterns have been reported in many other mammals, such as the giant panda, yak, and snub-nosed monkey (Qiu et al, 2015; Zhao et al, 2013; Zhou & Pawlowski,2014).Global glaciations are the most probable cause of sudden change in the global climate and can directly affect species populations.Indeed, we found that after theNaynayxungla Glaciation (780-500 kya), northern cattle experienced a long-term bottleneck period until 70 kya.In contrast, the effective population size (Ne) of southern cattle recovered rapidly after the Naynayxungla Glaciation(Figure 1F), consistent with previous studies (Chen et al,2018; Lan et al, 2018; Mei et al, 2018); this could be explained by the improved living environment in southern areas during glaciation (Murray et al, 2010).At ～60 kya, HN and YN cattle showed differentNetrends.TheNeof HN cattle increased rapidly (Figure 1F), likely due to the geographical location of Hainan, a small and comparatively isolated island that lacks natural predators, which promoted the survival and reproduction of cattle.According to mitochondrial DNA haplotypes,B.taurus(northern cattle) andB.indicus(southern cattle) were both derived from extinct wild aurochs(B.primigenius), with divergence between the two species dating back 250 kya (Bradley et al, 1996).

Figure 1 Population genetic analysis

Genomic scan of selective sweeps

To identify genetic modifications that occurred under different temperatures, we analyzed selective sweeps between the cattle groups: i.e., northern (MG and YB) and southern (YN and HN) cattle.Selective sweep analysis was performed for whole genomes based on the distribution ofFSTvalues.First,we identified highly differentiated regions usingFST, and then determined the top 5% in 50 kb windows with 25 kb steps.Final candidate genes were then determined and ranked using Fisher’s exact test (q＜0.01).In total, 197 candidate genes were identified with strong selective sweep signals(Supplementary Table S5).The most significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of the candidate genes (e.g.,SCP2,Cpt2, andAPOA5) was the PPAR signaling pathway (P=1.6×10-2) (Supplementary Table S6).SCP2 expression significantly alters the structure of lipid droplets (Atshaves et al, 2001) and affects the functionof BAT inCpt2A-/-mice, thereby hindering their ability to adapt to temperature changes (Lee et al, 2015).Furthermore,Cpt2A-/-interscapular BAT fails to induce the expression of thermogenic genes such asUCP1andPGC1-ain response to adrenergic stimulation (Lee et al, 2016).APOA5treatment can also increase the expression of theUCP1gene in adipocytes(Zheng et al, 2017).Furthermore, many fatty acids positively affect thermogenesis by activating BAT (Heeren & Scheja,2018; Li et al, 2018; Quan et al, 2020; Takato et al, 2017).We also found many candidate genes (e.g.,PDE3B,CPT2, andALDOB) involved in fatty acid, fructose, and mannose metabolism and associated with signaling pathways, such as the insulin signaling pathway (Supplementary Table S7).Knockout ofPDE3Bin mice has demonstrated that this gene is involved in the formation of BAT in epididymal WAT depots(Guirguis et al, 2013).ALDOBis involved in insulin biosynthesis and secretion, as well as insulin receptor signaling (Gerst et al, 2018).Insulin pathways and fat metabolism are inseparable and can affect the development of BAT, leading to obesity and insulin resistance (Lynes et al,2015; Montanari et al, 2017).Consistently, in our study, Gene Ontology (GO) enrichment analysis revealed two candidate genes (PRDM16andASXL1) related to fat cell differentiation(GO:0045598), brown fat differentiation (GO:0050873), and white fat cell differentiation (GO:0050872) (Figure 2A;Supplementary Tables S8, S9).

Among genes with selective sweep signals, two candidate genes (PRDM16andCPT2) were involved in thermogenesis;PRDM16was of the most interest as it is known to increase thermogenesis by promoting the expression of the key geneUCP1(Seale et al, 2007) (Figure 2B, C).Analysis indicated that there was no strong LD among thePRDM16SNPs(Figure 2C).PRDM16had the lowestP-value (P=3.8×10-11)and highestFST(0.52) among genes related to thermogenesis(Figure 2D, E).In addition, although nucleotide diversity (Pi)(0.8×10-3) ofPRDM16was similar to that of other thermogenesis-related genes, Tajima’sDanalysis supported the idea thatPRDM16was under selection (D=–1.661)(Figure 2D, E).ThePRDM16genotypes found in the northern and southern cattle were well distinguished and consistent with the phylogenetic tree created using the SNPs of this gene(Figure 3A).We discovered five non-synonymous single nucleotide variants (SNVs), one of which (c.2336 T＞C,p.L779P) was found at a higher level (93%) in southern cattle than in northern cattle (Figure 3B, C; Supplementary Table S10).

Next, we compared the PRDM16 protein sequences to other species (Figure 3C;Supplementary Figure S2), and found that the substitution at Leu779in thePRDM16gene in northern cattle was the same as that in species with complete BAT function (e.g., mouse, rat, and hamster) (Figure 3C).In rodents, BAT is intact and persists throughout their lifetime,and thermogenesis activity is complete (Kirov et al, 1996;Scarpace et al, 1994).However, in many large mammals,such as humans and sheep, BAT function is available during infancy but can only be activated under certain conditions inadults (Lidell et al, 2013; Nahon et al, 2020).Conversely, the proline substitution in southern cattle was the same as that in species with incomplete or null BAT function (e.g., sheep, pig,whale, horse, platypus, elephant, sirenian, marsupial, human,and rabbit) (Figure 3C).Moreover, we explored the genetic pattern of these substitutions (c.2336 T＞C, p.L779P) across cattle genomes worldwide, and found that cattle in cold regions had a higher frequency of the c.2336 C＞T mutation,consistent with the pattern in China (Figure 3D).Thus, we hypothesized that the substitution of residue 779 in thePRDM16gene is probably related to BAT function, and this locus is likely to play a role in cold tolerance.

Figure 2 Selection feature of thermogenic candidate gene

Figure 3 Genetic polymorphism of PRDM16 across cattle populations

Mutation (c.2336 T＞C) effects of PRDM16

To determine the biochemical function of the substitution inPRDM16, 3T3-L1 cells (preadipocyte cell line) ectopically expressing the cattlePRDM16andPRDM16MU (c.2336 T＞C, L779P mutation ofPRDM16) coding sequences were generated and induced to differentiate towards beige adipocytes (Figure 4A).The overexpression efficiency was kept at equivalent levels (Figure 4B, E).After fulldifferentiation, no differences in morphological characteristics between thePRDM16andPRDM16MU groups were observed (Figure 4C).In addition, we did not find significant differences in the mRNA and protein expression levels of PPARγ, a key adipogenesis-regulating gene, between thePRDM16andPRDM16MU groups (Figure 4D, E).However,the differentiation efficiencies of PPARγ mRNA and protein expression were lower in the control group (cells infected with an empty vector) than in thePRDM16andPRDM16MU groups, supporting the idea thatPRDM16loss significantly impedes brown adipocyte differentiation, andPRDM16overexpression significantly increases brown adipocytes(Seale et al., 2007).Despite the similar differentiation efficiency between the two ectopicPRDM16-overexpressing groups, the mRNA expression levels of four BAT-selective genes (i.e.,UCP1,C/EBPβ,PGC1-α, andCIDEA) were significantly lower in thePRDM16MU group than in the PRDM16 group (Figure 4F).Moreover,PRDM16overexpression increasedUCP1expression to a much greater degree than that found inPRMD16MU (Figure 4F, G).These results indicate that the L779P mutation significantly impaired normal PRDM16 function in the formation of brown adipocytes in southern cattle, which live in warmer areas relative to northern cattle.

Figure 4 PRDM16 779P allele reduced brown adipogenesis

DlSCUSSlON

We compared the whole genomes of northern and southern cattle in China, which live in extremely cold and warm environments, respectively.We identified a total of 197 candidate genes with selective sweep signals.However, these genes should be subjected to further validation given the many challenges in accurate detection of selective sweeps across genomes.For example, the current methodology could be confounded by many processes, such as recombination and drift, and the effects of changing demography over time(Horscroft et al, 2019).Nevertheless, we found that one candidate gene,PRDM16, is a forceful genome effector that facilitates cold adaptation.PRDM16is a key transcriptional regulator in beige adipocyte formation, which stimulates authentic brown fat cells (Seale et al, 2007).In previous research, althoughPRDM16was introduced before cell differentiation, nearly all adipocytes were activated to express BAT-selective genes (Seale et al, 2007).In this study, we found that BAT-selective genes were up-regulated in PRDM16-overexpressing 3T3-L1 cells compared to controls,indicating that thePRDM16mutation influences gene function in brown adipogenesis.PRDM16 regulates thermogenic genes by forming complexes with various transcription factors,includingC/EBPβ,PGC-1α,PPARα, andPPARγ(Kajimura et al, 2010).Here, although the same differentiation efficiency was induced, suppression ofC/EBPβandPGC-1αmRNA expression levels in thePRDM16MU group indicated reduced transcription complex formation and thermogenesis-related gene expression, e.g.,UCP1, compared to thePRDM16group.Functional differences inPRDM16caused by sequence variation could explain why northern cattle are more cold-tolerant than southern cattle.For example,B.indicusmay experience higher mortality thanB.taurusin cold conditions (Carstens, 1994), possibly due to exhausting their post-natal BAT lipids (Smith et al, 2004).Therefore, on the one hand, well-functioningPRDM16is required for northern cattle to resist extreme cold, and on the other hand, functional inactivation ofPRDM16impairs beige adipocyte formation,which is beneficial for the environmental adaptability of southern cattle.These findings help improve our understanding of adaptive genetic variations in cattle and other livestock species living in different temperature regions.

DATA AVAlLABlLlTY

This whole-genome shotgun project was deposited in the NCBI under BioProjectID PRJNA737584 and in GSA under accession No.subCRA008925 and in Science Data Bank under DOI: 10.11922/sciencedb.01524.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETlNG lNTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRlBUTlONS

C.G.Y., X.M.Z., and W.Z.J.designed the research, analyzed data, and revised the manuscript.C.L.Y., J.L., and Y.Y.H.performed experiments, analyzed data, and wrote the manuscript.Q.S.G.and Z.Y.P.collected samples, performed experiments, and analyzed data.S.L.Y.and X.R.analyzed data.L.C.revised the manuscript.R.C.Y, M.D., H.L.Z., H.Q.Z.,and X.X.J.collected samples.All authors read and approved the final version of the manuscript.

ACKNOWLEDGEMENTS

We thank Dr.Inge Seim for value suggestions and comments.