Ancient polymorphisms and divergence hitchhiking contribute to genomic islands of divergence within a poplar species complex.

How genome divergence eventually leads to speciation is a topic of prime evolutionary interest. Genomic islands of elevated divergence are frequently reported between diverging lineages, and their size is expected to increase with time and gene flow under the speciation-with-gene-flow model. However, such islands can also result from divergent sorting of ancient polymorphisms, recent ecological selection regardless of gene flow, and/or recurrent background selection and selective sweeps in low-recombination regions. It is challenging to disentangle these nonexclusive alternatives, but here we attempt to do this in an analysis of what drove genomic divergence between four lineages comprising a species complex of desert poplar trees. Within this complex we found that two morphologically delimited species, Populus euphratica and Populus pruinosa, were paraphyletic while the four lineages exhibited contrasting levels of gene flow and divergence times, providing a good system for testing hypotheses on the origin of divergence islands. We show that the size and number of genomic islands that distinguish lineages are not associated with either rate of recent gene flow or time of divergence. Instead, they are most likely derived from divergent sorting of ancient polymorphisms and divergence hitchhiking. We found that highly diverged genes under lineage-specific selection and putatively involved in ecological and morphological divergence occur both within and outside these islands. Our results highlight the need to incorporate demography, absolute divergence measurement, and gene flow rate to explain the formation of genomic islands and to identify potential genomic regions involved in speciation.

Understanding how and why genomes diverge during speciation is fundamental to an understanding of how species evolve (1). With the advance of high-throughput sequencing technologies, considerable progress has been made in documenting the genomic landscape of divergence between recently evolved species (2–23). Most genomic speciation studies have revealed a highly heterogeneous pattern of genomic divergence, with peaks of elevated differentiation scattered across the genome. This “mosaic” pattern of elevated genomic divergence has been hypothesized to be associated with the presence of genomic regions that are resistant to gene flow (1, 23), most likely because they contain barrier variants related to local adaptation and/or reproductive isolation (24). Under this speciation-with-gene-flow model, the size and number of these regions, generally called genomic islands (1, 2, 25), are predicted to expand due to divergence hitchhiking. Moreover, both are expected to grow with divergence time and greater gene flow between incipient lineages (26, 27).

However, genomic islands with elevated divergence (usually quantified by Wright’s fixation index, FST) may also occur in the absence of recent gene flow and can be created by sorting of ancient divergent haplotypes, recent selection at ecologically relevant loci regardless of gene flow and recombination rate, and/or recurrent background selection and selective sweeps in regions of low recombination (12, 21, 28–31). These alternative deterministic and stochastic processes are not mutually exclusive and are therefore difficult to disentangle. Dxy, an absolute measure of genetic divergence between incipient lineages, is less affected than FST by genomic reductions in genetic diversity (12). Thus, Dxy is expected to be elevated in regions with reduced gene flow but unchanged or decreased in regions under the effect of recurrent background selection or selective sweeps (21). This is because recurrent background selection or selective sweeps tend to reduce genetic diversity in the ancestral population, which leads to reduced Dxy due to decreased time to the most recent common ancestor between individuals representing the two descendent lineages (1, 23). Such a scenario should be more apparent in genomic regions of low recombination (12, 15, 20, 32), resulting in a general positive correlation between genetic diversity and recombination rate (33). Nonetheless, divergent sorting of ancient divergent haplotypes will also elevate Dxy and cause island-like patterns of genomic divergence to originate between diverging lineages as in the speciation-with-gene-flow model (21, 34). Thus, analyses of Dxy will only discriminate between some of the alternative factors that cause genomic islands of divergence identified by FST (12, 21, 23, 29, 35, 36). These problems can be alleviated by (i) reconstructing the past demography of the lineages concerned and (ii) comparing lineages with contrasted gene flow and times of divergence along the speciation continuum (22, 23). This is especially critical for plants in which interspecific gene flow occurs more frequently than in animals (1).

Populus euphratica Oliv. (Pe) and Populus pruinosa Schrenk (Pp) are two sister species in section Turanga (37–39), both of which are well adapted to the extreme desert environment (). The natural range of P. euphratica extends discontinuously from Morocco in North Africa through the Middle East to central and western Asia, whereas P. pruinosa is restricted to western China and adjacent regions (40) (Fig. 1). Despite their close relationship, these two species are morphologically well differentiated and grow in distinctly different desert habitats. Leaves of P. euphratica are glabrous and vary in shape from linear or lanceolate on juvenile plants and basal twigs of mature plants, to rounded, ovate, or obovate on older twigs of mature plants. In contrast, leaves of P. pruinosa are always ovate or kidney-shaped with thick hairs (). Although their flowering periods partially overlap, P. euphratica flowers earlier than P. pruinosa in western China (41). It is also evident that where they co-occur, P. euphratica occupies sites with deep underground water, whereas P. pruinosa occurs where underground water is close to the surface, near ancient or extant rivers. Both species form clones by vegetative propagation, mainly through root suckers (38, 40). Since asexual reproduction leads to more extensive linkage disequilibrium, it can play a part in speciation and the maintenance of species. A previous study of nuclear sequence variation in samples from western China has shown that these two species diverged during the Pleistocene in the presence of strong gene flow (42). Therefore, P. euphratica and P. pruinosa provide an excellent system for evaluating how various evolutionary forces could have shaped patterns of genomic divergence during speciation, especially for tree species.

Fig. 1.

Phylogenetic and population genetic analyses of P. euphratica and P. pruinosa. (A) The biogeographic regions of whole-genome–sequenced individuals for P. euphratica (three lineages identified here: PeC, PeE, and PeW) and P. pruinosa (a single lineage identified here: Pp). The colored areas indicate biogeographic distribution of P. euphratica (gray) and P. pruinosa (pink), and the frequencies of each lineage in potential hybrid and/or backcross populations are shown using larger pie charts. (B) Phylogenetic network inferred using the Neighbor-Net method based on genome-wide SNPs. The ancestral state was identified by genotyping the segregating sites with 10 additional poplar species (see for detailed information). (C and D) PCA plots of SNP data for P. euphratica (+) and P. pruinosa (x). (E) Population structure bar plots. Each vertical bar represents a single individual, and the height of each color represents the probability of assignment to that cluster. The number of individuals in each lineage is also shown.

We resequenced the genomes of both species across their biogeographical distributions to examine genomic patterns of divergence, and in particular to infer the relative roles played by gene flow, divergence time, lineage sorting, recombination, and selection during the speciation process. We resolved a paraphyletic pattern of divergence between the two species with four lineages originating during two continuous phases of divergence. Furthermore, we found that divergent sorting of ancient polymorphisms and divergence hitchhiking contributed mainly to the heterogeneous patterns of genomic divergence between the lineages detected. Levels of gene flow and time since divergence were found not to affect the general formation, size, and number of genomic islands that distinguish lineages. Our study not only provides insights into the speciation histories and the corresponding genomic changes of the two desert poplars but also identifies important lineage-specific genes under positive selection that could be responsible for the ecological, physiological, and morphological diversification between these two species.

Results

Genome Resequencing and Genetic Relatedness.

We resequenced the genomes of 122 P. euphratica and 34 P. pruinosa individuals spanning their worldwide geographic distributions (mean sequencing coverage of 11.36×; Fig. 1 and ). Sequences were first aligned to the reference genome of P. euphratica (43), and population genetic statistics based on the site frequency spectrum were then estimated using methods that accounted for uncertainty in the assignment of genotypes (44, 45). Additional tools (46, 47) were used to increase the stringency of variant calls when required. In total, 4,285,372 and 2,460,093 high-quality SNPs were identified in P. euphratica and P. pruinosa, respectively, of which 1,688,107 were shared between the two species (39.4% and 68.6% of the total SNPs in each species, respectively).

Given that poplar trees are known to form clones by vegetative propagation, mainly through root suckers (38, 40), we estimated the genetic relatedness between individuals on the basis of pairwise comparisons of SNP data. We found extensive clonal growth in almost all populations, with P. euphratica populations in Spain, Algeria, and Israel most extreme in this regard (). The latter likely reflects strong founder effects by a few or even a single individual, followed by extensive clonal propagation. Therefore, we discarded individuals that were more related than third-degree relationships, retaining a total of 99 individuals () for subsequent population genetic analyses.

Population Structure.

To examine genetic relationships between individuals, we constructed a phylogenetic network based on nuclear genome polymorphisms. This revealed an evolutionary history distinctly different from the phenotype-based classification (Fig. 1). P. euphratica comprises three lineages (PeE, PeW, and PeC) and is paraphyletic with P. pruinosa; that is, P. pruinosa individuals belong to a single lineage (Pp) sister to one of the three P. euphratica lineages, PeW. The other two lineages, PeC and PeE, are closely related. PeW is mainly distributed in western Iran and extends to Israel, while PeE occurs mainly in eastern Iran. Both PeC and Pp mainly occur in western China and adjacent regions (Fig. 1). This paraphyletic pattern was also confirmed by our maximum-likelihood phylogenetic analyses based on complete chloroplast sequences (), although, in this case, a few P. pruinosa individuals grouped with the lineage PeC, which may reflect recent gene flow between P. pruinosa and PeC in western China. The results of a principal components analysis (PCA) of SNP genotypes across all individuals clearly reflected a history of lineage divergence (Fig. 1 ). The first principal component (PC1; variance explained = 25.70%; Tracy–Widom test, P = 1.99e-14) separated PeC and PeE from the two other lineages, while the second one (PC2; variance explained = 12.73%; Tracy–Widom test, P = 7.94e-44) separated PeW from Pp. Lineages PeC and PeE were clearly differentiated by PC3 (variance explained = 3.52%; Tracy–Widom test, P = 1.87e-134), confirming them as genetically distinct. The PCA also showed some individuals with intermediate positions among PeW, PeE, and PeC, indicating population admixture and recent hybridization and backcrossing between these lineages. Population structure analyses using NGSadmix (48) (Fig. 1) further indicated that potential hybrids and backcrosses (with >10% of the genome from another lineage) between PeW and PeE are distributed in central Iran, while those between PeC and PeE occur in western China (Fig. 1). These admixed individuals were excluded from downstream analyses (). Moreover, this analysis showed that individuals from Spain and Algeria are hybrids between PeW and Pp. Because such hybrids occur disjunctively from the current distributions of both PeW and Pp lineages, it is feasible they were introduced as clones to these areas by humans. Alternatively, both lineages may have been much more widely distributed in the past, with the hybrids in Spain and Algeria representing remnants of their previous co-occurrence and contact in these areas.

The genetic diversity (π) of each lineage ranged from 5.07 to 8.30 × 10−3 (). Lineage Pp exhibited the highest genetic diversity and showed much more extensive genome-wide linkage disequilibrium (LD) than PeC and PeW (). This could be partly explained by a more recent bottleneck in Pp and strong asymmetric genetic introgression from PeC into Pp as previously suggested (42). Pairwise genome-wide averages of genetic divergence (FST) between lineages ranged from 0.19 to 0.36, indicating increased divergence times among lineages. This pattern was also evident from the mean pairwise nucleotide difference in interlineage comparisons (Dxy) ().

Demographic History.

To investigate the evolutionary history of the four lineages, we first examined the historical changes in effective population size (Ne) for each lineage using the Pairwise Sequentially Markovian Coalescent method (49). This showed that all populations/lineages experienced a period of dramatic decline in Ne starting ∼1 Mya (Fig. 2 and ). Algerian and Spanish populations of P. euphratica exhibited the smallest historical Ne, consistent with a clonal origin. The lineage PeW underwent a more severe bottleneck than lineages PeE and Pp ∼50–300 kya, while PeC maintained a relatively large Ne between 100 and 300 kya followed by a sharp decrease. To evaluate alternative divergence models (), we further analyzed pairwise joint site frequency spectra using a composite likelihood approach as implemented in fastsimcoal2 (50). The best-supported model (Fig. 2 and ) suggests that the common ancestor of PeC and PeE split from the common ancestor of PeW and Pp at least 1.86 Mya [95% highest posterior density (HPD) = 1.73–1.90 Mya], while PeW and Pp diverged 0.88 Mya (95% HPD = 0.80–0.94 Mya) and PeC and PeE diverged more recently, ∼80 kya (95% HPD = 75–84 kya). In addition, our simulations suggest that the two continued phases of divergence gave rise to the four lineages and these were accompanied by different rates of gene flow. Gene flow was highest from PeC to Pp (>PeE to PeC > PeE to PeW > Pp to PeC), as also suggested by TreeMix analyses (51) (Fig. 2 and ). Moreover, the much more extensively shared identical-by-descent haplotypes of PeC and Pp indicate that gene flow occurred recently between these two lineages (Fig. 2). Thus, the contrasting patterns and rates of gene flow between lineages provide a good model for examining the potential effect of gene flow on patterns of genomic divergence and speciation.

Fig. 2.

Inferred demographic history for P. euphratica (PeC, PeE, and PeW) and P. pruinosa (Pp). (A) Changes in effective population size (N) through time inferred by the Pairwise Sequentially Markovian Coalescent model. (B) Schematic of demographic scenario modeled using fastsimcoal2. Split times (kya), population size, and migration rates correspond to 95% CIs obtained from this model are shown in . Estimates of gene flow between populations are given in the migration fraction per generation. (C) The maximum likelihood tree inferred using TreeMix with two allowed migration events. P. ilicifolia was assigned as an outgroup. The scale bar shows 10 times the average SE of the entries in the sample covariance matrix, and migration edges are depicted as arrows colored by migration weight. (D) Estimated haplotype sharing between individuals. Heat-map colors represent the total length of IBD blocks for each pairwise comparison. IBD, identical-by-descent.

Genomic Islands of Divergence Between Lineages.

The shape of the FST distribution varied between population pairs (Kolmogorov–Smirnov test, P < 2.2e-16 for all pairwise comparisons) (). Consistent with a divergence scenario that includes a burst of recent gene flow, the FST values between PeC and Pp showed an l-shaped distribution with only 9,546 fixed differences (). This contrasts with comparisons between PeC and PeW, which exhibit similar divergence times but reduced introgression with a greater proportion of highly divergent loci (75,467 fixed differences, accounting for 1.99% of the total number of polymorphisms in these two lineages). We next surveyed the genomic landscape of divergence using a sliding window approach (Fig. 3 and ). As previously found in other species (9, 15, 16), genetic divergence along the genome was highly heterogeneous, irrespective of the level of genome-wide differentiation and geographical separation. In addition, coalescent simulations using msms (52) under the best-fitting demographic model showed that the distributions of simulated polymorphisms and divergence were generally in agreement with the observed patterns across the genome (). This indicates that much of the genomic heterogeneity in divergence can simply be attributed to the effects of demographic processes and genetic drift (31, 53).

Fig. 3.

The heterogeneity of genomic divergence in P. euphratica (PeC, PeE, and PeW) and P. pruinosa (Pp). (A) Pairwise genetic divergence (FST) in 10-kb sliding windows along an example chromosome (chromosome 1) for all comparisons (see , for all chromosomes). Genomic islands of divergence are shown in red. (B) Correlation of 10-kb window-based estimates of genetic divergence (FST), absolute divergence (Dxy), nucleotide diversity (π), and population-scaled recombination rate (ρ) for all six pairwise comparisons. Red circles indicate a positive correlation, while blue ones indicate a negative correlation. The color intensity and circle size are proportional to Spearman’s correlation coefficient.

To further understand the evolutionary forces shaping the heterogeneous landscape of genomic divergence, we detected outlier windows for each pair of lineages. These were defined as those falling within the top 1% of the empirical FST distribution that exhibited significant (false discovery rate < 0.01) divergence compared with a null distribution obtained using a permutation approach (16). Across all comparisons, a total of 1,869 extreme 10-kb outlier windows were identified (). A few of these outlier windows were shared between pairwise comparisons (). After examining their genomic distribution and overlap relationships, we combined adjacent outlier windows and established that these genomic islands of divergence were rather small and scattered across the genome (Fig. 3 and ). In addition, these islands were characterized by strongly reduced levels of nucleotide diversity (π), skewed allele frequency spectra toward an excess of low-frequency variants (i.e., more negative Tajima’s D), elevated levels of LD, and significantly higher proportions of positively selected genes [PSGs; 108 PSGs of 1,537 genes in these regions; χ2 test, χ2 (1) = 210.62, P < 2.2e-16], indicating that linked selection acted on these regions ().

We further examined what factors contributed to the formation of genomic islands based on values of Dxy (12, 21). We compared the level of Dxy in genomic islands to that of genome background in each pair of lineages. We recovered significantly elevated Dxy in genomic islands for all pairwise comparisons (), regardless of rate of gene flow between lineages under comparison. This observation suggests that haplotypes that include genomic islands may have become genetically isolated before the rest of the genomes of the lineage pairs under comparison. No significant difference was observed in the number and size of genomic islands among comparisons for pairwise lineages with contrasting amounts of gene flow and divergence times [number of genomic islands: χ2 test, χ2 (5) = 5.35, P = 0.375; size of genomic islands: nonparametric Kruskal–Wallis rank sum test, χ2 (5) = 5.53, P = 0.355]. We additionally examined changes of recombination rates and found that the population-scaled recombination rate (ρ = 4Nec) was significantly reduced in genomic islands compared with the rest of the genome for all paired comparisons (). However, the correlation between recombination rates ρ and absolute divergence Dxy across the whole genome, although significant, was very low (close to zero; Fig. 3 and ).

Genes Under Positive Selection.

Despite being closely related, P. euphratica and P. pruinosa differ in ecology, physiology, and morphology (41, 42). To explore the genetic basis of this differentiation, we conducted HKA (Hudson-Kreitman-Aguadé) (54) and PBS (Population Branch Statistic) (55) tests to identify genes with high interlineage divergence under recent natural selection following their paraphyletic divergence. Among 28,148 genes analyzed (), a total of 451 were identified under positive selection (188, 204, and 74 in lineages PeC, PeW and Pp, respectively; ). In addition, 15 of these PSGs were shared between any two lineages, a proportion that is significantly higher than expected by chance (for simulation, we randomly sampled 188, 204, and 74 out of 28,148 genes without replacement for lineages PeC, PeW, and Pp, respectively, and a maximum of 12 overlaps between any two lineages were observed after 100,000 replications). This suggests the presence of diverged haplotypes in their ancestral population possibly followed by strong divergent selection. Importantly, only 108 PSGs were located inside identified genomic islands, while the rest were distributed outside them. Gene ontology enrichment analyses of all PSGs yielded several candidate genes associated with organ development, local adaptation, and reproductive isolation (). Among these genes, four (GIS, MYB66, CPC, and RSL2) are associated with trichome differentiation and branching, whereas eight (TCP13, LCR, RDR6, CRF5, YUC4, TRN2, NAC083, and IAA24) are involved in regulation of leaf development and morphology (). GIS encodes a putative C2H2 transcription factor that plays a key role in trichome initiation and branching in Arabidopsis thaliana. In addition, GIS loss-of-function mutations exhibit a premature decrease in trichome production on successive leaves, stem internodes, and branches, while overexpression has an opposite effect and causes significant delay in flowering compared with wild-type (56). TCP13 is a member of the CINCINNA-like TCP (TEOSINTE BRANCHED1, CYCLOIDEA, PCF) group of transcription factors, which play a pivotal role in signaling pathways that generate different leaf forms via negative regulation of the expression of boundary-specific genes (57, 58). LCR (LEAF CURLING RESPONSIVENESS) encodes an F-box protein (SKP1-Cullin/CDC53-F-box) involved in the regulation of leaf curling-related morphology. Overexpression of LCR leads not only to abnormal leaf development but also to a salt-tolerant phenotype of A. thaliana in an abscisic acid-dependent manner (59, 60). Moreover, we found that several genes subjected to selection encoded proteins associated with reproductive function, including flowering time (e.g., EFM, HLP1, CIB1, and GAS41), anther development (RPK2), pollen tube growth and guidance (GLR5, PIP5K4, and DOK1), and the development of male and female gametophytes (ARID1, PIN8, DHAD, and RH36) (). These reproductive proteins have diverged rapidly across lineages and emerged as candidates involved in postzygotic isolation between lineages.

Discussion

Our population genomic analyses of P. euphratica within its vast but fragmented range across Eurasia and North Africa, and additional analyses of P. euphratica and P. pruinosa within the sympatric area in west China, together suggest a scenario of paraphyletic divergence of these two morphologically delimited species. Four lineages were identified: PeW, PeE, and PeC representing P. euphratica, and Pp representing P. pruinosa (Fig. 1). Coalescent-based simulations indicate that these lineages diverged during the Pleistocene. Although the three P. euphratica lineages show limited morphological differentiation, the possibility of incipient speciation occurring within this species is supported by the finding that divergence between PeW and P. pruinosa occurred later than between their ancestral clade and the ancestral clade of the other two lineages, PeE and PeC. While P. euphratica is genetically divided into three lineages, there is little congruence between the current geographical distributions and the genetic distances between them, suggesting a complex expansion and colonization history of the total complex, with some populations possibly established or affected by humans (e.g., Spanish populations and those from North Africa and Central Asia).

Few speciation studies have explicitly reconstructed the demographic history and estimated the amount of gene flow between diverging lineages (23, 24). In this study, we assessed the timing and amount of gene flow and changes in population size of the recovered lineages, revealing pairs of lineages with either large or limited amounts of gene flow occurring between them. In particular, a relatively high level of recent gene flow was shown to have occurred between PeC and Pp. The level of divergence between pairs of lineages was highly heterogeneous along the genome, which could be largely explained by historical demographic processes and genetic drift as suggested by our simulations (). In addition, we identified numerous small genomic islands (quantified by FST) scattered across the genome. Dxy was elevated in these islands (), which is consistent with a model in which these islands were derived from divergent sorting of ancient divergent haplotypes (21). An alternative explanation is the speciation-with-gene-flow model in which restricted gene flow at genomic islands would be expected to increase both absolute and relative measures of genomic divergence (12). Under this alternative hypothesis, genomic islands are predicted to be more pronounced with increased gene flow and divergence time (26, 27). However, our findings indicate that Dxy was of the same magnitude across genomic islands, and no significant difference in size and number of genomic islands was observed among interlineage comparisons with contrasting levels of gene flow and divergence times. Therefore, our observations suggest that recent gene flow is unlikely to be a major factor in generating genomic islands of divergence in these desert poplars.

Several selection signals, including elevated levels of LD and reduced levels of genetic diversity and Tajima’s D, were observed across genomic islands. In addition, we observed a strong association between islands of divergence and low recombination rates. Such an association suggests the possible occurrence of divergence hitchhiking, or background selection and/or recurrent selective sweeps within regions of low recombination (12, 21), both of which tend to reduce genetic diversity and increase FST. If genomic islands are mainly produced by low recombination rate and linked selection, levels of Dxy within them will not be elevated (12, 21), which is inconsistent with our observations (). In conjunction with the covariation of population genomic parameters (Fig. 3 and ), genomic islands of divergence observed here are best explained by divergent sorting of ancient polymorphisms predating lineage divergence (12, 15, 20, 21, 61). The link between low recombination rates and high Dxy of genomic islands most likely resulted from “divergence hitchhiking” because divergent sorting of anciently diverged haplotypes tends to reduce gene flow in surrounding linked regions (62).

Interestingly, we found that genes under positive selection and exhibiting high interlineage divergence were detected both within and outside genomic islands (). These genes are putatively functionally related to organ development, local adaptation, and reproductive isolation, which may underlie ecological and morphological divergence between lineages (). Several genes showed high divergence and selection signatures between P. pruinosa and the P. euphratica lineages. Two of these, TCP13 and IAA24, are involved in the regulation of leaf development and morphology, of which a major difference exists between P. euphratica and P. pruinosa. Similarly, the fixation of different alleles for the gene GIS could be a cause of dense trichome development on leaves of P. pruinosa, while high allelic divergence and fixation of the gene GAS41 in P. pruinosa may contribute to its later flowering relative to PeC (40, 41). In addition, numerous single-nucleotide differences outside these islands or undetectable here may have also contributed to dramatic differentiation between lineages (or species) through an accumulation of small individual effects underlying polygenic barriers (63). What drove cryptic differentiation within P. euphratica and between it and P. pruinosa remains to be investigated in the future.

In conclusion, our results revealed lineage divergence within P. euphratica and a paraphyletic relationship of this species with P. pruinosa. Genomic islands of divergence that distinguish the four lineages of this poplar complex most likely derived from divergent sorting of ancient polymorphisms and divergence hitchhiking but not heterogeneous gene flow. Our findings highlight the need to integrate information on demographic history, gene flow, and absolute genomic divergence to distinguish evolutionary processes that generate genomic islands and identify genomic regions potentially involved in speciation.

Materials and Methods

Silica gel-dried leaves were collected from western China, Kyrgyzstan, Iran, Israel, Spain, and Algeria (). Total genomic DNA was extracted using a standard cetyltrimethylammonium bromide protocol and sequenced on the Illumina HiSeq 2000 and 2500 platforms. After data-quality control, reads were mapped to the reference genome using BWA-MEM. Population genetic structure was conducted using the Neighbor-net method, PCA, and admixture estimation. Population genetic statistics π, FST, Dxy, and Tajima’s D were estimated based on the site frequency spectrum inferred by analyses of next generation sequencing data. Detailed information on materials and methods with associated references are available in the .