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

Environmental temperature serves as a major driver of adaptive changes in wild organisms. To discover the mechanisms underpinning cold tolerance in domestic animals, we sequenced the genomes of 28 cattle from warm and cold areas across China. By characterizing the population structure and demographic history, we identified two genetic clusters, i.e., northern and southern groups, as well as a common historic population peak at 30 kilo years ago. Genomic scan of cold-tolerant breeds determined potential candidate genes in the thermogenesis-related pathways that were under selection. Specifically, functional analysis identified a substitution of PRDM16 (p.P779L) in northern cattle, which maintains brown adipocyte formation by boosting thermogenesis-related gene expression, indicating a vital role of this gene in cold tolerance. These findings provide a basis for genetic variation in domestic cattle shaped by environmental temperature and highlight the role of reverse mutation in livestock species.

INTRODUCTION

Temperature is one of the most important environmental factors driving evolutionary change in organisms (Parsons, 2005). Mammals require a constant body temperature to ensure 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 non-contractile 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. taurus and B. indicus (Lai et al, 2006; Lei et al, 2006). Archaeological studies support the claim that B. taurus was imported into northern China and northeast Asia from north Eurasia between 5 000–4 000 BP (Cai et al, 2014), and that B. indicus migrated 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 identified PRDM16 as a candidate gene under selection, which is responsible for the modification of BAT function and underpins cold-tolerance in northern cattle.

MATERIALS 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 the B. taurus genome 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: 1118091400014).

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 best K determined using cross-validation (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 (r) 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−8 and g=5, respectively. The mutation rate was estimated using baseml in the PAML package.

Selective sweep analysis

The population-differentiation statistic (FST) (using VCFtools) (Danecek et al, 2011) and nucleotide diversity (Pi) and Tajima’s D (using VariScan v2.0) were estimated using 50 kb sliding windows with a 25 kb step size along each chromosome. Windows in the top 5% of FST values 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, the Q-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 high-capacity cDNA reverse transcription kit (Promega, USA). qRT-PCR 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_10717965; R&D Systems, USA), anti-UCP1 (Cat# ab209483, RRID: AB_2722676; Abcam, UK), anti-PPARγ (Cat# 2430; RRID: AB_823599; CST, USA), anti-HSP90 (Cat# 4874; RRID: AB_2121214; CST, USA) and anti-β-actin (Cat# A5441, RRID:AB_476744, Sigma, USA). The dilution ratio of anti-PRDM16, anti-UCP1, anti-PPARγ and anti-HSP90 was 1:1000 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 ChemiTM 580 (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’s t-test. Statistical significance was set to P<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 of r) were 18.3 kb (r=0.37), 12.9 kb (r=0.26), and 6.3 kb (r=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, including K=2, clearly indicated that the cattle samples could be classified into northern and southern groups (Figure 1E).

Figure 1

Population genetic analysis

A: Geographical distribution of selected cattle in this study. Climate layer (annual mean temperature) with a spatial resolution of 2.5 arc-min was obtained from the WorldClim database (v.2; Fick and Hijmans, 2017; https://worldclim.org/data/worldclim21.html). MG, Mongolia cattle; YB, Yanbian cattle; YN, Yunnan cattle; HN, Hainan cattle. B: Principal component (PC) analysis, PC1 against PC2. C: Neighbor-joining tree of relationships of four cattle breeds. D: Genetic structure of cattle breeds using ADMIXTURE. E: LD decay in HN, YN, southern, and northern cattle. F: Demographic history inferred by PSMC model. Xixiabangma Glaciation (XG, 1 170–800 kya), Naynayxungla Glaciation (NG, 780–500 kya), and Last Glacial Maximum (LGM, ~20 kya) periods are colored.

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 the Naynayxungla 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 different Ne trends. The Ne of 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) and B. 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).

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 of FST values. First, we identified highly differentiated regions using FST, 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, and APOA5) 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 function of BAT in Cpt2A−/− 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 as UCP1 and PGC1-a in response to adrenergic stimulation (Lee et al, 2016). APOA5 treatment can also increase the expression of the UCP1 gene 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, and ALDOB) involved in fatty acid, fructose, and mannose metabolism and associated with signaling pathways, such as the insulin signaling pathway (Supplementary Table S7). Knockout of PDE3B in mice has demonstrated that this gene is involved in the formation of BAT in epididymal WAT depots (Guirguis et al, 2013). ALDOB is 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 (PRDM16 andASXL1) 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).

Figure 2

Selection feature of thermogenic candidate gene

A: Top 10 enriched GO terms in BP (biological process), CC (cellular component), and MF (molecular function). B: KEGG pathway of thermogenesis. Phosphorylation is represented by green arrow lines, expression is represented by black arrow lines, and indirect effect is represented by dashed lines. C: LD analysis of exon SNPs of PRDM16 between northern and southern cattle, and genotype heat map of non-synonymous and synonymous mutation sites in PRDM16 exon. D: Top thermogenic candidate genes under selection. E: Tajima’s D, Pi, and FST values for PRDM16.

Among genes with selective sweep signals, two candidate genes (PRDM16 andCPT2) were involved in thermogenesis; PRDM16 was of the most interest as it is known to increase thermogenesis by promoting the expression of the key gene UCP1 (Seale et al, 2007) (Figure 2B, C). Analysis indicated that there was no strong LD among the PRDM16 SNPs (Figure 2C). PRDM16 had the lowest P-value (P=3.8×10−11) and highest FST (0.52) among genes related to thermogenesis (Figure 2D, E). In addition, although nucleotide diversity (Pi) (0.8×10−3) of PRDM16 was similar to that of other thermogenesis-related genes, Tajima’s D analysis supported the idea that PRDM16 was under selection (D=–1.661) (Figure 2D, E). The PRDM16 genotypes 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).

Figure 3

Genetic polymorphism of PRDM16 across cattle populations

A: SNP tree of PRDM16. B: Percentage of each PRDM16 genotype (c.2336) in northern and southern cattle, respectively. C: Five non-synonymous SNVs in southern cattle and information on PRDM16 779 site in different species. D: Genetic pattern of PRDM16 (c.2336) in cattle genomes worldwide.

Next, we compared the PRDM16 protein sequences to other species (Figure 3C;Supplementary Figure S2), and found that the substitution at Leu779 in the PRDM16 gene 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 in adults (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 the PRDM16 gene is probably related to BAT function, and this locus is likely to play a role in cold tolerance.

Mutation (c.2336 T>C) effects ofPRDM16

To determine the biochemical function of the substitution in PRDM16, 3T3-L1 cells (preadipocyte cell line) ectopically expressing the cattlePRDM16 and PRDM16 MU (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 full differentiation, no differences in morphological characteristics between the PRDM16 and PRDM16 MU 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 the PRDM16 and PRDM16 MU 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 the PRDM16 and PRDM16 MU groups, supporting the idea that PRDM16 loss significantly impedes brown adipocyte differentiation, and PRDM16 overexpression significantly increases brown adipocytes (Seale et al., 2007). Despite the similar differentiation efficiency between the two ectopic PRDM16-overexpressing groups, the mRNA expression levels of four BAT-selective genes (i.e., UCP1, C/EBPβ, PGC1-α, and CIDEΑ) were significantly lower in the PRDM16 MU group than in the PRDM16 group (Figure 4F). Moreover, PRDM16 overexpression increased UCP1 expression to a much greater degree than that found in PRMD16 MU (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

A: Schematic of in vitro differentiation of brown adipocytes. B: mRNA level of PRDM16. C: Oil-red O staining of 3T3-L1 cells at full differentiation. D: mRNA level of PPARγ at full differentiation. E: Protein expression levels of PRDM16 and PPARγ. F: mRNA expression levels of brown fat-selective genes. G: Protein expression level of UCP1. Data are mean±SE. n=5–6/group (A, B, E); n=3/group (C, F). Groups were compared using one-way ANOVA with Tukey post-hoc test. *: P<0.05; **: P<0.01; ***: P<0.001.

DISCUSSION

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. PRDM16 is a key transcriptional regulator in beige adipocyte formation, which stimulates authentic brown fat cells (Seale et al, 2007). In previous research, although PRDM16 was 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 the PRDM16 mutation influences gene function in brown adipogenesis. PRDM16 regulates thermogenic genes by forming complexes with various transcription factors, including C/EBPβ,PGC-1α, PPARα, and PPARγ (Kajimura et al, 2010). Here, although the same differentiation efficiency was induced, suppression of C/EBPβ and PGC-1α mRNA expression levels in the PRDM16 MU group indicated reduced transcription complex formation and thermogenesis-related gene expression, e.g., UCP1, compared to the PRDM16 group. Functional differences in PRDM16 caused by sequence variation could explain why northern cattle are more cold-tolerant than southern cattle. For example, B. indicus may experience higher mortality than B. taurus in cold conditions (Carstens, 1994), possibly due to exhausting their post-natal BAT lipids (Smith et al, 2004). Therefore, on the one hand, well-functioning PRDM16 is required for northern cattle to resist extreme cold, and on the other hand, functional inactivation of PRDM16 impairs 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 AVAILABILITY

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 to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

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.