1,135 Genomes Reveal the Global Pattern of Polymorphism in Arabidopsis thaliana

1Arabidopsis thaliana remains at the forefront of modern genetics. Decades of work have not only established much of what we know about the physiology and development of plants but also provided insight into how wild populations adapt to biotic and abiotic environments. Few systems share the key advantage of A. thaliana for GWAS or complementary forward genetics approaches: the ready availability of a large collection of naturally inbred lines (accessions) that are products of natural selection under diverse ecological conditions. This makes it possible to link genotypes and phenotypes to fitness effects in the laboratory and the field. By adding molecular data for genetically identical individuals—e.g., RNA expression or epigenetic marks— underlying mechanisms can be elucidated much more easily than in other multicellular species.

2The fundamental resource for this research program is a set of accessions with complete genome sequences, collected from different locales. Systematic characterization of genome-wide polymorphism in A. thaliana, paralleling efforts in humans, began with a description of linkage disequilibrium and population structure in 96 accessions. This was followed by a whole-genome map of deletions and SNPs in 20 global accessions, which in turn was the basis of a 250k SNP array with multiple markers in each haplotype block. This array was used to genotype the RegMap collection of 1,307 diverse accessions. Concurrent with the application of short read sequencing to human genomes, the first A. thaliana genomes were resequenced, soon followed by the analysis of larger collections. Similar efforts have led to large-scale surveys of sequence diversity in other plants, mostly crops.

3We extend these efforts with 1,135 A. thaliana accessions from a worldwide hierarchical collection. There were several motivations for the current study: to quantify genome variation in a larger and more representative sample of accessions; to investigate the demographic history of the species; to identify features that make specific geographic or genetic subsets particularly well suited for forward genetics, field experiments and selection scans; and to provide a powerful GWAS platform. Previous studies had shown that the ability to detect footprints of selection depended greatly on the sample (e.g.,,, ). Similarly, while GWAS have identified common alleles with major effects from as few as 96 accessions, a much larger sample is required for most traits. The SNP-genotyped RegMap panel provided such a collection but did not efficiently capture all SNPs and structural variants. Fully sequencing this collection would be of limited benefit, as one could accurately impute the missing data by sequencing a subset. We therefore assembled a set of accessions that sufficiently overlap the RegMap panel for imputation of variants in all lines. The combined collection constitutes a first-rate resource for determining how genetic variation translates into phenotypic variation.

4We selected accessions for Illumina short read sequencing with several objectives in mind. We sought to cover the global distribution of A. thaliana more evenly than the RegMap panel while including large regional collections of particular interest from ecological and evolutionary perspectives, notably from Sweden and the Iberian Peninsula (Figure 1A). We also wanted to better sample interesting regions based on prior knowledge of the population structure of the species, such as North America and Central Asia. Our collection is hierarchical, with a range of geographic distances between nearest neighbors, and a few very densely sampled locales. Most accessions had been genotyped with 149 genome-wide intermediate-frequency SNP markers to avoid sequencing identical individuals.

5After filtering (described below), we retained sequences of 413 RegMap and 722 new lines, for a total of 1,135 accessions with whole-genome information (see the Data Release section). These 1,135 lines are the focus of this paper; the imputed RegMap set will be described in another paper. Together, the RegMap and 1001 Genomes samples include 2,029 natural A. thaliana accessions with high-quality polymorphism data (Figure 1B).

6The genomes presented here integrate previously published subsets (,,,,,; Figure 1B). All accessions are available from the stock centers, and we have generated an accession list (see Data Release section) that unifies previous naming schemes and provides provenance information. Our intention is for this collection to remain actively curated as ever more accurate genomes are produced and a wide range of phenotypic data are generated (not only by us, but also by the community— see www.1001genomes.org for information on how to contribute).

7A range of Illumina platforms were used across several sequencing centers and over several years, so we instituted stringent quality controls to pare an initial set of over 1,200 sequenced genomes to a final set of 1,135 (see Data Release section). The data are the intersection of the MPI (SHORE) and GMI (GATK) pipelines, independently validated in our pilot studies. An average of 100 Mb (84%) per line were called against the TAIR10 reference genome (119 Mb). The missing positions differ greatly between lines, such that only 2% of the reference genome lack calls entirely. Based on comparisons with one long read (Pacific Biosciences) and three short read (Illumina) based de novo genome assemblies, we estimate that fewer than 3% of SNP calls are erroneous (i.e., should be reference instead) independently of dataset coverage, with the vast majority being singletons. Over 98% of genotype calls were correct at SNP sites, and only 1.5% of SNPs were mistakenly called as reference (Supplemental Experimental Procedures). We emphasize, however, that this calculation ignores SNPs missed because they are in the vicinity of structural variants, which are difficult to assess with short read technology.

8After filtering, the nuclear genomes contained 10,707,430 biallelic SNPs and 1,424,879 small-scale indels (up to 40 bp). This represents one variant on average every 10 bp of the single copy genome, which is the densest variant map for any organism, including the most recent release of the 1000 Genomes Project for humans. 2,842 biallelic SNPs were called in chloroplast genomes and 824 in mitochondrial genomes. The complete data are available as VCF files and as FASTA pseudogenomes (see Data Release section). We also developed web applications to subset the full genome VCF or pseudogenome files and extract data on a selection of genomes and/or specific loci as well as a “Strain ID” application, with which users can identify the genomes in our sample that are most closely related to a newly sequenced genome (see Data Release section). As with all short read data, we advise caution in using our pseudogenomes with applications in which the contiguity of DNA sequence is critical, for example, in the generation of PCR primers. Finally, we are committed to supporting the community in developing additional applications that make use of these data.

9A major motivation for sequencing a large collection of accessions is to enable GWAS with nearly complete genotype information. For comparison with the RegMap data, we measured flowering time under different environments (10°C and 16°C) in our collection and performed GWAS. We note first that there is little reason not to use full genome data, as permutation-based multiple-comparison thresholds can be used to minimize the statistical cost of additional markers (Figure 2A). The chromosome 5 peak at ∼23.25 Mb nicely illustrates the advantage of the full genome data (Figure 2B). Although clearly visible using the 192,498 biallelic variants from the 250k SNP array, not a single SNP reaches genome-wide significance, and the peak might well have been ignored, were it not for the fact that the most significant SNP lies in an intron of the flowering time regulator VIN3. In contrast, the full data clearly reveal a significant peak. Notably, the lead SNP in this peak is not in linkage disequilibrium with the tag SNP from the 250k SNP array, even though the two variants are only 60 bp apart.

10The remaining peaks (Figure 2A) contain the flowering regulators FT, SVP, FLC, all previously linked to flowering time variation, and the dormancy regulator DOG1, recently shown to affect also flowering time. As previously noted, linkage disequilibrium is normally too extensive to directly pinpoint the causative genes or variants with GWAS alone. For example, the peak on chromosome 2 that contains SVP also includes At2g22590, which codes for an UDP-glucosyltransferase, a family of proteins linked to the control of FLC expression.

11GWAS provide insights into the genetic basis of natural variation. To interpret such variation, it is essential that we understand the evolutionary history of a species. For an organism such as A. thaliana, the simplest population genetics model is strict isolation by distance (IBD), under which the genetic distance between individuals reflects only geographic distance. This model does not fit, as the peaks of pairwise differences do not reflect geography (Figure 3A). One extreme encompasses groups of (nearly) identical individuals corresponding to inbred lineages, a result of selfing. This includes 78 North American accessions, with several smaller clusters of three to seven members, and 40 pairs of accessions that differ by fewer than 1k SNPs (Figure 3B). 60 additional pairs differ by fewer than 50k SNPs, much less than the median of 439,145 for all comparisons. Excluding North American and British accessions, ∼80% of these nearly identical pairs were collected within 1 km of each other, most within a few meters. It remains unclear whether the remaining pairs represent true long distance migration, rather than mis-assignment or mix-ups after collection (see Supplemental Experimental Procedures).

12Both North America and the British Isles show evidence of recent long-range dispersal. While North America harbors a single lineage due to very recent colonization, the British Isles contains numerous widely spread genotypes, suggestive of a more ancient and gradual colonization. The median geographic distance between nearly identical British pairs is 303 km, and only 1 of 40 nearly identical pairs was collected from the same site. While some pairs may reflect labeling errors after collection, close genetic relationships are also observed among more diverged but still rather similar pairs of British accessions, supporting that they are the product of recent gene flow.

13At the other end, extreme pair-wise divergences (Figure 3A) are seen with 26 accessions, including 22 from the Iberian Peninsula, and one line each from the Cape Verde Islands, Canary Islands, Sicily, and Lebanon (see Data Release section). We refer to these accessions as “relicts.” The 22 Iberian relicts are no more different from each other than are pairs of non-relicts (Figure 3C). The remaining four relicts stand apart from each other and from all other accessions.

14By genetic distance, our 1,135 accessions thus comprise six diverged groups: four relict groups with a single line each; one relict group of 22 Iberian accessions; and the majority group of 1,109 accessions. It should be noted that accessions Mr-0 from Italy and Tnz-1 from Tanzania also were extremely diverged, but their sequences failed quality controls and were not included in the final 1,135 accessions. Re-sequencing confirmed that Mr-0 (closely related to Sicilian relict Etna-2) and two further Tanzanian accessions, Tanz-1 and Tanz-2, are relicts. Their sequences will be available in the next data release.

15The geographic distribution of relicts and non-relicts (Figure 3B) confirms that a naive IBD model cannot hold. For example, Iberian non-relicts are more closely related to accessions from Kazakhstan than to Iberian relicts. Moreover, while relicts show strong IBD on all geographic scales, non-relicts have a similarly clear pattern only over short distances (Figure 3D), as expected if they had spread rapidly to occupy their current, global range, while the relicts had largely stayed put. On a regional scale, there is considerable geographic variation in the strength of IBD among the non-relicts, indicating that the history of colonization is complex (Figure 3E). The existence of outlier accessions, such as Cvi-0 and Mr-0, has been noted before; it is now clear that there are many such accessions, and that they can be common locally. Our data also confirm that the colonization of North America was recent and rapid. In addition to groups of nearly identical individuals, 47% of North American pairs exhibit extensive haplotype sharing (total identity-by-descent length over 85 Mb, as inferred using Beagle and GERMLINE, Figure S1), indicating recent mixing among a limited number of initial immigrants. Conversely, European accessions have low genetic relatedness, and the extent of haplotype sharing generally decays with geographic distance.

16To examine population structure in greater detail, we used ADMIXTURE to cluster the accessions. In addition to identifying most of the relicts as a genetically distinct group, this analysis breaks non-relicts into eight clusters that broadly correspond to geography (see Data Release section). We defined nine groups based on these clusters and assigned each individual to a group if more than 60% of its genome derived from the corresponding cluster. The 135 individuals not matching this criterion were labeled “Admixed.” There is evidence for admixture between the relict and non-relict groups, as two accessions initially identified as relicts, from Sicily and Lebanon, were found to be admixed. These ADMIXTURE classifications were used in all subsequent analyses.

17The ADMIXTURE groups do not correspond to idealized randomly mating populations. There is a regional and variable pattern of IBD (Figure 3E). Similarly, geographic locality prediction using SPA demonstrates the existence of population structure both within and between groups (Figure S2) and highlights the variability in IBD (Figure S3).

18To elucidate the historical processes that have shaped extant diversity, we estimated the distribution of coalescence times for the different populations using MSMC. The results suggest that glacial refugia are largely responsible for present population structure (Figure 4A). Coalescence rates are an indication of relatedness, with higher rates indicating closer average relatedness (and smaller effective population size). Since the last glaciation, coalescence rates within non-relict ADMIXTURE groups were much higher than for Iberian relicts, or between members of different non-relict groups, and coalescence rates between relicts and non-relicts were essentially zero. The rate of coalescence between relicts and non-relicts was also lower than the other rates during the last glaciation, indicating that they were isolated from each other during this period. At the same time, the rate of coalescence among Iberian relicts was high, indicating a local bottleneck, with only slight differences in coalescence rate within and between non-relict groups, consistent with these groups being the product of post-glacial expansion.

19In contrast, current population structure is not reflected in the rate of coalescence before the last glaciation; there has since been sufficient migration and time to erase all traces of earlier population structure. The distribution of highly diverged haplotypes at individual loci in the genome is thus independent of present population structure. The first polymorphism study in A. thaliana, with ADH, noted already the presence of surprisingly diverged haplotypes, and interpreted it as evidence for balancing selection. Many realized that the phenomenon was common and that deep population structure must be responsible. However, the population structure required to account for pairwise sequence divergence of several percent at individual loci, compared to a genome-wide average of 0.5%, must be far older than the most recent glaciation. Thus, while recent coalescence times, as reflected in low pairwise sequence divergence, are more common in within-group comparisons (Figure 4B), the tails of extreme values look very similar for within- and between-group comparisons (Figure 4B, inset).

20It is important not to exaggerate the divergence between relicts and non-relicts. In a survey of four-sample gene genealogies between the Col-0, Ler-0, and two random Iberian relicts, 26% place Col-0 closer to one relict than to Ler-0, rather than the two non-relicts together (Table S1). Two additional four-sample analyses, one with the outgroup A. lyrata, Col-0, and two relicts, and the other with A. lyrata, one relict, Col-0 and another non-relict, produced similar patterns. While genes with the expected topology were significantly more common (43% and 53%, p < 0.001), many genes supported the alternative topology of a non-relict being closest to a relict.

21The last glacial maximum was followed by pronounced expansion of the global A. thaliana population. It is therefore natural to search for footprints of selection related to adaptation to new environments, especially to climate, which varies considerably across the species range. We examined correlations of genetic variants with six climate variables that capture variability in current temperature and precipitation using a mixed model that controls for genome-wide relatedness across samples. Twenty SNPs are significantly correlated with precipitation-related variables, at a False Discovery Rate of 5% (Table S2). Three associations are characterized by a much higher derived allele frequency in the Iberian relicts than in the general population, possibly indicative of local adaptation from new mutations. One affects the ERF1 drought response regulator. ERF1 is also involved in resistance to several fungal pathogens, as is MLO11, which is located near two of the other variants. The connection to drought response and fungal defense suggests that selection could be due to tradeoffs between abiotic and biotic stress.

22Strong local adaptation may create abrupt geographical changes in allele frequency. We used SPA to search for SNPs showing this pattern, and intersected these SNPs with climate GWAS hits to identify variants with local adaptation to climate. Spatial and climatic distributions of genetic variants are intertwined with population structure (Figure S3), so almost no significant variants remain once population structure is taken into account. A single variant associated with precipitation in the wettest quarter also shows a significant geographic gradient (Figure 5A). This variant is in a genomic region densely populated with genes that have been implicated in root development and metabolism, flowering time and flower development, salt tolerance, and detoxification (Figure 5B).

23Because the strategies above attempt to eliminate false positives from population structure, it is difficult to detect variants under population-specific selection. To identify such genes, we calculated FST between admixture groups for all SNPs. The most diverged region is on chromosome 2 at 6.401 Mb, overlapping the gene encoding AGP9, which has not been linked to adaptive processes before (Figures 5C and S4A). Regions adjacent to centromeres exhibit the lowest FST values. In agreement with previous results, these regions contain excessive linkage disequilibrium (ω, Figures 5D and S4B), which suggests that they have been shaped by selective sweeps or background selection.

24In addition to these global patterns, we identified loci that may contribute to adaptive differences between Iberian relicts and non-relicts. We paired each relict with the geographically closest non-relict (Table S3). Over 100 variants have diverged between the two groups, including several in or near EIN2, a development and stress regulator, and AP2, which is important for flower and seed development. Additional genes with differentially fixed variants are LUG and SLK1, which encode transcriptional co-repressors that interact biochemically and genetically with each other and with AP2. Finally, a deeply diverged region around 18.796 Mb on chromosome 2 includes two flowering time regulators, AGL6 and SOC1. As expected from these candidates of selection, the top biological processes (GO terms) strongly overrepresented in these results are “flower development” and “ABA-related activities” (Table S4). Consistent with differentiation in flowering time, relicts flower in 10°C long days on average 21 days later than their nearest non-relicts (t = 4.69, df = 41; p = 3 × 10−5), suggesting that life-history differences contributed to the spread of non-relicts.

25Demographic history can affect the efficacy of selection, and mutations that are likely to be deleterious are common in A. thaliana, especially in marginal populations. We therefore predicted the impacts of coding sequence variants in different genetic groups using SNPeff. Most genes, 27,525, contained at least one variant likely to change protein function, with 17,692 having at least one high-impact variant. On average, 440 genes per accession, for a total of 15,060 genes, had at least one variant predicted to inactivate the gene, although this is likely an overestimate, as it does not account for compensatory mutations or different transcript isoforms. Relicts have proportionally the most genes with potentially deleterious mutations, consistent with a reduced efficiency of selection in the relicts due to small effective population size, with the caveat that mapping to a non-relict reference may again lead us to overestimate such variants (Figure S5).

26The exquisite detail with which we have characterized the spatial pattern of polymorphism in A. thaliana has clarified prior hypotheses and revealed surprising aspects of the species’ history. In particular, the crucial importance of the last ice age has come into much sharper relief. The picture that emerges is that modern A. thaliana is a complex mixture of survivors from multiple glacial refugia, with population expansion having strongly favored the descendants of a specific refugium, possibly as a result of human activity. Under this model, the “relict” accessions are simply those that survived this expansion/invasion. Several lines of evidence support this interpretation.

27The pattern of isolation-by-distance suggests that relict populations have been relatively stationary while the non-relicts’ range rapidly expanded (Figure 3D). Consistent with this model, the climate at the relict locations has changed much less since the last glacial maximum than where modern non-relicts are found (Figures 6A and S6). The Iberian Peninsula is especially interesting, given the presence of a large number of relicts interspersed with non-relicts (Figure 3C). Although relicts are widely distributed there, they are restricted to a very specific environment characterized by old oak and pine forests, high climate seasonality, high temperatures, and low rainfall (Figure 6A). Iberian relicts correspond to a genetic lineage that has been previously identified in the southwestern Mediterranean region, supporting the idea that they survived in a glacial refugium in North Africa. In contrast, non-relicts are found more often in agricultural and urban areas, consistent with expansion of non-relicts having been associated with human activity, and with the relative rarity of relicts reflecting destruction of undisturbed habitats.

28The source of the non-relicts, which comprise most modern A. thaliana individuals, remains obscure. The Iberian Peninsula has the largest regional diversity, and Mediterranean regions tend to be more diverse than other regions (Figure 6B). There is a gradient of decreasing diversity from south to north (Figure 6C), as expected after a range expansion from southern glacial refugia. However, this pattern is likely due to admixture between relicts and invading non-relicts in these regions (high diversity in the Iberian Peninsula almost certainly is) and does not reveal the origin of the invaders. Indeed, omitting relict and admixed accessions, it would be easy to come to the conclusion that the center of diversity is southern Sweden and that diversity decreases from north to south across the entire range. The relatively high diversity seen in southern Sweden (Figure 6D), and also in Russia (Figure 6B) may similarly be the result of admixture, in this case between the original post-glacial colonizers and more recent weedy varieties that accompanied the spread of agriculture, perhaps giving rise to the higher values of Tajima’s D in these regions. Resolving this issue using only contemporary collections will be difficult.

29One pattern that does seem clear is that longitudinal gradients of regional diversity are much weaker than latitudinal ones (Figure 6C), most likely reflecting the relative ease with which organisms in Eurasia can move along the east-west axis. The spread of A. thaliana and other weedy species may have been further enhanced by the rapid expansion of agriculture along this axis. Of particular interest in this respect are unusual populations, such as those in North America, which was colonized only a few centuries ago, and in Central Asia, where reduced genetic differentiation suggests rapid expansion over large geographic distances from a few very small and remote glacial refugia.

30Questions one can address with natural accessions of A. thaliana include how patterns of genetic and epigenetic diversity arose and which forces drive adaptation to the environment. In addition, our knowledge of fundamental molecular processes can be greatly improved through the study of natural change-of-function alleles. Crucial for these purposes is a well characterized, curated, and publicly available collection of accessions. We provide such a collection. Using it as a starting point, increasingly detailed information about (epi)genomes and molecular and non-molecular phenotypes can now be generated. This is a sharp distinction from similar efforts in outcrossing organisms, in which immortalized genotypes are either only available as cell lines, or do not represent adapted genotypes sampled from nature.

31The selection of accessions for ecological field studies and laboratory experiments should take into account their full genetic backgrounds. No subset will be optimal for all purposes. A more diverse sample will contain more genetic heterogeneity, which reduces genetic mapping power but captures more variants. The other extreme is represented by the lineage that has recently colonized North America; while it is phenotypically quite uniform, its low diversity provides an opportunity to study the role of de novo variation in adaptation. One should also take into account the natural history of accessions, including local ecology and climate, which may enable informed decisions about phenotypic variation that is likely to reflect adaptation. For example, temperature and precipitation vary greatly across the species’ range and between groups (Figure S6), and one would expect differences in physiological and developmental responses of Spanish and Swedish accessions.

32Few, if any, systems offer the benefits of A. thaliana: a myriad of sequenced, clonal lineages from a range of ecologically diverse habitats, with patterns of linkage disequilibrium favorable for GWAS, all in an experimentally tractable organism. Another dimension can now be added to traditional functional genomics databases: adaptive variation. The 1001 Genomes collection provides an outstanding opportunity to decipher how genetic variation translates into phenotypic variation and to study the many ways in which plants respond—and have responded—to environmental challenges.

33We initially selected 1,227 worldwide accessions based on genotyping and geographic diversity. They were sequenced by Weigel (MPI), Nordborg (GMI), Ecker (Salk), Mott (Oxford), and Monsanto. Bergelson (University of Chicago) generated the bulk of the seed and tissue used. Paired-end (PE) sequencing employed several generations of the Illumina platform: 1.3+ (80 accessions), 1.5+ (396 accessions), and 1.8+ (751 accessions).

34Variants were called with MPI-SHORE and GMI-GATK (v1.6-5, ) pipelines, validated in our pilot studies. We generated intersection VCF files with high quality in both pipelines. A series of quality checks resulted in a final set of 1,135 accessions, used for further analyses unless mentioned differently. Variant calls were benchmarked using whole-genome alignments of one long read (Pacific Biosciences) and three short read (Illumina) de novo assemblies against the TAIR10 reference. The average true positive rate (TPR) was 98%, the average false negative rate (FNR) 1.5%, the false discovery rate (FDR) 3%, independent of coverage depth used for the variant calls (Table S5). Pseudogenomes were generated by combining reference and variant calls, including indels.

35Please see Supplemental Experimental Procedures for details.

36Data and tools are available at http://1001genomes.org. We uploaded raw reads in FASTQ format for 1,135 final accessions to NCBI SRA (SRP056687). We are releasing the following files at http://1001genomes.org/data/GMI-MPI/releases/v3.1: full VCF variant files for each accession, VCF files with quality reference calls, a combined Full Genome VCF file for all genomes, a standard merged group VCF file without invariant positions, a variant annotated SnpEff VCF file, and individual pseudogenome files. Several tools to facilitate the use of this data are available under http://tools.1001genomes.org, including a strain ID web application, a viewer pf ADMIXTURE group membership, and a tool to retrieve specific regions of pseudogenomes in FASTA. Accession metadata, including group membership, are available under http://1001genomes.org/tables/1001genomes-accessions.html. See supplemental data release for additional tools and datasets.

37The members of The 1001 Genomes Consortium for this project are Carlos Alonso-Blanco, Jorge Andrade, Claude Becker, Felix Bemm, Joy Bergelson, Karsten M. Borgwardt, Jun Cao, Eunyoung Chae, Todd M. Dezwaan, Wei Ding, Joseph R. Ecker, Moises Exposito-Alonso, Ashley Farlow, Joffrey Fitz, Xiangchao Gan, Dominik G. Grimm, Angela M. Hancock, Stefan R. Henz, Svante Holm, Matthew Horton, Mike Jarsulic, Randall A. Kerstetter, Arthur Korte, Pamela Korte, Christa Lanz, Cheng-Ruei Lee, Dazhe Meng, Todd P. Michael, Richard Mott, Ni Wayan Muliyati, Thomas Nägele, Matthias Nagler, Viktoria Nizhynska, Magnus Nordborg, Polina Yu. Novikova, F. Xavier Picó, Alexander Platzer, Fernando A. Rabanal, Alex Rodriguez, Beth A. Rowan, Patrice A. Salomé, Karl J. Schmid, Robert J. Schmitz, Ümit Seren, Felice Gianluca Sperone, Mitchell Sudkamp, Hannes Svardal, Matt M. Tanzer, Donald Todd, Samuel L. Volchenboum, Congmao Wang, George Wang, Xi Wang, Wolfram Weckwerth, Detlef Weigel, Xuefeng Zhou.

38J.B., J.R.E., M.No., M.S., and D.W. coordinated the project. C.A.-B., C.B., J.B., J.C., E.C., T.M.D., J.R.E., A.M.H., S.H., M.H., A.K., P.K., N.W.M., M.Na., T.N., M.No., P.N., F.X.P., B.A.R., K.J.S., F.G.S., M.M.T., D.T., W.W., and D.W. selected and generated the samples. C.B., J.B., J.C., X.G., C.L., B.A.R., M.J., R.A.K., T.P.M., R.M., V.N., R.J.S., F.G.S., M.S., S.L.V., and X.Z. generated and handled sequence data. J.A., F.B., A.F., D.G.G., D.M., P.Y.N., A.P., F.A.R., A.R., C.W., and X.W. performed primary sequence analyses and variant annotation. M.E.-A., W.D., A.F., A.M.H., S.R.H., M.H., A.K., C.-R.L., M.No., H.S., and G.W. performed population genetic analyses. J.F., P.K., A.K., A.P., Ü.S., and C.W. curated the online resources and databases. K.M.B., D.G.G., A.K., M.No., and P.A.S. performed GWAS. M.E.-A., A.F., F.B., A.M.H., M.H., A.K., C.-R.L., M.No., A.P., H.S., C.W., G.W., and D.W. wrote the manuscript.