Selective Loss of Diversity in Doubled-Haploid Lines from European Maize Landraces.

Maize landraces are well adapted to their local environments and present valuable sources of genetic diversity for breeding and conservation. But the maintenance of open-pollinated landraces in ex-situ programs is challenging, as regeneration of seed can often lead to inbreeding depression and the loss of diversity due to genetic drift. Recent reports suggest that the production of doubled-haploid (DH) lines from landraces may serve as a convenient means to preserve genetic diversity in a homozygous form that is immediately useful for modern breeding. The production of doubled-haploid (DH) lines presents an extreme case of inbreeding which results in instantaneous homozygosity genome-wide. Here, we analyzed the effect of DH production on genetic diversity, using genome-wide SNP data from hundreds of individuals of five European landraces and their related DH lines. In contrast to previous findings, we observe a dramatic loss of diversity at both the haplotype level and that of individual SNPs. We identify thousands of SNPs that exhibit allele frequency differences larger than expected under models of neutral genetic drift and document losses of shared haplotypes. We find evidence consistent with selection at functional sites that are potentially involved in the diversity differences between landrace and DH populations. Although we were unable to uncover more details about the mode of selection, we conclude that landrace DH lines may be a valuable tool for the introduction of variation into maize breeding programs but come at the cost of decreased genetic diversity.

Maize is an outcrossing species and has been cultivated for millenia in open-pollinated populations known as landraces. Mass selection in these populations has been highly successful, allowing maize landraces to adapt to a breadth of environments and a wide array of cultural preferences (Bellon ). Over the last century, an inbred-hybrid system has replaced landraces in modern agriculture due to its higher yields and increased stability (Troyer 2001). But the inbred-hybrid system has focused on an ever-decreasing pool of germplasm, restricting genetic variation compared to landraces (van Heerwaarden ).

Though lower-yielding in industrial conditions, landraces continue to serve as an important genetic resource for future crop improvement and adaptation (Sood ; Gates ). But the conservation of landrace diversity imposes a number of challenges. In situ conservation by practicing farmers has been very successful (Bellon ), but is vulnerable to changing economic considerations and does not provide easy access for breeders. Conservation of germplasm ex-situ provides straightforward and safe long-term access to plant breeders, but genebank accessions have to be maintained as large populations to prevent the loss of diversity due to drift (Ellstrand and Elam 1993).

Recently, Melchinger suggested using doubled-haploid (DH) technology as a means of conserving landrace diversity in a homozygous form that would simplify germplasm conservation and be more readily usable by plant breeders (Sood ; Gorjanc ; Melchinger ; Mayer ). The directed induction of DH lines has been developed for several crops to accelerate breeding (Smith ; Gomez-Pando ; Dunwell 2010). The technology permits the instantaneous development of homozygous lines within a single generation instead of six to ten generations of conventional recurrent self-pollination (Prigge ). While not all landraces produce DH lines with equal success, Melchinger concluded that genetically stable DH line libraries of landrace accessions could be used for ex-situ conservation of maize without major loss in genetic diversity.

Theoretical considerations, however, suggest that the instantaneous inbreeding associated with DH production may impact genetic diversity and fitness (Keller and Waller 2002; Charlesworth and Willis 2009). Inbreeding reduces the effective size of a population, resulting in increased genetic drift and a loss of genetic diversity (Charlesworth and Willis 2009). Inbreeding is also thought to impact deleterious alleles and overall genetic load. As an outcrossing species, maize harbors a substantial number of deleterious, partially recessive alleles (Yang ), mostly maintained at low frequencies (Mezmouk and Ross-Ibarra 2014; Yang ) likely at mutation-selection balance (Eyre-Walker and Keightley 2007). Depending on the population history, inbreeding can have opposing effects on deleterious alleles. On one hand, inbreeding exposes recessive deleterious alleles, which can then be purged from the population (Keller and Waller 2002; Henn ). But purging of deleterious alleles is most efficient when inbreeding occurs gradually over several generations (Keller and Waller 2002), rather than instantaneously as expected during DH line production. Instead, instantaneous homozygosity likely decreases the efficacy of selection and increases genetic load. These processes likely contribute to the highly reduced efficiency of DH production in outcrossing maize landraces compared to modern breeding germplasm that has already experienced conventional inbreeding (Böhm ; Melchinger ).

Given these considerations, here we re-evaluate the effects of DH production in maize landraces. We quantify the changes in genetic diversity due to DH production and investigate the role of drift and selection in creating the observed patterns. Combining published genotype and phenotype data from a number of sources (Melchinger ; Mayer ; Brauner ), we analyze and compare samples from five populations of European maize landrace accessions and their derived DH lines. In contrast to previous reports (Melchinger ), we find that landrace genetic diversity is not fully captured by DH line libraries. Although we are unable to pinpoint the causes underlying allele frequency changes at individual outlier loci, we find evidence suggesting that selection against recessive deleterious alleles causes reduced genetic diversity in DH populations. We conclude that DH technology is not suited to conserve maize landraces and its use would result in the loss of potentially important diversity in germplasm collections.

Material and Methods

Data preparation

We used genomic data of five European maize landraces and their derived DH lines to study the effect of instantaneous homozygosity. For the landrace derived DH lines (DH) we used data from Melchinger from a total of individuals genotyped on the Illumina MaizeSNP50 BeadChip (Ganal ). The genotypes were derived from five accessions: Bugard (BU, ) from France, Gelber Badischer (GB, ), Schindelmeiser (SC, ) and Strenzfelder (SF, ) from Germany, and Rheinthaler (RT, ) from Switzerland. For the landrace population samples (LR) of the same accessions, we used data of individuals (, , , , , respectively) from Mayer that were also used in the analyses of Melchinger , genotyped on the 600k Affymetrix Axiom Maize Genotyping Array (Unterseer , Table S1). The LR individuals represent heterozygous subsamples of genetic material maintained as outcrossing populations in genebanks, DH individuals homozygous samples of that source population (Melchinger ).

After combining the two datasets based on physical positions (AGPv2), we removed SNPs that were monomorphic across all accessions in the LR and DH. For all further analyses, we then used updated positions (AGPv4) for the SNPs obtained from supplementary data in Mayer . We also removed SNPs that violated Hardy-Weinberg equilibrium in all LR accessions using exact tests with mid-p adjusted , as well as low quality SNPs that matched the quality criteria ’off-target variant’, ’CallRateBelowThreshold’ and insertion-type SNPs of the Affymetrix Axiom 600k genotyping chip (classifications followed Table S6 in Unterseer ; for details see Table S2 and supplementary note), removing a total of 83011 SNPs. In total, 64930 genotypes (on average 0.6%) in DH individuals remained heterozygous and were set as missing data. Finally, we removed sites that were missing in all individuals by filtering out sites with missing data above 0.99 using plink 1.9 (Chang ). This resulted in 533190 SNPs in the LR, 37967 of which overlapped between the LR and DH. We refer to the smaller set of data hereafter as the ’50k’ dataset.

Unless otherwise specified, all statistical analyses described below were conducted using R (R Core Team 2018).

Phasing and imputation

To project all 533190 SNPs of the LR dataset onto the DH data, we used a two-step approach to phase and impute the data using BEAGLE 5.0 (Browning and Browning 2009) for each LR-DH combination separately. First, we phased and imputed the LR data, then we used this data as a reference to impute the DH lines. Parameters used for BEAGLE were ne = 100000, phase-states = 200, nsteps = 14. To assess the quality of the imputation and establish optimal parameters for the algorithm, we dropped 10000 known SNPs randomly in the DH and calculated imputation error rates for the ith SNP as , with z matches and m missing genotypes out of k individuals. The mean error rate was used to establish optimal imputation and phasing parameters for the algorithm after several runs with different parameters. We compared estimated error rates to diversity and recombination rate to exclude potential imputation biases. Estimated error rates varied from 10.6% in GB to 15.9% in BU, but are not correlated with recombination rate (Figure S1) and appear to be randomly distributed across the genome (Figure S2). Minor alleles are associated with higher error rates, however major alleles are mostly imputed correctly (Figure S3, see supplementary note). We refer to the set with 533190 SNPs with imputed DH individuals and LR genotypes as the ’600k’ dataset (Table S3).

Genetic analyses

To compare our findings to published results (Melchinger ), we calculated nucleotide diversity (π) on a per-site basis for the 50k dataset using vcftools 0.1.17 (Danecek ). We removed monomorphic sites within each LR-DH pair and calculated π of the remaining n polymorphic sites within each landrace accession. We then compared in each LR-DH pair two vectors (LR, DH) with n elements of π using Mann-Whitney-Wilcoxon tests (unpaired, two-sided; Mann and Whitney 1947).

We used the R package SNPRelate (Zheng ) to conduct a principal component analysis (PCA) for the 50k dataset to investigate the relationship between LR and DH. Furthermore, we calculated principal components in windows in a region around a putative inversion on chromosome 3 of the 600k LR set of accession BU using the R package lostruct (Li and Ralph 2019) with 500 SNPs per window and genome-wide using SNPRelate (Zheng ).

We compared allele frequencies between DH and LR populations. First, we defined minor alleles using the pooled set of all DH and LR accessions and classified them as alternative alleles. Then, we calculated allele frequencies as counts of the alternative allele for each population using plink 1.9 (Chang ).

We determined haplotypes and their respective frequencies in each population by concatenating SNPs in 50kb non-overlapping windows for the 600k dataset. After removing windows with SNPs, haplotypes contained on average 21.2 SNPs. We repeated this analysis with windows based on genetic distance of 0.2 cM for the 50k data without filtering for short haplotypes due to low genomewide SNP density in the 50k data, resulting in windows with on average 8.7 SNPs. We identified the most abundant haplotype in each window in the LR, and classified haplotypes as ’lost’, ’fixed’ or ’segregating’ according to their frequency in the DH. Haplotype diversity was calculated for all windows with at least two haplotypes in LR as where is the haplotype frequency determined in each bin and N the sample size, after Nei and Tajima (1981).

Ancestral allele frequencies and joint probabilities

To distinguish between drift and selection as causes of allele frequency change during the process of DH production, we used a maximum likelihood grid search. We estimated ancestral frequencies and confidence intervals for the joint frequency spectra of DH and ancestral frequencies. Genotyped landrace individuals were sampled from the accession independently of those individuals that gave rise to the DH lines (Mayer ; Melchinger ). Because the landrace individuals are thus not the direct parents of the DH lines, a simple binomial sampling from the LR to DH to estimate confidence intervals would not be appropriate. Therefore, we considered for each accession three binomial sampling events in our estimation. From an ancestral landrace population, a first set of samples was genotyped (LR), and a second set used to produce DH lines from which another subsample was genotyped (DH) (Figure S4). This sampling scheme represents a conservative approximation of a complex production workflow for DH lines described in Melchinger . For each accession and site of the 50k data, we estimated likelihoods across a grid of 100 possible ancestral frequencies ranging from 0.01 to 0.99 as the product of these three binomial probability mass functions, defined as , with n trials, k successes and representing probabilities for three different sampling events, namely (1) the surviving DH lines from haploid induction until genotyping stage, created from the accession (), (2) the genotyped DH samples () and (3) the genotyped landrace samples ().

For one accession this is a matrix with elements for the j-th ancestral frequency and i-th surviving DH individuals and can be estimated withwith d as the DH-allele count and l as the LR-allele count of the q-th site, and and as DH and LR chromosome counts. By maximizing the surface we obtain the maximum likelihood for the ancestral frequency’s probability .

Similarly, we computed a 95% confidence interval by estimating a vector of probabilities for ancestral frequencies for the s-th DH line allele count and the i-th surviving DH individual byfor each site (Figure S4). We used the central 95% probability density of this distribution to define confidence intervals and defined SNPs outside of this confidence interval as allele frequency outliers (aSFS outliers). In brief, the aSFS test informs if a given DH-allele frequency is outside of a 95% confidence interval calculated from LR and DH-line allele frequencies.

We computed a second test statistic to infer outlier SNPs based on the joint probability for a given allele frequency in each population. Here, we computed the joint probability of landrace genotyping, DH line survival and DH line genotyping for each site. We model simple binomial sampling from an ancestral population with unknown allele frequency x which follows a beta distribution with parameters . Integrating over this unknown frequency and using the notation above, the joint probability of observing d and l becomes:with and for each site and accession. We defined SNPs with joint probability in the top 5% of the as outliers in each accession. In essence, outlier alleles identified with the joint probability test are maximally deviating from expected frequencies without selection.

Functional characterization of outlier SNPs

We hypothesized that outlier SNPs are deleterious and recessive and thus should be mostly found in heterozygous genotypes. To test this, we first asked whether outlier SNPs are deleterious and recessive by investigating their level of heterozygosity in the LR. We estimated genotype frequencies for all SNPs of the 50k dataset and removed non-outlier SNPs in linkage disequilibrium () using plink 1.9 (Chang ). We then compared equally-sized samples of outlier SNPs and non-outlier SNPs in 10 frequency bins between 0 and 1.

To study whether outliers affect functional phenotypic trait variation more than random sites, we computed trait effect sizes using a BayesB (Meuwissen ) genomic prediction model implemented in GCTB (Zeng ). We used arithmetic means over four locations of published phenotypes from 351 individuals of six European landrace DH line libraries (GB, RT, SF, Campan Galade, Walliser, Satu Mare) and 53 elite flint lines (Brauner ). We calculated effect sizes based on the pooled dataset of these populations and additional parameters–chain-length, 30000–burn-in 5000. We used BEAGLE 5.0 (Browning and Browning 2009) to impute missing data after filtering based on the same cutoffs as the 50k dataset, resulting in 37884 SNPs for 404 individuals. We then performed an ANOVA using absolute effect sizes for seven traits (shoot vigor, female flowering, Fusarium ear rot resistance, plant height, oil content, protein content and grain yield) separately with the following linear model:with s as the effect of the j-th SNP-type (outlier/non-outlier SNP), α as the effect of the binned (10 bins) k-th frequency in LR, as the interaction effect between frequency bin and SNP-type, and ϵ as the residual effect.

To investigate the potential fitness consequences of outlier loci, we used published genomic evolutionary rate profiling (GERP) scores (Davydov ; Rodgers-Melnick ) estimated by Wang from a multi-species whole-genome alignment of 12 plant genomes and the maize B73 reference genome (AGPv3) and corrected for reference genome bias (Wang ). We used CrossMap 0.2.8 (Zhao ) to update positions of GERP scores to version four of the maize B73 reference genome (AGPv4) (Jiao ).

We compared outlier SNPs to an equally-sized set of non-outlier alleles stratified to match the ancestral allele frequency distribution of outlier SNPs by sampling from bins of 0.1 ancestral allele frequency. Then we calculated genetic distances in 1 Mbp windows based on a maize genetic map (Ogut ) across the genome and partitioned SNPs into recombination quantiles. Lastly, we sampled equally sized fractions of outliers and non-outliers in recombination quantiles to account for differences in expected genetic load in genomic regions with different recombination rates. We calculated load based on observed genotypes under an additive and a recessive models for each SNP of the 50k dataset. We then calculated the sum of overlapping GERP scores in 1 cM windows around each SNP for both models.

Data availability

Scripts used for this project are stored in a public GitHub repository (https://github.com/LZeitler/eurodh-scripts). Supplemental material available at figshare: https://doi.org/10.25387/g3.12382577; https://doi.org/10.6084/m9.figshare.12456308.

Results

DH lines show decreased genetic diversity compared to landrace populations

We first evaluated the population structure of our samples using a genome-wide principal component analysis (PCA). Groupings largely follow overall expectations, with doubled-haploid (DH) and landrace (LR) individuals clustering well by accession on the first two principal components. The third principal component, however, separates a subgroup of the DH individuals from the main RT cluster (Figure 1A and S5).

Figure 1

(A) Principal component analysis for the DH and LR of the 50k dataset. (B) Average per site nucleotide diversity at polymorphic sites for unimputed 50k data in each of the five LR-DH pairs. Means, represented by diamonds, are significantly different from each other within accessions (). (C) The joint site frequency spectrum (jSFS) for DH and LR populations. Allele counts are based on the published filtered dataset (Melchinger ).

Using a set of genotyped SNPs (the 50k dataset, see Methods) from Melchinger , we compared per-site nucleotide diversity (π) in individual accessions (LR-DH pairs). We found that average π is significantly different between the LR and DH (two-sided Mann-Whitney-Wilcoxon for all comparisons). In virtually all comparisons the nucleotide diversity is higher in LR than in DH (Figure 1B). In the accession RT, however, π is lower in the LR compared to DH, likely reflecting the population substructure observed among the DH (Figure 1A).

To investigate differences in diversity in more detail, we compared the joint site frequency spectrum (jSFS) of the two populations for each accession (Figure 1C). The jSFS shows substantial variation around the expected 1:1 line, highlighting differences in allele frequencies between the two populations. In particular, the jSFS reveals a number of alleles segregating in the LR that have been lost in the DH lines, as well as a smaller number that have been fixed. In some cases, we found segregating SNPs in the DH lines that were monomorphic in the respective LR population, suggesting that the frequency of the minor allele in the ancestral population must have been low and was not sampled in the genotyped set of LR individuals. Comparison of our jSFS with results from Melchinger reveals striking differences (Figure S6B), as these authors appear to have filtered out a substantial proportion of alleles at low frequency in one or both populations, removing much of the observed signal of allele frequency change.

Genome-wide pattern of allele frequency distortion

The jSFS highlights the dramatic difference in allele frequency between LR and DH populations. But, because DH lines were not the direct offspring of genotyped LR individuals (Figure S4), a direct comparison between LR and DH populations is not straightforward. To circumvent this difficulty, we first inferred allele frequencies in the ancestral population of both LR and DH (see Methods). Under neutrality, the expected allele frequency in the DH population is identical to that of the ancestor, and comparison of the site frequency between these populations (ancestral site frequency spectrum test; aSFS test) allows identification of alleles with unusual shifts in frequency (Figure 2). The number of outlier SNPs in the aSFS test (those outside of the 95% confidence interval around the ancestral frequency) varied between 1769 (4.66%) in accession SF and 6364 (16.78%) in accession RT. In total, we identified 12345 distinct outlier loci in the five accessions. Of these, 9305 were detected in only one of the accessions and 1877, 539, 317, 307 SNPs overlapped in 2, 3, 4, 5 accessions, respectively (Table S4). A substantial fraction (14.71%) of these outlier loci were segregating in LR but lost in DH, while only a very small minority of outlier loci (1.45%) were fixed during DH production. In addition, sites that were outliers in multiple accessions had a higher mean frequency in LR and lower mean frequency in DH populations compared to outliers that were unique to one population (Figure S7A-B). Similarly, we found that a frequency reduction was more common in aSFS outliers shared among accessions while on average unique outliers changed only little (Figure S7C).

Figure 2

(A) Estimated ancestral and DH allele frequencies for all accessions of the 50k dataset show significant outliers, with the 95% confidence interval represented by blue lines in the joint frequency spectrum (aSFS test). Percentages indicate the proportion of SNPs above, below and inside the interval. (B) Joint probability test along chromosome 3 in the Bugard landrace (BU). Colored dots represent the top 5% -values which we defined as outliers. Colors represent the number of accessions in which a given locus is an outlier. The dashed line indicates the centromere position.

As a second means to identify outlier loci, we calculated the joint probability of the LR and DH genotypes (see Methods), yielding a p-value for each SNP (Figure S8). The test identified 9458 outliers across accessions. We compared these outliers among different accessions and found that highly significant outlier SNPs were often shared among several populations (Figure S9).

The two approaches yielded largely similar results. Although the aSFS test identified nearly twice as many outliers as the joint probability test, virtually all of the outliers (97%) identified in the joint probability test were also found in the aSFS comparison (Figure S10). Some clusters of co-located outlier sites were evident in both tests (Figure 2B). For instance, we detect a large region close to the centromere of chromosome three with multiple loci that were detected as outliers in multiple accessions.

Haplotype tests reveal selection hotspots

Tests based on individual SNPs yielded outliers and indicated that these are often shared between accessions. However, such tests are limited in their power to detect changes in low frequency alleles. Haplotype based comparisons are more sensitive to rare alleles, as every new allele creates an additional haplotype. We imputed the 600k SNPs in the DH data using genotypes from the LR populations, and identified haplotypes in non-overlapping 50kb windows.

We observed a reduction in the number of haplotypes in the DH compared to the LR population: in the 600k data, we identified on average 7.65 segregating haplotypes per window in the LR and 4.36 in the DH. Haplotype diversity was significantly reduced in the DH compared to the LR populations (0.40 compared to 0.63; Figure 3A). The difference in haplotype diversity between LR-DH pairs was less pronounced in the 50k dataset, but haplotype diversity of the DH populations was also reduced in every accession compared to the LR (Figure S11 and S12A).

Figure 3

(A) Comparison of average haplotype diversity in 50kb windows for the 600k dataset in different accessions. Haplotype diversity of the imputed DH dataset is reduced compared to the LR dataset. (B) LR haplotype frequency in BU along chromosome 3, colored by the fate of the haplotype in the DH population. The centromere is shown as a vertical dashed line. Percentages listed in the legend correspond to genome-wide proportions of fixed, lost and segregating haplotypes.

We tracked the frequency change of the most common haplotype in each window between LR and DH populations in the 600k dataset and classified haplotypes according to their fate in the DH. While the majority (13607 or 58.35%) of the most common haplotypes in the LR were still segregating in the DH populations across accessions, a substantial minority (5113, 21.93%) were lost and another large fraction (4600, 19.72%) were fixed (Figure S13, Table S5).

Similar to the SNP outlier tests, we found several genomic regions with multiple consecutive windows exhibiting fixation or large changes in haplotype frequency (Figure S13). Windows with losses of major haplotypes coincided with highly significant joint probability outliers (Figure S14). In particular, the same region of chromosome 3 identified by SNP outliers showed strong signals of haplotype change in BU and RT, where even haplotypes with intermediate frequencies in LR were fixed in the DH (Figure 3B). To further investigate this region, we conducted a local PCA (Li and Ralph 2019) in the BU landrace data between positions 65 Mb and 95 Mb, revealing three distinct clusters across multiple consecutive windows (Figure S15A). This contrasts with genome-wide PCA (Figure S15B), and is consistent with previous reports of a segregating inversion polymorphism in this region (Romero Navarro ).

To exclude potential effects of imputation errors, we performed the same analysis with the 50k data and windows based on genetic distance. The estimated fraction of fixed alleles was reduced and the putative inversion on chromosome 3 was no longer detected (Figure S12B), likely due to lower SNP density and imputational quality for minor alleles (Figure S3, supplementary note). Nonetheless, lost alleles occurred at similar rates (10.03%, Figure S12B), confirming the results of the 600k data.

Outliers are more heterozygous than random alleles

To characterize potential changes due to selection, we first asked whether outlier SNPs are more likely to be recessive compared to random sites. We hypothesized that recessive deleterious sites should show higher observed heterozygosity in the outbred LR as selection should effectively remove homozygous genotypes. Indeed, outliers had a significantly higher frequency of heterozygotes in the LR compared to a frequency-matched sample of non-outliers for all populations (; Figure 4 and S16).

Figure 4

Violin plots for the frequencies of heterozygous genotypes of LD-pruned non-outlier SNPs and outlier SNPs in the LR accessions for aSFS outliers. Comparisons with asterisks have significantly different means (1000 bootstraps, ).

Next, we investigated whether outlier regions were enriched for functional variants. To test this, we estimated polygenic effect sizes for seven traits from a DH line panel with 404 individuals from different landraces using a BayesB prediction model (Meuwissen ). We then tested whether effect sizes differed between outlier and non-outlier windows across a range of allele frequencies. Allele frequency was highly significant for all traits, and while outliers had significantly different absolute effect sizes for only two traits (plant height, oil content) we found significant interactions between SNP-type and allele frequency for all traits (Table S6). The relationship between effect size and allele frequency was negative for outlier SNPs, but positive for non-outlier SNPs across all traits after fitting linear regression models (Figure S17).

Lastly, we aimed to characterize changes in genetic load due to selection using published GERP estimates of evolutionary constraint at each SNP calculated from a phylogeny of 13 species (Wang ). Previous studies in maize have shown that GERP scores correlate with estimated SNP effects on yield and are thus a quantitative proxy for the fitness effects of a locus (Yang ). We analyzed whether outliers contribute higher genetic load than random SNPs by summing GERP scores in the 1 cM centered around each SNP. Under a recessive model, in three out of the five LR and four DH populations outlier windows showed lower genetic load than random windows, while there was no significant difference in BU in both the LR and DH (t-test, ; Figure S18). Under an additive model, the mean load of outliers was significantly lower compared to non-outlier in all accessions in the LR and in all but BU in the DH (t-tests, ; Figure S18).

Discussion

Reduction in genetic diversity between landrace and DH populations

We find a significant reduction in genetic diversity during DH production in four out of five accessions (Figure 1B and 3A). The increase in π seen in RT may have resulted from the observed population sub-structure in the DH population, perhaps due to the use of seed from distinct rounds of regeneration ex situ (Chebotar ). The RT landrace itself also exhibited higher homozygosity than samples of other LR populations, suggesting a history of inbreeding during conservation. Whatever the cause, our estimates of diversity at the haplotype level reinforce these overall findings, showing even greater losses of diversity than seen at individual SNPs (Figure 3A and S11). Altogether, these results closely follow theoretical predictions regarding the consequences of inbreeding (Charlesworth and Willis 2009; Schnable and Springer 2013).

The loss of diversity we observe in DH populations stands in contrast to previous findings (Melchinger ) using the same data. For a detailed comparison, we reconstructed the jSFS using the original data from the previous study, revealing that Melchinger had filtered the data in such a way as to remove sites with extreme allele frequencies in either population (Figure S6B, Melchinger ). While minor allele frequency filters are often applied in quantitative genetic studies, the removal of rare alleles can strongly influence results of population genetic analyses (Weale 2010; Linck and Battey 2019). Moreover, such alleles are of particular interest for the conservation of genetic diversity. Therefore, we limited our filtering to data quality but did not remove rare alleles (Table S2).

DH production creates selection hotspots

To understand the effect of DH line production on the change in diversity across the genome, we employed two outlier tests comparing the allele frequency changes between the LR and DH. These tests identified loci for which the allele frequency shifted more than expected by random drift. Outliers exist in all five European landrace accessions (Figure 2A). While the aSFS test resulted in a larger set of outliers than the joint probability test, both tests identified a largely overlapping set of outliers (Figure S10). And while many outliers were shared among accessions (Figure 2B and S9), indicating some shared signal resulting from DH production, the majority of outliers were accession specific. In other crops like potato, it has been shown that the genomic signals of inbreeding are largely specific to individual lines (Zhang ). The increased strength of selection due to the instantaneous homozygosity during DH production and the shared history of European maize landraces might have caused the increased signal of shared outliers among accessions.

While we found outliers distributed across the whole genome, we also observed clustering in specific genomic regions (Figure 2B). The distribution of the fate of major haplotypes in windows along the genome revealed regions enriched for outliers that go to fixation or loss (Figure 3 and S11). One of the most pronounced signals was on chromosome 3 in BU and in RT. This approximately 25-Mb region (70-95 Mb) overlaps with a previously identified putative 6-Mb inversion that is associated with flowering time in maize (Romero Navarro ). Further testing using local principle components indicated the presence of this inversion in the landrace sample of accession BU (Figure S15). In this landrace from Southern France, the inversion may be involved in flowering time adaptation. Alternatively, unconscious selection on flowering time might have occurred during haploid induction to synchronize landrace flowering with the inducer line or subsequent cultivation in northern latitudes. Other regions where outliers clustered in longer windows in all accessions (Figure S13) were mostly located in low recombination regions around centromeres (Ogut ). Weakly deleterious alleles are likely to accumulate in such regions (Rodgers-Melnick ; Yang ), and if most fitness-affecting mutations are at least partially recessive (Yang ) such regions might be expected to show selection when made homozygous during DH production.

Imputation of the 600k data strongly increased the marker density and therefore the power of our study. Yet, minor alleles were more affected by imputational error (Figure S3), which might introduce a potential bias toward a lower number of estimated haplotypes per window and therefore on haplotype diversity. Hence, we also performed the haplotype analysis using the 50k dataset. Despite the lower power due to the sparse marker density, the 50k data showed mostly similar patterns for haplotype diversity and major haplotype change. One difference was the lack of evidence for the putative inversion (Figure S12B). While the 600k data were designed with a more diverse sample of maize lines that included both haplotypes at this inversion, the region on the MaizeSNP50 chip does not include SNPs segregating for the inverted haplotype (see Pyhäjärvi ). Low frequency polymorphisms should have only little effect on the fate of the most abundant haplotype during DH production, except for fixing major haplotypes in windows where all haplotypes were rare in the LR and minor alleles were incorrectly imputed. The higher values, i.e., outliers (based on unimputed data) overlapping with lost major haplotypes confirm that the fate of the major haplotype was little affected by imputation errors (Figure S14). We expect that most rare alleles in the DH should have been correctly imputed as we used the LR that they were derived from as reference for the imputation algorithm. In the absence of selection, alleles in the DH populations should exist at similar frequencies in their respective LR. Hence, potential imputation errors rather lead to an underestimation of the selective loss of alleles.

Differences between DH and LR are potentially due to recessive deleterious load

Doubled-haploid lines show particularly poor fitness compared to outcrossing lines and modern inbred lines (Strigens ; Böhm ). The observed inbreeding depression in DH and inbred lines in maize is likely due to accumulation of deleterious alleles as a result of inbreeding (Bataillon and Kirkpatrick 2000; Charlesworth and Willis 2009). Recent work has shown that the observed decrease of heterozygosity during inbreeding of maize landraces is slower than expected, suggesting that the exposure of recessive deleterious alleles removes certain haplotypes and maintains heterozygosity (Roessler ). While multiple cycles of inbreeding allow for recombination and purging of genetic load, DH production induces instantaneous homozygosity which reduces the possibility for effective purging. Moreover, heterozygosity can be regenerated by intercrossing DH lines only in regions of the genome where multiple haplotypes were captured during DH production. Our results show that this is not the case in >20% of the genome, where even the previously most common haplotype was lost (Figure 3B and S12B). For lower frequency haplotypes this effect is expected to be even more severe. Hence, a substantial number of haplotypes will not be available for further selection in intercrossed DH line populations.

Consistent with the assumption of recessive deleterious mutations, we find evidence that outlier SNPs are more likely to be heterozygous in LR populations (Figure 4 and S16). While outlier sites show strong shifts in allele frequency, they are unlikely to be directly causal themselves, but instead linked to deleterious sites. The comparison of genetic load around outliers and non-outliers, identified significantly lower load around outlier SNPs in both LR and DH (Figure S18), which might indicate that the causal deleterious SNP in the region was not genotyped in our analysis. The relatively low marker density of the SNP chips, in addition to the ascertainment bias towards inbred breeding material, likely prevented assaying the causal deleterious alleles, but allowed us to capture the loss in diversity around causal sites. Furthermore, the estimated genetic load differences between outlier and non-outlier are mainly driven by the number of loci with annotated GERP scores overlapping with the investigated SNPs in each 1 cM window, i.e., the number of SNPs the load is calculated with. Because we stratify the SNPs for proper comparison based on recombination rate and allele frequency we can rule out uneven representation of outlier/non-outlier in genetic load rich pericentromeric regions, compared to more distal, gene-rich chromosomal regions. Consequently, a genomewide underrepresentation of outlier alleles in conserved, gene rich regions can be another interpretation of reduced outlier load in LR.

Finally, we searched for evidence that outlier loci were particularly likely to contribute to phenotypic variation. We find evidence of selection in plant height and oil content (Table S6), despite our relatively simple additive GWAS model and the fact that most loci show at least partially recessive effects on yield (Yang ). Negative slopes (Figure S17) in the regression models of allele frequency and absolute effect size are consistent with the deleteriousness of outlier SNPs and match previous findings in maize that deleterious mutations explain a significant proportion of phenotypic variance (Yang ). The significant effect of outliers on plant height effect sizes might be the signature of natural or unintended selection on height or a related phenotype in maize landraces. While oil content is sometimes used as a means to identify haploid seed in the creation of DH lines (Prigge and Melchinger 2012), a different approach was used to create the DH lines used here (Melchinger ), and we have no explanation for the signal of selection on oil content observed here.

Overall, our results suggest that the observed reduction in diversity within different populations is not caused by a few large-effect loci, but rather by a polygenic effect of partially recessive, mildly deleterious mutations (Bataillon and Kirkpatrick 2000).

Conservation of landrace diversity

Landraces are an invaluable source of adaptive diversity (Bellon ; Gates ), and their conservation should remain a high priority for future generations. Here, we showed that DH line libraries from landraces do not capture the full diversity present in the landrace. Therefore, while DH line libraries present a valuable tool to introgress known alleles into breeding programs, we conclude they can not replace ex situ and in situ conservation efforts. To preserve landraces and their full genetic diversity, they should be reproduced in large populations to prevent inbreeding and the consequent shift of allele frequencies. An improved understanding of inbreeding and the underlying genomic changes it induces will help to conserve these genetic resources and harness their diversity to breed improved crop varieties.