Introduction

Transposable elements (TEs) are an intrinsic part of eukaryotic genomes and their evolution [1–12]. In addition to modulating genome size, the ability of TEs to create genetic diversity through insertion and excision events can lead to new phenotypes on which selection can act. TEs can alter phenotypes through various mechanisms, including the functional disruption of genes [1, 2], large-scale changes in the regulatory apparatus [3, 4], alteration of epigenetic landscapes [5, 6], ectopic recombination and structural rearrangements [7, 8]. In plants, the dynamics of TE loss and proliferation play a major role in genome evolution [e.g., 9-12]. TEs therefore constitute potentially important drivers of plant evolution, both in nature and during domestication [13].

Beyond their influence on genome structure, and given that their transpositional activity can be stress-inducible [for review 14], TEs are often regarded as more likely than classical point mutations to produce the diversity needed for individuals to respond quickly to challenging environments [15–17]. For instance, punctual TE polymorphisms can lead to gains of fitness and evolve under positive selection [2, 18–23]. TE polymorphisms can even induce more extreme changes in gene expression than single nucleotide polymorphisms (SNPs) in plants [24, 25].

Despite such evidence, whether TE polymorphisms are major contributors to adaptation to changing environments is still debated. Indeed, TE transposition can be disruptive, and purifying selection has been shown to play an important role in TE evolution [e.g., 26, 27]. Based on simulations, it has been suggested that the persistence of TE polymorphisms within a genome without an uncontrolled accumulation, can only be achieved if weak purifying selection is the main force governing TE evolution [28–31]. The uncertainty surrounding this important question in evolutionary genomics results from the limited number of studies that comprehensively tested the extent to which selection shapes TE allele frequencies, both in plants [23, 32] and animals [26, 33–36] and characterized the distribution of fitness effects of new TE insertions. To clarify this question, we used the plant model system Brachypodium distachyon [37] to disentangle the effects of purifying and positive selection on TE polymorphisms in natural populations.

B. distachyon is a wild annual grass endemic to the Mediterranean basin and Middle East. Recent genetic studies based on more than 320 natural accessions spanning from Spain to Iraq (hereafter referred to as the B. distachyon diversity panel) revealed that B. distachyon accessions cluster into three main genetic lineages (the A, B and C genetic lineages), which further divide into five main genetic clades that display little evidence for historical gene flow (Fig. 1A; [38, 39]). Niche modeling analyses suggest that the species moved southward during the last glacial period and recolonized Europe and the Middle East within the last five thousand years [39]. Consequently, while some B. distachyon genetic clades currently occur in the same broad geographical areas (Fig. 1A), natural accessions are adapted to a mosaic of habitats [38, 39]. These past and more recent shifts in the species distribution led to clear footprints of positive selection in the genome [39, 40] and make B. distachyon an ideal study system to investigate the contribution of TEs to the adaptation of plants in the context of environmental changes.

Fig 1.

In B. distachyon, TEs are exhaustively annotated and account for approximately 30% of the genome [37]. Recent TE activity has been reported for many families, but despite past independent bottlenecks and expansions experienced by the different genetic clades, no lineage-specific TE family activity has been observed [32]. Rather, TE activity tends to be homogeneous throughout the species range and across genetic clades, indicating a high level of conservation of the TE regulatory apparatus [32]. While purifying selection shapes the accumulation patterns of TEs in this species [32], some TE polymorphisms have been observed in the vicinity of genes [32], potentially affecting gene expression [41]. These early studies, based on a relatively small number of accessions originating exclusively from Spain and Turkey, suggested that TE polymorphisms could contribute to functional divergence and local adaptation in B. distachyon [32].

To test this hypothesis, we used the B. distachyon diversity panel to identify TE polymorphisms in a large set of 326 natural accessions spanning the whole species distribution. We combined a set of population genomic analyses to assess the proportion of TE polymorphisms associated with positive or purifying selection as well as neutral evolution. We also quantified the strength of purifying selection through forward simulations. Altogether, our work provides the first quantitative estimate, to our knowledge, of the adaptive, neutral, and disruptive potential of TEs, while accounting for changes in TE activity, in a plant harboring a relatively small genome. Altogether, our result advocate against an extended role of TEs in recent adaptation.

Results

Genetic variation in Brachypodium distachyon

Using the B. distachyon diversity panel (Fig. 1A), we identified 97,660 TE polymorphisms in our B. distachyon dataset, of which 9,172 were retrotransposons, 52,249 were DNA-transposons and 36,239 were unclassified. We also identified 9 million SNPs across the 326 samples, including 182,801 synonymous SNPs. A Principal Component Analysis (PCA) performed either with SNPs or TE polymorphisms reflects the previously described population structure of B. distachyon [38, 39], with the first two components of the PCA splitting the data according to the demographic structure (Additional file 1: Fig. S1). Investigating the genetic variation caused by retrotransposons and DNA-transposons revealed that the observed diversity in retrotransposons strongly correlated with the demographic structure (Mantel test; r = 0.79, p value = 0.001), while the observed diversity in DNA-transposons only had a weaker correlation (Mantel test; r = 0.36, p value = 0.001) with the demographic structure (Additional file 1: Fig. S2).

From the initial TE and SNP dataset, we could estimate the time of origin in generations (age) of 50,891 TE polymorphisms and 108,855 synonymous SNPs based on pairwise differences in identity by descent (IBD) regions around the focal mutation (see Materials and methods). The results of the age estimate analysis were checked by contrasting the observed correlation between allele age and frequency of synonymous SNPs to the theoretical predictions of Kimura and Ohta [42] for neutrally evolving mutations. We found that, the observed correlation matched expectations (Additional file 1: Fig. S3), with older alleles found on average at higher frequencies than younger ones. Furthermore, most TE polymorphisms in our dataset were young and only a few were very old (Additional file 1: Fig. S4).

The overall contribution of TEs to clade differentiation and adaptation is limited

To examine the overall contribution of TEs to evolution and adaptation in B. distachyon, we first identified regions of the genomes that were likely affected by recent selective sweeps. The fast increase in the frequency of a beneficial allele is expected to lead to a longer than average haplotype around the mutation under positive selection. Such events (known as selective sweeps) can be identified by computing the integrated haplotype score (iHS) around focal mutations [43]. We therefore computed iHS along the genome for the four derived genetic clades. Regions of the genomes with significantly higher iHS than average are expected to harbor mutations that were under positive selection during evolution and adaptation. We hypothesized that if TEs constitute an important part of the genetic makeup that led to adaptation in a given genetic clade, then they should be more frequently fixed or at higher frequencies in regions with high iHS than in the corresponding regions that did not experience recent selective sweeps in other clades.

First, we tested if more TE polymorphisms were fixed in a specific region of the genome if a genetic clade had a high iHS, and presumably experienced a selective sweep, than in other genetic clades. An analysis of covariance (ANCOVA) revealed that the number of fixed TE polymorphisms per clade did not significantly differ between high iHS regions and the same regions in other clades (Table 1). These results indicate that there is no correlation between the overall number of fixed TE polymorphisms per clade in a region and recent selective sweeps. However, the number of fixed TEs in genomic regions along the genome was significantly affected by the total number of TEs in the region, the TE superfamily, the TE age, the genetic clade and the overall genetic features of the region (e.g., recombination rate, see Materials and methods) but not by the iHS itself (Table 1). Similarly, we tested if the allele frequency of TE polymorphisms was significantly higher in a specific region of the genome if a genetic clade had a high iHS than in other genetic clades. A second ANCOVA revealed that the allele frequency of TE polymorphisms was significantly influenced by the TE superfamily, TE age, clade and overall genetic features of the region but not by the iHS (Table 2). These results indicate that TEs in high iHS regions did not experience a significant increase in their frequency and that TEs in high iHS regions are experiencing the same selective constraints as other TEs.

Table 1

Table 2

A complementary approach to explore the impact of positive selection on TEs consists in investigating their genetic differentiation among populations. Using the five genetic clades as focal populations, we computed X^tX values, a standardized measure of genetic differentiation corrected for the neutral covariance structure across populations [44, 45], for each TE polymorphism. Mutations affected by positive selection are expected to be over-differentiated between clades and display significantly higher X^tX values than other mutations [45]. In contrast, a low X^tX value implies that the mutation is less differentiated than other mutations and potentially evolves under balancing selection, whereas purifying selection and a neutral evolution are not expected to impact the differentiation of a mutation among populations [44]. We contrasted the observed X^tX values computed for each TE polymorphism to a simulated pseudo-observed dataset (simulated observations under the demographic model inferred from the covariance matrix of the SNP dataset, for more details see [45]) and found that only a small fraction of the TE polymorphisms (0.06%) displayed X^tX values higher than the 97.5% quantile of the simulated values (Fig. 1B). This indicates that only a few TE polymorphisms are over-differentiated among genetic clades and might have been affected by positive selection. However, a relatively larger portion of the TE polymorphisms (4.3%) displayed X^tX values smaller than the 2.5% quantile of the simulated values (Fig. 1B), indicating that balancing selection might also shape TE frequency in B. distachyon.

To further examine the contribution of TEs to adaptation, we tested whether and how many TE polymorphisms were significantly associated with environmental factors. If the presence of a TE provides an advantage in a certain environment and contributes to adaptation, we expected a correlation between the environment and the presence/absence of this TE. In this context, we performed genome-environment association analyses (GEA) using all TEs and SNPs identified across the 326 samples and 32 environmental factors associated with precipitation, solar radiation, temperature, elevation and aridity (see in Materials and methods for the full list). The GEA revealed that only nine of the 97,660 TE polymorphisms were significantly associated with some environmental factors (Additional file 2: Table S1), confirming that TEs only had a limited contribution to adaptation in B. distachyon. Importantly, two of these nine TEs were found in a gene, and three were in the vicinity of genes (less than 2 kilobase (kb) away, Additional file 2: Table S1).

Purifying selection dominates the evolution of TE polymorphisms in B. distachyon

To further characterize the forces governing the evolution of TE polymorphisms in B. distachyon, we examined the genome-wide frequency distribution of TEs. We first computed the folded site frequency spectrum (SFS) and found that the folded SFS of TE polymorphisms was shifted toward a higher proportion of rare minor alleles compared to neutral sites in all genetic clades (Fig. 1C-G). Splitting the TE data into DNA-transposons and retrotransposons resulted in similar folded SFS and shifts in both TE classes (Additional file 1: Fig. S5 and S6).

These shifts could be the result of purifying selection as the analyses presented above indicate that positive selection has a negligible effect on TE polymorphism frequencies in B. distachyon. However, in contrast to SNPs, TEs do not evolve in a clock-like manner, as their transposition rate is known to vary between generations [46, 47]. Changes in transposition rate and purifying selection can lead to similar shifts in the SFS but can be disentangled using age-adjusted SFS [48]. In brief, if TE polymorphisms are evolving neutrally, they are expected to accumulate on average at the same rate in a population as neutral SNPs of the same age. Hence, Δ frequency, the difference between the average frequency of TE polymorphisms and neutral sites in a specific age bin, will remain close to 0 regardless of the polymorphisms’ age. In contrast, if TE polymorphisms evolve under purifying selection, they will tend to occur at lower frequencies than neutral SNPs of the same age, as selection will prevent them from accumulating in the population. Consequently, the Λ1 frequency will reach negative values for older TE polymorphisms [48].

Because this model does not allow for back mutations, as typically observed for DNA-transposons that can excise from the genome, we primarily investigated the age-adjusted SFS of retrotransposons in the four derived clades. This analysis revealed that retrotransposons are indeed prevented by natural selection from randomly accumulating, as older retrotransposons are significantly less frequent than neutral SNPs of the same age (Fig. 2; one-sided Wilcoxon test, Bonferroni corrected p value < 0.01).

Fig. 2

As previous studies showed that the distance between TE polymorphisms and the next gene can impact the strength of selection affecting TEs [6, 49, 50], we further split our retrotransposon polymorphisms into three categories based on their distance to the next gene: retrotransposons (i) in and up to 1 kb away from genes, (ii) between 1 kb and 5 kb away and (iii) more than 5 kb away. The age-adjusted SFS of all three categories displayed the same pattern as that observed for the whole retrotransposon polymorphism dataset: older retrotransposon polymorphisms were significantly less frequent than neutral sites of the same age regardless of their distance to genes (one-sided Wilcoxon test, Bonferroni corrected p value < 0.01), indicating that retrotransposons more than 5 kb away from genes are also affected by purifying selection (Fig. 2).

Retrotransposon polymorphisms tended to be more deleterious than SNPs predicted to have a high impact on fitness. Indeed, the age-adjusted SFS of retrotransposons resulted in a larger deviation of Δ frequency from 0 than for non-synonymous SNPs and high effect SNPs (Fig. 2). In addition, Δ frequency in the oldest (last) age bin was significantly smaller than in all other age bins in the A_East, B_East and B_west clades (one-sided Wilcoxon test, Bonferroni corrected p value < 0.01). In the A_Italia clades the oldest age bin was not significantly different form the second oldest age bin (two-sided Wilcoxon test, Bonferroni corrected p value N.S.). While older non-synonymous SNPs and high effect SNPs were generally less frequent than neutrally evolving SNPs at the same age, the negative Δ frequency trend was reversed for the oldest non-synonymous SNPs and high effect SNPs (Fig. 2). In all clades, Δ frequency in the oldest age bin was significantly higher than at least the lowest Δ frequency observed in the other age bins for non-synonymous SNPs, as well as high effect SNPs (one-sided Wilcoxon test, Bonferroni corrected p value < 0.01). This might be because not all predicted non-synonymous SNPs and high effect SNPs might result in fitness effects. Those SNPs can therefore evolve neutrally or nearly neutrally and persist as polymorphic SNPs much longer in a population than those affecting fitness negatively. Hence, the last age bin of the non-synonymous and high effect SNP age-adjusted SFS likely harbors mainly neutrally and nearly neutrally evolving mutations, and consequently, Δ frequency in not the smallest in the last age bin in these age-adjusted SFS.

To assess whether similar forces may drive retrotransposon and DNA-transposon evolution, we repeated the analysis for DNA-transposons. The age-adjusted SFS of DNA-transposons revealed very similar patterns, with Δ frequency showing significant deviations from 0 in older age bins (one-sided Wilcoxon test, Bonferroni corrected p value < 0.01) and DNA-transposon polymorphisms being more deleterious than non-synonymous SNPs and high effect SNPs (Additional file 1: Fig. S9).

Forward simulations allow us to quantify the strength of purifying selection

To evaluate to what extent the proportion of neutrally evolving mutations in the focal group of mutations affects the shape of the age-adjusted SFS, we ran forward simulation with mutations under multiple selective constraints, and we tested what ratio of neutral to selected mutations can lead to an age-adjusted SFS similar to that observed for retrotransposons in B. distachyon. Specifically, we investigated the conditions under which we observed a Δ frequency in the oldest age bin significantly smaller than Δ frequency in all other age bins. Our simulations revealed that the shape of the age-adjusted SFS of retrotransposons could only be reproduced if less than 10% of the mutations were neutrally evolving for most of the selective constraint investigated (Additional file 1: Fig. S7 and Additional file 2: Table S2).

Finally, we used the results from our simulations to narrow down the selection strength affecting retrotransposons in B. distachyon by investigating the age of the oldest retrotransposons in our dataset. The main difference between the age-adjusted SFS of mutations evolving under weak and strong purifying selection is that the oldest mutations are much older in the simulation with weak purifying selection than in the simulation with strong purifying selection. This age difference arises because mutations under strong purifying selection are removed from the population more effectively and, therefore, cannot persist as long in the population. Examining the age of the last retrotransposon bins in the age-adjusted SFS revealed that the ages of the oldest retrotransposons were the most similar to the expected ages of the oldest mutations in our simulations, with a scaled selection coefficient (S) of –5 and –8 (Fig. 3), indicating that retrotransposons in B. distachyon are under moderate purifying selection. In simulations with a nearly neutral selection coefficient (S = –1), the simulated mutations were much older than the oldest observed retrotransposons (Fig. 3). Conversely, in simulations with a strong purifying selection coefficient (S < –10), they were much younger than the oldest observed retrotransposons (Fig. 3).

Fig. 3

Discussion

B. distachyon is a widely used model species in evolutionary genomics, molecular ecology, developmental biology and crop functional genomics [for review 51, 52] with past and ongoing TE movements in its genome [32]. In this study, we used a diversity panel containing next-generation sequencing data from over 320 individuals sampled across the whole geographical range of B. distachyon to examine the role of TEs during evolution and adaptation. We investigated the frequency with which positive selection led to an increase in the frequency and fixation of TEs and quantified the strength of purifying selection on TE polymorphisms. Accounting for population structure and fluctuant transposition rates, we demonstrate that TEs are rarely part of the genetic makeup that was positively selected during environmental adaptation in B. distachyon. Furthermore, we show that the majority of TE polymorphisms found in the natural population of this model species are under weak to moderate purifying selection, with only a small minority of TE polymorphisms evolving neutrally.

Rare instances of positive selection on TEs

By combining complementary approaches, we were able to demonstrate that TEs are rarely the target of positive selection in B. distachyon. We first probed for footprints of positive selection on TE polymorphisms using the five genetic clades as focal populations. In conducting this analysis, we did not find TE polymorphisms to be at high frequencies or fixed at higher rates than expected, in regions of the genome presumably harboring selective sweeps in at least one of the genetic clades (high iHS regions). This suggested that TEs were rarely the target of positive selection, which we confirmed with a genome-wide scan for overly differentiated TE polymorphisms using X^tX analysis. Indeed, this approach revealed that only a very small proportion of TE polymorphisms are more differentiated than expected under a neutral scenario.

Importantly, the X^tX analysis also revealed that a non-negligible fraction of the TE polymorphisms is less differentiated than expected and are shared among genetic clusters. This could be the result of selection favoring the same TE polymorphisms in different accessions to adapt to similar environmental constraints across genetic clades. To test this scenario, we performed GEA with 32 environmental factors, and found only nine TE polymorphisms significantly associated with any of these, and representing a very small proportion (< 0.01%) of all the TE polymorphisms we identified. Interestingly though, these nine TE polymorphisms were associated with environmental variables pertaining to precipitation, temperature and altitude, which are known to drive adaptation in B. distachyon [39]. Some insertions were found within or in close proximity of genes, making these polymorphisms very good candidates for future functional validation.

Single TE insertions can have a drastic impact on phenotypic variation and be affected by positive selection [for review, see 16, 53, 54]. For instance, TEs have increased in frequency through positive selection in humans [23] or during range expansion in Arabidopsis [25] and D. melanogaster [20, 22]. Evidently, B. distachyon exhibits a different pattern, as causal mutations for adaptation in this grass species are rarely TEs. Only a few studies have thoroughly quantified the extent to which positive selection influences the evolution of TEs [25, 26, 28, 30, 31]. But two of these drew similar conclusions to us, in the green anole Anolis carolinensis [30] and in the invasive species Drosophila Suzukii [31]. In addition, a large number of candidate genes for adaptation were identified with a similar approach focusing on SNPs [39], indicating that population structure or demographic events are not limiting factors for the methods we used. Altogether, these observations call for a closer investigation of which forces, e.g., purifying selection or neutral evolution, are important in shaping TE allele frequency in natural populations.

Moderate purifying selection is the dominant force during TE evolution

Our results suggest that purifying selection is an important factor limiting the ability of TE polymorphisms to fix and increase their frequency in B. distachyon. Indeed, one of the significant explanatory variables in our ANCOVA models was the genetic clade, a proxy for the effective population size (N_e), which affects the efficiency with which selection can fix beneficial mutations and purge deleterious ones. In B. distachyon, the number of fixed TE polymorphisms per clade and the frequency of TE polymorphisms were negatively correlated with N_e, indicating that the accumulation of TEs is significantly lower in genetic clades with a larger N_e, potentially because of a greater efficacy of purifying selection.

It is widely accepted that most new TE insertions have a deleterious or no effect on the fitness of the host [22, 26, 28, 30, 33–36, 55]. To properly quantify the effect of purifying selection on TE evolution in B. distachyon, we used age-adjusted SFS analyses to evaluate the selective constraint experienced by TE polymorphisms while accounting for previously reported changes in their activity [32]. While this method can only be applied to retrotransposons (because the model does not allow back mutations), it provided a first clue on the importance of purifying selection on TE evolution and revealed that overall, retrotransposons evolved under purifying selection in all four derived genetic clades. Indeed, the Δ frequency was significantly smaller than 0, especially for older retrotransposons, meaning that old retrotransposons are less common than neutrally evolving SNPs at the same age. This further demonstrates that even after accounting for the different genetic clades and using a large sample size, retrotransposons evolve under purifying selection in B. distachyon.

We also revealed that only a minority of retrotransposons evolved neutrally, as the observed shape of the Δ frequency curve could only be reproduced in our simulation if the proportion of neutrally evolving mutations in our focal mutations was below 10%. This estimate gives a first glimpse into the distribution of fitness effects of new TE insertions, a fundamental parameter in genetics that describes the way in which new TE insertions can contribute to evolution and adaptation [56]. Here, we show for the first time that new TE insertions have a less than 10% chance to insert into the genome of B. distachyon in a way that will allow them to evolve neutrally, advocating for a large potential of TEs to create, through their movement, new phenotypic variation on which selection can act on. PCAs based on TE polymorphisms allowed us to recover the population structure of B. distachyon, implying that demographic history and hence neutral processes may indeed partially explain the differences in the TE distribution we observed between genetic clades, as shown in Arabidopsis thaliana and Arabidopsis lyrata [26, 57], Drosophila melanogaster [58], humans [59] and the green anole (Anolis carolinensis; [30]). However, and in line with our simulations, the first two axes of the PCA explain less than 7% of the variance, indicating that neutrally evolving TEs contribute only mildly to overall TE diversity in our system.

Because TEs can cause phenotypic variation through new insertions [1–8], it is not surprising that most new insertions interfere with the function of the genome, especially in a species with a small genome, such as B. distachyon (272 Mb) [37]. The proportion of neutrally evolving TE polymorphisms is expected to be very small in genes, as insertions in genic regions are likely to result in loss-of-function [1, 2]. Similarly, TE insertions in close proximity to genes are expected to be highly disruptive, as regulatory elements such as cis-regulatory elements are predominantly located in the proximity of genes. In A. thaliana, for instance, TEs located in the vicinity of genes (less than 2 kb) globally result in downregulation [60]. Although only specific families alter gene expression in B. distachyon [41], the observed Λ1 frequency for retrotransposon polymorphisms in genes and in their 1 kb surroundings matched our expectations. The fact that TE polymorphisms located more than 5 kb away from genes are also evolving under purifying selection was more surprising. That said, little is known about the distance between cis-regulatory sequences and genes in B. distachyon. In plants, TEs are believed to affect gene expression in trans through the production of small-interfering RNA [61–65].

Hence, the fact that only a small proportion of TEs can accumulate neutrally indicates that, in a gene-dense genome such as that of B. distachyon (42.5% of the genome are genes) [65], TE insertions in any genomic compartment may result in some cis– or trans-regulatory effects visible to selection.

To further ascertain the strength of purifying selection, we used forward simulation and showed that simulations assuming a moderately weak selection pressure (S = –5 or S = –8) against TE polymorphisms best fitted our observed data. In theory, no TE polymorphisms under strong purifying selection should be present in a natural population, as such mutations are expected to be quickly lost, especially in a predominantly selfing species where most loci are expected to be homozygous. Therefore, it is not surprising that TE polymorphisms in B. distachyon are under weak to moderate selection, as also shown, for example, for the L1 retrotransposons in humans [27] or the BS retrotransposon family in Drosophila melanogaster [58].

While some of the parameters we chose for our simulations, such as the dominance or selfing rate, can affect the efficiency of TE purging, it is unlikely that discrepancies in the true and assumed values for these parameters would have led to drastically different results. For example, we assumed codominance for all mutations, which might not hold true for each TE polymorphism. However, because of the high selfing rate observed in B. distachyon [38], heterozygous loci are expected to be rare, and dominance is unlikely to have a strong impact on our observations. Similarly, with a higher selfing rate, deleterious TE polymorphisms should be removed more efficiently by purifying selection. To check whether a lower selfing rate could allow a higher proportion of TE polymorphisms to evolve neutrally, we reran the simulations assuming fully outcrossing individuals. This also resulted in simulation with weak to moderate

selection strength on TE polymorphisms best fitting the observed data, further strengthening our results.

While the analyses of positive selection and GEA were based on both DNA-transposons and retrotransposons, we only used retrotransposons to assess the strength of selection on TE polymorphisms, as the age-adjusted SFS was developed with the assumption of no back mutations [48]. Yet, DNA-transposons do not solely transpose through cut and paste mechanisms as they would otherwise not be so abundant in Eukaryotic genomes. DNA-transposons can also create extra copies of themselves by transposing during chromosome replication or repair from a position that has already been replicated, or repaired [66]. We therefore repeated the age-adjusted SFS analyses using DNA-transposons to evaluate whether DNA-transposons were affected by similar selective constraints. The folded SFS of DNA-transposons and retrotransposons display similar shifts toward high proportions of rare alleles and Δ frequency deviations from 0 in the age-adjusted SFS of DNA-transposons and retrotransposons are comparable. Hence, we argue that the conclusion drawn for retrotransposons also holds for DNA-transposons, and that purifying selection affect TEs broadly.

Conclusion

Adaptation to different environmental conditions is a complex process that involves various mutation types. Here, we show that the vast majority of TE polymorphisms are under purifying selection in the small genome of B. distachyon. Conversely, only a very small proportion of TEs seem to have contributed to adaptation. The observed lack of neutrally evolving TE polymorphisms in B. distachyon advocates for a large potential of TE polymorphisms to contribute to the genetic diversity and phenotypic variation on which selection can act and highlights the need to consider TE polymorphisms during evolutionary studies. Finally, our work shows that the ability of TEs to cause phenotypic variation does not necessarily translate into being favored during evolution and adaptation over other mutations with more subtle effects, such as SNPs.

Materials and methods

Whole-genome resequencing data

In this study, we analyzed a total of 326 publicly available whole-genome sequencing data from Brachypodium distachyon accessions sampled around the Mediterranean Basin (Fig. 1A; Additional file 2: Table S3). Our B. distachyon dataset consisted of 47 samples published by Gordon et al. [8], 57 samples published by Skalska et al. [67], 65 samples published by Gordon et al. [68], 86 samples published by Stritt et al. [38] and 71 samples published by Minadakis et al. [39], covering all five genetic clades previously described in this species [38, 39]. Each sample was assigned to a genetic clade based on previously published results [39].

Data processing

Raw reads were trimmed using Trimmomatic 0.36 [69] and mapped to the B. distachyon reference genome version 3.0 [37] using bowtie2 [70] and yaha [71], and TE polymorphisms were identified using the TEPID pipeline [72] and the recently updated TE annotation by Stritt et al. [73] and Wyler et al. [65]. TE polymorphisms include both TE insertion polymorphisms (TIPs; insertions absent from the reference genome but present in at least one natural accession) and TE absence polymorphisms (TAPs; insertions present in the reference genome but absent from at least one natural accession). The class, superfamily and family of each TE call were assigned based on the TEPID results and the TE annotation from the reference genome. TIPs that were less than 100 base pairs (bp) apart in different samples and assigned to the same TE family were merged.

Single nucleotide polymorphisms (SNPs) were called using GATK v.4.0.2.1 [74] using HaplotypeCaller [75] following Minadakis et al. [39]. The SNP calls were hard filtered using the following conditions: QD < 5.0; FS > 20.0; SOR > 3.0; MQ < 50.0; MQRankSum < 2.5; MQRankSum > –2.5; ReadPosRankSum < 2.0; ReadPosRankSum > –2.0. Because B. distachyon displays a high selfing rate [38], most genetic variants are expected to be homozygous within an individual. Hence, all TE calls were treated as homozygous, and heterozygous SNP calls were removed from our dataset to reduce false variant calls. Additionally, all sites with multiallelic TE and SNP calls were removed. SNPs were classified as synonymous, non-synonymous and of high fitness effect using SnpEff [76]. SNPs and TE polymorphisms were merged into a single vcf file using custom scripts provided in github (see section Availability of data and materials).

To estimate the age of each SNP and TE polymorphism, the SNPs and TEs found in the A_East, A_Italia, B_East and B_West clades were polarized using the C clade, which was identified as the most ancestral B. distachyon clade [38] and used as the outgroup throughout this study. An estimate for the time of origin of all SNPs and TE polymorphisms was calculated with GEVA, a nonparametric approach that relies on pairwise differences in identity by descent (IBD) regions around the focal mutation to estimate the time of origin [77]. GEVA was run separately for each clade using the genetic map produced by Huo et al. [78] and a mutation rate of 7×10^-9 substitutions/generation. The theoretical prediction of the correlation between allele age and allele frequency of neutrally evolving mutations based on N_e [42] was compared to the observed correlation between allele age and frequency of synonymous SNPs to check the sanity of the age estimates.

The observed SNP and TE diversity was first examined using a principal component analysis (PCA), and correlations between TE diversity and genetic clades were tested with a mantel test using the ade4 package version 1.7-22 [79] in R version 4.1.2 [80]. The folded site frequency spectrum (SFS) was computed for TE polymorphisms and SNPs using the minor allele frequency in R version 4.1.2 [80]. Finally, the map of the geographical distribution of the used accessions was done in R using the rnaturalearth package 0.3.3 [81].

Analyses of positive selection

Regions of the genome affected by positive selection were identified using the integrated haplotype score (iHS), a measure of the amount of extended haplotype homozygosity along the ancestral allele relative to the derived allele for a given polymorphic site [82]. iHS was calculated using the SNP dataset, and regions displaying longer haplotypes and hence high iHS were identified in R using the rehh package [83, 84]. The threshold to distinguish between regions of high iHS and other regions was selected such that less than 5% of the B. distachyon genome was classified as high iHS regions in each clade (Additional file 2: Table S4). Candidate regions under positive selection were defined as all regions that were found to have high iHS in each clade separately.

A first ANCOVA was used to model the number of fixed TE polymorphisms in each clade found in the candidate region under positive selection based on the following genetic features: total number of TEs, TE superfamily, TE age (split into three categories: young: age < 10,000 generations; intermediate: age between 10,000 generations and 60,000 generations; old: age > 60,000 generations), clade, genomic region (a unique ID for each candidate region under positive selection) and iHS classification of the regions in each clade (high or average). A second ANCOVA was used to model the allele frequency of TE polymorphisms found in the candidate region under positive selection based on the following genetic features: TE superfamily, TE age, clade, genomic region and iHS classification of the regions in each clade. The TE superfamily was included to account for different evolutionary behaviors of TEs from different superfamilies. Age accounted for differences in the fixation rate and frequency distribution between young and old TEs. The clade was included to account for clade-specific differences such as differences in N_e. Finally, a unique ID for each candidate region under positive selection was included to account for region-specific differences such as differences in the recombination rate and GC content. In the end, regions that were found to have a high iHS in some clades were compared to the same regions in the other clades. All ANCOVAs were run in R using the car package [85].

The standardized allele frequency of a mutation across populations (X^tX) values [44] were computed for the combined TE and SNP dataset using Baypass version 2.3 [45, 86]. The X^tX values were used to identify over– and under differentiated TE polymorphisms between clades. A pseudo-observed dataset (POD) of 100,000 SNPs was simulated under the demographic model inferred from the covariance matrix of the SNP dataset. The POD was then used to determine the 97.5% (over-differentiated polymorphisms) and 2.5% (under differentiated polymorphisms) quantiles.

Genome-environment association analyses

We identified TE polymorphisms significantly associated with environmental factors using genome-environment association analyses (GEA) following Minadakis et al. [39]. GEAs were run with GEMMA 0.98.5 [87] using the combined TE and SNP vcf file against the following 32 environmental factors extracted by Minadakis et al. [39]: altitude, aridity from March to June, aridity from November to February, annual mean temperature, mean temperature of warmest quarter, mean temperature of coldest quarter, annual precipitation, precipitation of wettest month, precipitation of driest month, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter, precipitation of coldest quarter, mean diurnal Range, isothermality, temperature seasonality, maximum temperature of warmest month, minimum temperature of coldest month, temperature annual range, mean temperature of wettest quarter, mean temperature of driest quarter, precipitation from March to June, precipitation from November to February, solar radiation from March to June, solar radiation from November to February, mean temperature between March and June, mean temperature between November and February, maximum temperature between March and June, maximum temperature between November and February, minimum temperature between March and June and minimum temperature between November and February. We applied a False Discovery Rate (FDR, [88]) threshold of 5% to control for false positive rates.

Age-adjusted frequency spectra and analyses of purifying selection

Footprints of purifying selection on TE polymorphisms were first evaluated using folded SFS. An age-adjusted site frequency spectrum (age-adjusted SFS) approach was used to further investigate the impact of purifying selection on retrotransposons while accounting for nonconstant transposition rates. Briefly, the age-adjusted SFS is a summary statistic that describes the difference between the average frequency of TEs at a specific age and the average frequency of neutral sites of the same age [48]. Therefore, the TE dataset was sorted by age and split into equally large bins with respect to the number of observations in each age bin. Neutral sites were then randomly down-sampled to match the number of observations in the TE dataset and its age distribution [48].

The difference between the average TE and neutral site frequency, or Δ frequency, was computed for each age bin [48]. This method allows for an unbiased comparison between the allele frequencies of TEs and neutral sites, and is robust to transposition rate changes and demographic changes [48]. However, the theory behind this method was developed assuming no back mutations and is therefore best suited for retrotransposons, as DNA-transposons can exit an insertion site [48]. We used the synonymous SNPs identified with SnpEff as the neutrally evolving sites. However, because estimating the population wide frequency of TEs is more challenging than estimating SNP frequencies, putative biases in frequency estimates need to be assessed before performing age-adjusted SFS analyses. To do so, the SNP dataset was resampled so that the SNP dataset used in the age-adjusted SFS had a frequency distribution that matched the observed TE frequency distribution. The age-adjusted SFS of retrotransposons was contrasted against the age-adjusted SFS of non-synonymous, as well as against high fitness effect SNPs. Therefore, 10,000 non-synonymous and high fitness effect SNPs were randomly selected for each clade to reach approximately the same number of retrotransposon polymorphisms, non-synonymous and high fitness effect SNPs for final comparisons. To estimate the variation in Λ1 frequency estimates, all age-adjusted SFS were computed 100 times. All Wilcoxon tests and Bonferroni p value corrections were done in R version 4.1.2 [80].

Forward simulation

We used SLiM 4.0.1 [89, 90] to run forward simulations and assess the proportion of neutrally evolving retrotransposons and the average selection strength affecting them. The simulations were designed to reflect the population size and demographic history of B. distachyon. The simulated genomic fragment was 1 megabase (Mb) long and included neutral (synonymous) mutations as well as focal mutations that evolved under different selective constraints. The focal mutations were a mix of neutrally evolving mutations and mutations evolving under a constant selection pressure. Therefore, the ratio (r) of focal mutations that evolved neutrally was either 0%, 5%, 10%, 25% or 50%. The scaled selection coefficient (S, defined as N_es, with s the strength of selection and N_e the effective population size) affecting the remaining focal mutations was set at the beginning of the simulation to be either –1, –5, –8, –10, –12, –15, –20 or –50 to cover effectively neutral (0 > S ≥ –1), intermediate (–1 > S ≥ –10) and strongly deleterious (–10 > S) selective constraints. The selfing rate was set to 70%, as B. distachyon is a highly selfing species with occasional outcrossing [38, 39]. In addition, a high recombination rate was chosen to minimize the effects of linked selection in the small genomic fragment simulated. Simulations for each combination of these two parameters were run 20 times to assess the variation in the resulting age-adjusted SFS. The shape of the resulting age-adjusted SFS was used to narrow down the ratio of neutrally evolving TE polymorphisms. Similarly, the age distribution of the mutations in the oldest bin of the age-adjusted SFS was used to narrow down the strength of selection affecting TE polymorphisms.

Acknowledgements

We would like to thank Jeffrey Ross-Ibarra and Mitra Menon as well as Fabrizio Menardo, Michael Thieme, Wenbo Xu, Jigisha, Lars Kaderli and Serafin Schefer for all the discussions and their comments on this project. We thank Emmanuelle Botté for professional editing.

Declarations

Consent for publication

All authors have read and approved the submission of this manuscript.

Availability of data and materials

The datasets supporting the conclusions of this article are publicly available on the European Nucleotide Archive (https://www.ebi.ac.uk/ena/browser/home) and National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/sra/), and the archive numbers of the accessions used are listed in the Additional file 1: Table S2. The scripts generated are available on GitHub https://github.com/Roberthorv/TE_in_Brachypodium/tree/main.

Funding

This project was funded by the Swiss National Science Foundation (project 31003A_182785) and the Research Priority Program Evolution in Action from the University of Zürich. Data analyzed in this paper were generated in collaboration with the Genetic Diversity Center (GDC), ETH Zürich.

Authors’ contributions

A.R. and R.H. conceived the study. R.H. carried out the study, did the TE and SNP calling, performed the age-adjusted SFS analyses, conducted the X^tX analyses and ran the forward simulations. N.M. performed the GEA. Y.B. ran the iHS analyses. A.R. acquired the fundings. R.H. wrote the manuscript and A.R. revised it. All authors discussed the results and commented on the manuscript.