A total of 1,350 tea accessions were collected at diverse sites from Lincang, Yunnan, and other regions of China during the period from 2019 to 2020. Among them, 31 samples were collected in Cangyuan, 115 in Fengqing, 78 in Gengma, 235 in Linxiang, 291 in Shuangjiang, 187 in Yongde, 217 in Yun, 81 in Zhenkang, and 115 in Others. The Others group was used as a control for the samples of the eight geographic populations, so the sampling range of the Others group was widely dispersed (Figure 1: black dot). Further, the accession KM6 (Camellia Cuspidata) was collected to use as an outgroup during phylogenetic analysis. The details about the sampled populations are presented in Supplementary Table 1 and the geographical distribution of these points is depicted in Figure 1.

The total DNA was extracted using the DNA secure plant kit (TIANGEN, Beijing), following the manufacturer’s protocol. Around 2 μg of the extracted DNA was used to construct the sequencing library for each accession using the NEBNext Ultra DNA Library Prep Kit (NEB Inc., America), following the manufacturer’s instructions. Paired-end sequencing libraries with an insert size of approximately 400 bp were sequenced on the Illumina NovaSeq 6000 platform. We used fastp v0.12.2 (Chen et al., 2018) for the removal of adaptor contamination, poly-N, and low-quality reads (reads having >40% bases and Phred score ≤20). Further, the paired-end reads with sequence lengths below 70 bp were filtered out. Thus, only high-quality cleaned reads were retained for downstream analysis.

Paired-end reads were mapped to the reference genome of Camellia sinensis (Shuchazao; Wei et al., 2018) through BWA v0.7.17-r 1,188 (Li and Durbin, 2009) using default parameters. Conversion of SAM to BAM and exclusion of unmapped and multi-mapped reads were performed through SAMtools v1.3.1 (Li et al., 2009). Further, the duplicated reads were filtered out using Picard v2.1.1.¹

After BWA alignment, the reads around indels were realigned. Realignment was performed with GATK 3.3-0-g37228af² (McKenna et al., 2010) in two steps. First, we used the RealignerTargetCreator package to identify regions where realignment was needed, followed by a realignment of reads to these regions using IndelRealigner, and created a realigned BAM file for each accession. Then we detected the variation of each sample and obtained the original variation set file (gVCF format) through GATK Haplotype Caller, and gVCF files were further integrated to obtain population variation data.

The SNP filter expression parameters were set as: QD<2.0 || MQ<40.0 || FS>60.0 || SOR>5.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0 || QUAL <30. The InDel filter expression parameters were set as: QD<2.0 || ReadPosRankSum < −20.0 || InbreedingCoeff < −0.8 || FS>200.0 || SOR>10.0 || QUAL<30 (DePristo et al., 2011). Only insertions and deletions shorter than or equal to 40 bp were considered. Indels and SNPs with none bi-allelic, >50% missing calls and MAF < 0.005 were removed, which yielded the basic set. SNPs with MAF < 0.05, none bi-allelic, >50% missing calls were further removed for phylogenetic tree structure, genetic diversity analysis, LD decay, PCA and population structure analyses (the core set). The annotation of SNPs and InDels was performed through ANNOVAR v2015-12-14 (Wang et al., 2010) using tea genome as a reference.

The populations were clustered to assess the pattern of variation among the sampled populations. We used the whole-genome SNPs to construct the maximum likelihood (ML) phylogenetic tree with 100 bootstraps using SNPhylo v20140701 (Lee et al., 2014). Camellia cuspidate was used as an outgroup. Color coding of the phylogenetic tree was done through the iTOL web server.³

The SNPs in LD were filtered out using PLINK v1.90b3.38 (Purcell et al., 2007) with a window of size 50 SNPs (advancing 5 SNPs at a time) and an r² threshold of 0.5 to determine a pruned SNP set to be used in the population structure analysis. LD-based pruning reduces the effects of ascertainment bias in a relatively efficient manner (Malomane et al., 2018). Principal component analysis (PCA) was performed with the Genome-wide Complex Trait Analysis (GCTA) v1.25.3 (Yang et al., 2011), and the first three eigenvectors were plotted. LD was calculated using PopLDdecay v3.41 (Zhang et al., 2019). The pairwise r² values within and between different chromosomes were calculated. The LD for each group was calculated using SNP pairs only from the corresponding group.

The population structure was analyzed using the ADMIXTURE v1.3 (Alexander et al., 2009) program with a block-relaxation algorithm. To explore the convergence of individuals, we predefined the number of genetic clusters K, from 2 to 9 and ran the cross-validation (CV) error procedure. Default methods and settings were used in the analyses.

The primary genetic diversity parameters like observed heterozygosity (H_O), expected homozygosity (H_E), inbreeding coefficient (F), the average pairwise diversity within a population (θπ), and Tajima’s D were calculated using the vcftools v0.1.13 (Danecek et al., 2011) with 100 kb sliding windows. In addition, overall genetic differentiation across populations measured by Weir and Cockerham’s estimator of F_ST (Weir and Hill, 2002) was also calculated using the same software.

To assess the potential correlation of environmental and geographic variation with the tea genetic structure, the use of variables was necessary to adequately capture the general environmental and geographical differences among the Lincang regions sampled. For this purpose, we calculated the pairwise genetic distance matrix using the PLINK (v1.90b3.38). Meanwhile, we downloaded 4 environmental variables and two geographic variables available at each ecotype site (resolution 30 arc seconds) through WorldClim⁴ (Hijmans et al., 2005).

Among the four environmental variables, we selected two temperature and two precipitation variables. The temperature variables (in °C*10, expressed to the nearest tenths) included the mean diurnal range (MDR) and the mean temperature of the wettest quarter (MTW). The precipitation variables (expressed to the nearest mm) included annual precipitation (AP) and precipitation of the wettest month (PWM; Supplementary Table 10). In addition, we derived 10 additional variables from the squares (4) and cross-products (6) of the original four environment variables to explore whether non-linear environmental effects may also affect the genetic structure.

The geographic variables included LONG and LAT. These variables tended to show lower correlations with each other. We intended to reduce the redundancy inherent among highly associated variables (Hall and Beissinger, 2014). For the geographical distance matrix, we calculated the great-circle distance (the closest distance between two points on the Earth’s surface) in miles from untransformed LONG and LAT values for each pair of accessions using the “geosphere” package in R.⁵ For the environment distance matrix, we calculated an environmental distance matrix from Euclidean distances between pairs of accessions using all 14 environmental variables.

To sort out the potential effects of environmental variables and geographical isolation on genetic variation in tea tree samples, the Mantel test (Diniz-Filho et al., 2013) of correlation among genetic, environmental and geographical distance matrices was applied by using the VEGAN⁶ package in R, in which significance testing of the correlations was performed with 10,000 permutations. Then, we explored whether there was a significant association between environmental distance and genetic distance by using a partial Mantel test, adjusted for any effects of geographical distances. Significance in this test was interpreted as meaning that genetic variation among the tea was influenced by environmental selection, whereas a nonsignificant result suggested a role for isolation by distance (genetic drift).

The historical relationship between Lincang tea geographical populations was estimated using TreeMix (Pickrell and Pritchard, 2012), which uses a Maximum Likelihood (ML) method based on a Gaussian model of allele frequency change. The topology of the ML trees changes depending on the number of migration events (m) allowed in the model. Here we use m = 1 to m = 5 (Fitak, 2021). The bootstrap values on the tree are based on 1,000 replicates. Arrows on the graph represent admixture events between different tea populations. The Other tea population was used for roots.

To uncover the evolution history of lincang tea subpopulations, we use MSMC2⁷ to infer population size of each group. The input files for MSMC2 were generated according to MSMC Tools.⁸ In brief, only sites with uniquely mapped reads and sites with coverage depths between 0.5-fold and 2-fold of mean depth were used in the analyses. The remaining genomic regions were masked using the script bamCaller.py. Then all segregating sites within each group were phased using SHAPEIT (Version: v2.r904; Delaneau et al., 2011). A mutation rate of 6.1 × 10^–9 per site per year was used.

RAiSD (Raised Accuracy in Sweep Detection, Version 2.9; Alachiotis and Pavlidis, 2018) was used to detect signatures of selective sweeps based on the μ statistics. The significant threshold for μ statistic score was set as top 0.1%. Then, We performed a genetic differentiation (F_ST) and nucleotide polymorphism (θπ) based cross approach to investigate the selection signals across the whole genome. A 50 kb sliding window with 10 kb step approach was applied to quantify F_ST and θπ by using VCFtools software (v0.1.13). The annotated genes living in these regions were considered candidate selected genes.

Resequencing of 1,350 tea samples yielded around 9.67 Tb data (64,462,516,344 paired-end raw reads). Mapping of these reads with the reference genome resulted in an average alignment rate of 98.9 ± 2.65% (63.79–99.54%, Supplementary Figure 1 and Supplementary Table 1). The information about 9 populations with reads is presented in Table 1. The average alignment rate of the population ranged from 96.75 to 98.75%; the average number of reads per sample of the population ranged from 45,974,894 to 56,685,556;

After applying basic filtering criteria (see “Materials and Methods” section), we identified 356,171,898 SNPs and 27,367,688 short genomic insertions and deletions (indels). The result of filtering is presented in Supplementary Figures 2, 3 and Supplementary Tables 2, 3. Further filtering yielded a core set of 27,550,879 SNPs and 1,139,750 indels (≤40 bp) with minor allele frequency (MAF) more than 0.05 and max missing less than 0.5. The information about SNPs and InDels is presented in Table 1 and Supplementary Table 4. Meanwhile, in the core data set, 91.36% (25,170,861) of the SNPs were located in the intergenic region, 4.82% (1,327,890) in the intronic region, 1.31% (36,049) in the 5′-UTR, and 2.95% (81,272) in the 3′-UTR. We observed 0.91% (251,002) and 1.01% (277,265) of the SNPs in the upstream and downstream regions of the genes, respectively. Further, 0.0151% (4,150) were located in the variable splicing region and 1.44% (395,863) were present in the exonic regions of the genes (Table 2 and Supplementary Figure 4). Among the exonic SNPs, the proportion of non-synonymous SNPs was found to be 56.31% (223,048) and that of synonymous SNPs was found to be 41.60% (164,772), with the non-synonymous to synonymous mutation ratio of 1.353. The total number of stop-gain SNP mutations was 7,652, whereas the total number of stop-loss SNP mutations was 557 (0.00202%). In addition, 950,567 indels were found to be located in the intergenic region (83.40%), followed by the intronic region (Intronic) with 118,677 (10.41%; Table 3 and Supplementary Figure 4).

The SNPs after filtering linkage disequilibrium sites (MAF > 0.05) were used to analyze the population structure and differentiation. According to the calculated CV error value, when K = 3, the CV error value is the smallest (Figure 2A). Therefore, there are three genetic stocks (genetically different populations) represented as red, green, and blue color (Figure 2B) in Lincang. Interestingly, Shuangjiang and Zhenkang populations were found to be almost genetically pure, whereas the other populations showed a substantial level of genetic admixtures (Figure 2B). The red genetic stock appears to be the most dominant followed by green, whereas the third genetic stock, blue, represents a very small proportion only in the Yun, Zhenkang, and Fengqing tea populations.

In the principal component analysis of tea accessions from Lincang (Figure 2), the three principal components of PC1, PC2, and PC3 represented 5.99, 2.28, and 0.61% of the total genetic variance, respectively. The accessions were grouped into three subgroups based on the STRUCTURE inferred clustering result with 1,105, 172 and 74 accessions came to gather for sub-populations 1, 2, and 3, respectively (When an individual has the highest proportion of red ancestors, it is classified as group 1; then, accessions with the highest proportion of yellow ancestors were assigned into the group 2; otherwise, group 3; Figure 2B). Among them, the first PC (PC1) distantly clustered varieties from groups 2 and the combination of PC1 (5.99%) and PC2 (2.28%) can distinguish groups 2 (red, Figure 2C) from other accessions.

Based on the core SNP set, the phylogenetic tree was constructed using the maximum likelihood (ML) method with KM6 (C. cuspidata) as the outgroup, and the bootstrap value of 100. The result of regional phylogenetic analysis showed that no precise geographic or regional clustering was observed in the phylogenetic tree (Supplementary Figure 5). However, when tea tree samples were analyzed into three subgroups based on the results of genetic stratification analysis (Figure 3A), the phylogenetic tree showed four distinct clustering situations, which recapitulates the same patterns in the principal component analysis (PCA) and model-based clustering (Figure 3B). Among them, group 1 was divided into two clusters, and group 2 and 3 were divided into a cluster. In group 1 (Figure 3A; blue), the samples collected in Cluster 1 were dominated by Yun and Linxiang, which had a close genetic relationship with the outgroup. It can be inferred that the border between Yun and Linxiang may be the origin of tea. The samples in Cluster 2 are dominated by Shuangjiang and Linxiang, which are genetically more distant from the outgroup, and it can be inferred that the area covered by Cluster 2 may have been introduced later and thus developed into numerous new branches. In group 2 (Figures 3A,B), tea samples from Yun County were dominant, and the tea samples collected in Yun (leaf nodes) were genetically closer to the proximal ancestors (the inner node). It can be further inferred that the region of Yun is the area where the initial origin has produced differentiation afterward. In group 3 (Figure 3A: green and Figure 3B), tea samples of Fengqing were dominated and were more distantly related to the ancestor (inner node), indicating that group 3 represented by Fengqing was more divergent compared to group 1 and group 2. Notably, we found that 15 tea accessions, collected from Fengqing, Other, and Yun of group 2, were clustered with the Cluster 2 of groups 1; Moreover, 4 tea accessions from Yun in subcluster 3 were clustered into Cluster 1 of groups 1. these samples may exist with the possibility of introgression in their clustered regions.

The parameters θπ, H_O, H_E, F, and Tajima’s D were calculated for tea accessions to estimate the patterns of genetic diversity. The primary genetic diversity parameters are presented in Supplementary Tables 5–9, Table 4 and Figure 4. For the subgroup, the expected heterozygosity (H_E) of the tea populations varied between 9.04 (group 2) and 10.25% (group 3), while the observed heterozygosity (H_O) of the tea populations ranged between 3.05% (group 2) and 3.43% (group 3). The inbreeding coefficient (F) of the tea populations varied between 65.52% (group 1) and 69.66% (group 3). It is worth noting that a degree of variability existed between subgroups in both H_E and F (Figure 4A). Nucleotide diversities (θπ) in the three subgroup were estimated at the individual level after the correction for sample size. The analysis found that group 1 (8.89 × 10⁻⁴) had the highest π value, while group 3 (3.31 × 10⁻⁴) had the lowest π value (Table 4 and Figure 4C). Meanwhile, tajamaD analysis showed positive Tajima’s D test values for all subgroup in Lincang (group 3 < group 2 < group 1), in agreement with the findings of π analysis, revealed that the Lincang tea populations may be experiencing group constriction, which may be related to the directional selection. For the regional population(Supplementary Figure 6; Supplementary Table 9), Cangyuan population had the highest nucleotide diversity (1.777 × 10⁻³), which is consistent with the expected heterozygosity and the observed heterozygosity. Otherwise, Tajama’D of Cangyuan is closest to 0 concerning the other eight populations, indicating that Cangyuan is the least selected in Lincang and preserves a large amount of tea germplasm, which is potential resources to expand the genetic resources of improvement. The lowest level of diversity was found for the Shuangjiang population (1.071 × 10⁻³), which resulted from a long history of breed formation and selective breeding more than in most other areas.

The pairwise F_ST between individual populations is presented in Figure 3C. Our analysis revealed that group 1 and 2 (0.092) are more closely distant genetically, while group 1 and 3 (0.226) are more distant from each other. It is inferred that the ancestral population of subpopulation 3 may have diverged prior to the ancestral populations of subpopulations 1 and 2, resulting in greater genetic divergence. ln addition, F_ST values between regions indicate a weak genetic differentiation between regional populations (Supplementary Figure 7).

In the Mantel tests, we can reject the null hypothesis that these three matrices, genetic distance and environment distance (r₁ = 0.1462), genetic distance, and geographic distance (r₂ = 0.07843), are moderately related with alpha = 0.01 (all p ≤ 0.0001). The observed correlation (r₁ = 0.1462, r₂ = 0.07843) suggests that the matrix entries are positively associated. This means that the larger geographical distances between tea accessions lead to greater genetic differences, and the higher differences in environmental variables result in more genetic distance.

In the partial Mantel test, the genetic and environmental distances still were moderate correlated (r₃ = 0.09541, p = 0.001) after considering geographic distance. This result suggests that environmental selection has a weak but non-negligible effect in shaping genetic variation in wild tea germplasm after individual control distance isolation. It can be inferred that there may be introductions of the tea accessions from each other resulting in the weak regional characteristics of tea trees in these populations, which is corresponding to the results of regional phylogenetic analysis (Supplementary Figure 5).

Linkage (LD) is a non-random combination of alleles at different positions in a given population, and usually expressed as D and r₂ values (Slatkin, 2008), which are mainly related to whether the same species have experienced domestication pressure, regional selection pressure, and nucleotide diversity. The LD analysis showed that the r₂ values of accessions of all subgroups are <0.4 (Figure 4C). When the r₂ value is <0.4, it is generally considered that there is no effective linkage or no linkage. At the same time, Linkage disequilibrium decay distance indicated that the degree of Linkage in Lincang is pretty low, which is consistent with the information that the tea samples is basically ancient tea trees with a short history of artificial cultivation.

Population history includes events such as population bottlenecks, expansions, migrations and admixtures, which have important implications for the formation of genetic polymorphic patterns in populations. To reveal the evolutionary history of Lincang tea, we applied the multiple sequentially Markovian coalescent (MSMC; Schiffels and Durbin, 2014) model to the analysis of phased SNP data from three subgroups. As displayed in Figure 5H, similar evolutionary patterns were observed for the three subgroups of Lincang. Group 1 (green line) and group3 (red line) manifested a slight Ne expansion around 10–1,000 Kya, 2–100 Kya and a subsequent Ne contraction (Ne ≈ 80,000 down to Ne ≈ 1,800 and Ne ≈ 100,000 down to Ne ≈ 2,000) around 0.5–10 Kya, 0.4–3 Kya, respectively. Interestingly, the range of expansion and contraction in group 2 is significantly greater compared to groups 1 and 2, inferring that group 3 is likely to be more heavily influenced by human and environmental influences.

Also, TreeMix was used for inferring historical segregation and admixture of populations based on genome-wide SNP allele frequency data. Regarding the changes in tea gene flow, Figures 5E,F both show that there is a Treemix vector connected Fengqing to Zhenkang, and the proportion of gene exchange is relatively large, 38.43 and 36.57%, respectively, indicating bidirectional gene flow. This indicates that hybridization likely occurs among the different geographic areas in nature.

To study the selective characterization of three subgroups of Lincang, 611, 1,328, and 2,090 selective sweep regions were both detected in different tea plant subpopulations, which harbored 65,145 and 201 candidate genes, respectively (Supplementary Tables 10–13, Figure 6, and Supplementary Figures 9–11). In group 1, the highest signal 88.62 was found on chromosome 7 at position 26,812,875 bp, and the most selective signals were detected on chromosome 7 (278 signals). In group 2, 344 selective signals were identified on chromosome 6, and the highest signal reached to 187.6 on chromosome 12 at 96,426,312 bp. For group 3, chromosome 15 (354 selective regions) has the most selective regions, and chromosome 12 at 96,038,055 bp has the most selective signal 247.9. GO enrichment of the 65,145 and 201 domesticated candidate genes of three subgroups showed significant functional representation in the GO categories of negative regulation of cellular process in group 1; peptidase S8/S53 domains, oxidoreductase activity, Golgi apparatus, obsolete oxidation–reduction process, cellular response to lipid in group2; and UDP-glucosyltransferase activity and glucosyltransferase activity in group3 (Supplymentary Figure 12). Interestingly, The peptidase S8/ S53 domains enriched in subpopulation 2 are functionally consistent with those enriched in previous studies of small-leaf tea (Wang et al., 2020), and it is inferred that subpopulation 2 may contain domesticated and hybrid species of small-leaf tea. On the drawback, the lack of functional clarity of many high-signal positional candidate genes makes it impossible to identify selective differential traits between different subpopulations.