Michael Thieme
email: michael.thieme@botinst.uzh.ch

Nikolaos Minadakis

Christophe Himber

Bettina Keller

Wenbo Xu

Kinga Rutowicz

Calvin Matteoli

Marcel Böhrer

Bart Rymen

Debbie Laudencia-Chingcuanco

John Vogel

Richard Sibout

Christoph Stritt

Todd Blevins

Anne C Roulin
email: anne.roulin@botinst.uzh.ch

Transposition of HOPPLA in siRNA-deficient plants suggests a limited effect of the environment on retrotransposon mobility in Brachypodium distachyon

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Transposable elements (TEs) are DNA sequences with the ability to form extrachromosomal copies and to reintegrate elsewhere into the host genome. In plants, TE-derived sequences are ubiquitous and can constitute more than 80 % of the genome [START_REF] Mirouze | Transposable elements, a treasure trove to decipher epigenetic variation: insights from Arabidopsis and crop epigenomes[END_REF] . In addition to playing a major role in genome size variation (e.g. [START_REF] Hawkins | Differential lineage-specific amplification of transposable elements is responsible for genome size variation in Gossypium[END_REF][START_REF] Piegu | Doubling genome size without polyploidization: dynamics of retrotranspositiondriven genomic expansions in Oryza australiensis, a wild relative of rice[END_REF][START_REF] Wang | Which factors contribute most to genome size variation within angiosperms?[END_REF]Yang et al., 2023) , TEs can alter gene expression by acting as promotors or by providing cis-regulatory elements to flanking regions [START_REF] Butelli | Retrotransposons control fruit-specific, cold-dependent accumulation of anthocyanins in blood oranges[END_REF][START_REF] Makarevitch | Transposable elements contribute to activation of maize genes in response to abiotic stress[END_REF][START_REF] Roquis | Non-coding RNA polymerases that silence transposable elements and reprogram gene expression in plants[END_REF] . TEs are therefore a major source of genetic change. Given that they are more likely than classic point mutations to cause extreme changes in gene expression and phenotypes [START_REF] Uzunović | Transposable Elements Are Important Contributors to Standing Variation in Gene Expression in Capsella Grandiflora[END_REF][START_REF] Thieme | Experimentally heatinduced transposition increases drought tolerance in Arabidopsis thaliana[END_REF][START_REF] Raúl | Transposons are important contributors to gene expression variability under selection in rice populations[END_REF] , they might be especially useful when the survival of an organism or its descendants depends on a quick response to new or challenging environmental conditions (for review [START_REF] Kawakatsu | Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions[END_REF][START_REF] Rey | Adaptation to Global Change: A Transposable Element-Epigenetics Perspective[END_REF][START_REF] Lanciano | Transposable elements: all mobile, all different, some stress responsive, some adaptive?[END_REF][START_REF] Dubin | Transposons: a blessing curse[END_REF] . Paradoxically, while population genomics data have revealed ongoing TE activity in natural plant populations (e.g. [START_REF] Stuart | Population scale mapping of transposable element diversity reveals links to gene regulation and epigenomic variation[END_REF][START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF][START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF] , only a handful of TE families have been reported to transpose in real-time. Therefore, how often or under which natural conditions TEs are activated in the wild remain open questions. In addition, while ongoing transposition is essential for TEs to survive, the presence of mobile and potentially 'selfish' DNA sequences requires the host to evolve robust silencing mechanisms to prevent an uncontrolled TE proliferation. TE activity thus remains a major puzzle in the field of evolutionary genomics.

In plants, the defence against TEs is multi-layered, comprising repressive histone modifications, DNA methylation and RNA interference [START_REF] Saze | DNA Methylation in Plants: Relationship to Small RNAs and Histone Modifications, and Functions in Transposon Inactivation[END_REF][START_REF] Zhang | Dynamics and function of DNA methylation in plants[END_REF][START_REF] Liu | The Epigenetic Control of the Transposable Element Life Cycle in Plant Genomes and Beyond[END_REF] . One of the main players of TE silencing is the RNA-directed DNA methylation (RdDM) pathway, which involves two plant specific RNA-polymerases derived from Pol II, namely Pol IV and Pol V. The largest subunits of each polymerase (NRPB1, NRPD1 and NRPE1, respectively) assemble with other proteins into enzymes with distinct RNA products and functions [START_REF] Ream | Subunit compositions of the RNA-silencing enzymes Pol IV and Pol V reveal their origins as specialized forms of RNA polymerase II[END_REF]for review Rymen et al., 2020) . As a core component of RdDM, Pol IV (including NRPD1) transcribes TE regions into the precursors of functionally specialized 24 nt small interfering RNAs (siRNAs) [START_REF] Sigman | The First Rule of Plant Transposable Element Silencing: Location, Location, Location[END_REF]Rymen et al., 2020;[START_REF] Liu | The Epigenetic Control of the Transposable Element Life Cycle in Plant Genomes and Beyond[END_REF] . Upon the base pairing of 24 nt siRNAs to scaffold transcripts produced by Pol V, the de novo DNA methyltransferase DRM2 [START_REF] Zhong | Molecular Mechanism of Action of Plant DRM De Novo DNA Methyltransferases[END_REF] is recruited to mediate the methylation and subsequent transcriptional repression of TEs. The essential role of RdDM in TE silencing has been shown in A. thaliana, where the knockout of Pol IV and resulting depletion of 24 nt siRNAs leads to a drastically increased heat-dependent transposition of the ONSEN family [START_REF] Tittel-Elmer | Stress-Induced Activation of Heterochromatic Transcription[END_REF][START_REF] Ito | An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress[END_REF] .

The case of the heat-responsive ONSEN elements not only illustrates the importance of epigenetic silencing in regulating TEs but also demonstrates that environmental factors may modulate the dynamics of TEs in plants. Since their discovery by Barbara McClintock, who linked the mobility of AC/DS elements in maize to the occurrence of a 'genomic shock' [START_REF] Mcclintock | The significance of responses of the genome to challenge[END_REF] , the activity of TEs has been frequently associated with the presence of biotic and abiotic stressors. In fact, certain TEs can sense specific physiological states of their host and use them to initiate their own life cycle [START_REF] Negi | Moving through the stressed genome: Emerging regulatory roles for transposons in plant stress response[END_REF] . Besides ONSEN in A. thaliana [START_REF] Cavrak | How a Retrotransposon Exploits the Plant's Heat Stress Response for Its Activation[END_REF] , the cold inducible Tcs1 element in blood oranges [START_REF] Butelli | Retrotransposons control fruit-specific, cold-dependent accumulation of anthocyanins in blood oranges[END_REF] or Tos17 that gets activated during tissue culture in rice [START_REF] Hirochika | Retrotransposons of rice involved in mutations induced by tissue culture[END_REF] are two other prominent cases of stress responsive TEs in plants. Mechanistically, stress can activate TEs via specific motifs allowing the binding of transcription factors and the subsequent recruitment of the transcription machinery to their promoter-like sequences [START_REF] Cavrak | How a Retrotransposon Exploits the Plant's Heat Stress Response for Its Activation[END_REF][START_REF] Grandbastien | LTR retrotransposons, handy hitchhikers of plant regulation and stress response[END_REF][START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF][START_REF] Zhang | Transposable elements orchestrate subgenome-convergent and -divergent transcription in common wheat[END_REF] . The small window of increased activity during well-defined physiological states suggests that some TEs have evolved a distinct lifestyle or 'niche' to successfully reproduce [START_REF] Kidwell | Transposable elements as sources of variation in animals and plants[END_REF][START_REF] Venner | Dynamics of transposable elements: towards a community ecology of the genome[END_REF][START_REF] Stritt | Rare transposable elements challenge the prevailing view of transposition dynamics in plants[END_REF][START_REF] Stitzer | The genomic ecosystem of transposable elements in maize[END_REF] .

Since transcription constitutes the initial step to transposition, the abundance of TE transcripts is often used as a proxy for TE activity [START_REF] Lanciano | Measuring and interpreting transposable element expression[END_REF] . However, the life cycle of TEs is complex [START_REF] Schulman | Retrotransposon replication in plants[END_REF] and the fate of TE transcripts depends on many factors. For instance, several transcriptionally active TEs have accumulated mutations that prevent the production of enzymes needed for their autonomous transposition [START_REF] Tanskanen | Life without GAG: the BARE-2 retrotransposon as a parasite's parasite[END_REF] . To selectively capture TEs that are not only transcriptionally active but also capable of transposing, several protocols have been developed, including ALE-seq [START_REF] Cho | Sensitive detection of pre-integration intermediates of long terminal repeat retrotransposons in crop plants[END_REF] , VLP-seq [START_REF] Lee | Arabidopsis retrotransposon virus-like particles and their regulation by epigenetically activated small RNA[END_REF] and mobilome-seq [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] . These recent approaches have been particularly successful when targeting long terminal repeat retrotransposons (LTR-RTs), which represent the largest fraction of TE-derived sequences in plant genomes [START_REF] Vitte | The bright side of transposons in crop evolution[END_REF] . Indeed, LTR-RTs transpose through a copy-and-paste mechanism which involves the reverse transcription of a fulllength RNA intermediate [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF][START_REF] Schulman | Retrotransposon replication in plants[END_REF] . As part of their life cycle and presumably through auto-integration, non-homologous and alternative end-joining, active LTR-RTs can also form extrachromosomal circular DNA (eccDNA) intermediates [START_REF] Flavell | Extrachromosomal circular copies of the eukaryotic transposable element copia in cultured Drosophila cells[END_REF][START_REF] Flavell | Role of reverse transcription in the generation of extrachromosomal copia mobile genetic elements[END_REF][START_REF] Garfinkel | Retrotransposon Suicide: Formation of Ty1 Circles and Autointegration via a Central DNA Flap[END_REF][START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF][START_REF] Lee | Arabidopsis retrotransposon virus-like particles and their regulation by epigenetically activated small RNA[END_REF]Yang et al., 2023a) , whose detection by mobilome-seq can be used as a proxy for their mobility [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] . For instance, mobilome-seq has been successfully used to track full-length eccDNA of mobilized autonomous RTs, containing both LTRs (2-LTR circles), in plants such as A. thaliana and rice [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF][START_REF] Thieme | Inhibition of RNA polymerase II allows controlled mobilisation of retrotransposons for plant breeding[END_REF] .

While the activity of LTR-RTs has been extensively studied in the model A. thaliana [START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF] , the interplay between genetic and environmental factors in other wild plant species remains poorly investigated. To clarify these questions, we exploit here the Brachypodium distachyon diversity panel (Minadakis et al., 2023a) and explore LTR-RT activity in a wild monocot.

B. distachyon is a Mediterranean grass with a compact diploid genome of ~272 Mb (International Brachypodium Initiative, 2010;[START_REF] Hasterok | Brachypodium: 20 years as a grass biology model system; the way forward?[END_REF] harboring 40 LTR-RT families [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] that constitute about 30 % of the genome (International Brachypodium Initiative, 2010) . In B. distachyon, LTR-RTs not only evolved varying insertion site preferences [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] but also significantly differ in terms of transposition dynamics [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF][START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . While we previously suggested an ongoing transposition of LTR-RTs based on population genomics data [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF][START_REF] Horvath | The evolution of transposable elements in Brachypodium distachyon is governed by purifying selection, while neutral and adaptive processes play a minor role[END_REF], here we aimed to clarify the contribution of genetic and environmental factors to LTR-RT activity in B. distachyon. To that end, we combined population genomics data available for 320 natural accessions with mobilome-seq under different stress conditions and asked: (i) does the accumulation of LTR-RTs in these natural accessions correlate with environmental variables, (ii) is LTR-RT mobility induced by specific stresses, and (iii) which genetic factors influence the accumulation of LTR-RTs ?

Results

Abundances of LTR-RT families differ but show limited association with bioclimatic variables

B. distachyon naturally occurs around the Mediterranean rim (Fig 1A) and groups into three main genetic lineages (A, B and C) that further split into five genetic clades [START_REF] Stritt | Migration without interbreeding: Evolutionary history of a highly selfing Mediterranean grass inferred from whole genomes[END_REF]Minadakis et al., 2023a): an early diverged C clade and four clades found in Spain and France (B_West), Italy (A_Italia), the Balkans and coastal Turkey (A_East) and inland Turkey and Lesser Caucasus (B_East).

We first aligned genomic reads of 320 accessions to the LTR-RT consensus sequences of B. distachyon obtained from the TRansposable Elements Platform (TREP). We then computed the abundance of TE-derived sequences, hence a proxy for copy number variation (pCNV), for the 40 annotated LTR-RT families using a coverage-based approach accounting for sample sequencing depth (see Materials and Methods). We favored this approach over an analysis based on transposon insertion polymorphisms (TIPs) because estimates based on TIPs are reference genome-dependent and biased by the phylogeny in our study system. We have, for instance, previously shown that accessions from the B_East clade harbor significantly less TIPs than accessions from the A_East clade due to the fact that the reference genome Bd21 belongs the B_East clade [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF] . In addition, whole-genome de novo assembly of 54 B. distachyon natural accessions and the subsequent pangenome analysis revealed that non-reference accessions display large genomic variations [START_REF] Gordon | Extensive gene content variation in the Brachypodium distachyon pan-genome correlates with population structure[END_REF] , which may further bias the estimates of TIP abundances.

The heatmap produced based on pCNV (Fig 1B) showed that LTR-RTs underwent different transposon accumulations. We found that Gypsy elements (RLG) harbor a higher pCNV than Copia To test whether the accumulation of LTR-RT sequences correlated with environmental factors, we retrieved bioclimatic variables comprising precipitation, temperature, aridity levels, solar radiation and atmospheric pressure at each locality. We then ran linear mixed models (LMM) where pCNV per LTR-RT family was entered as the response variable, the bioclimatic variables entered separately as fixed factors and the clade of origin as random factors to account for population structure. Marginal R 2 extracted for each LMM did not exceed 10% even for the putatively most recently active LTR-RT families (Fig 1E), indicating that albeit significant, the association between pCNV per family and the environment was mild in our study system. We observed similar associations between pCNVs and bioclimatic variables when not accounting for population structure

Recently active LTR-RT families produce eccDNAs

The coverage-based approach for estimating LTR-RT CNVs considers all TE-derived sequences and does not take into account that individual families differ in their age structures, turnover times and proportion of full-length, potentially autonomous mobile elements [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . We therefore complemented our in silico analysis by experimentally testing whether LTR-RT families were indeed still active in planta and whether the mild but significant correlation we observed between global pCNVs and bioclimatic variables may be due to a stress-specific activity.

To cover a wide range of the genetic, geographic and bioclimatic diversity of B. distachyon, we selected nine natural accessions belonging to the five genetic clades and originating from contrasting habitats (Fig 1A). Considering the role of Pol IV in the silencing of TEs in plants, we also included two independent Pol IV mutant lines, a sodium azide mutagenized line (hereafter Bd nrpd1-1) and a T-DNA line (hereafter Bd nrpd1-2), both carrying a homozygous mutation in the largest subunit of Pol IV (NRPD1; Bradi2g34876) in the Bd21-3 background. We exposed plants to eight different stresses, namely cold, drought, heat, salt, submergence, infection with Magnaporthe oryzae, treatment with glyphosate and chemical de-methylation and performed mobilome-seq on all resulting 105 samples (see Material and Methods).

Following the removal of organelle-derived reads, we first assembled mobilome reads and aligned the resulting ten longest contigs of each sample to the reference assembly of B. distachyon.

We then screened genomic regions for which at least three assembled mobilome contigs of different samples aligned and further assessed in which genotypes or stresses those contigs occurred. We only retained circle-forming regions with a specificity above 50 % (i.e., regions for which more than half of the contigs belonged to a certain stress or genotype). In addition, we also kept recurrently active regions, present at a high frequency independently of the stress or the genotype in at least ten samples. In total, we retained 15 circle-forming regions, all of which contained TE sequences (Fig 2). Eight of these corresponded to the Angela family (RLC_BdisC024, RLC_BdisC022, RLC_BdisC030, RLC_BdisC209), four to CRM elements (RLG_BdisC039, RLG_BdisC102), and one each to the SIRE (RLC_BdisC026), the Alesia (RLC_BdisC010) and the non-autonomous and unclassified RLG_BdisC152 family. We hereafter refer to RLC_BdisC024 as HOPPLA (German allusion for the the surprising finding of a jumping element). We did not find stress specificity in the formation of eccDNAs (Fig 2A). However, our results pointed to a genotype-dependent formation of eccDNAs for the RLC_BdisC209, RLC_BdisC026 and HOPPLA families (Fig 2B). In particular, two contigs containing HOPPLA elements were exclusively detected in the two pol IV mutants (Fig 2B).

HOPPLA activity is increased in the pol IV mutants regardless of the stress applied

Fragmented eccDNAs or circles containing only of the two LTRs (1-LTR circles) can be formed following reverse transcription by auto-integration, alternative end-joining in the virus-like particles [START_REF] Garfinkel | Retrotransposon Suicide: Formation of Ty1 Circles and Autointegration via a Central DNA Flap[END_REF]Yang et al., 2023a) or by a recombination of the two LTRs of genomic copies [START_REF] Smith | Small polydisperse circular DNA of HeLa cells[END_REF][START_REF] Gaubatz | Extrachromosomal circular DNAs and genomic sequence plasticity in eukaryotic cells *[END_REF] . Hence, 1-LTR circles do not necessarily imply LTR-RT mobility. In contrast, recent work indicates that 2-LTR circles are formed following the complete reverse transcription by non-homologous end-joining of an intact full-length linear RT copy that is capable of integrating into the genome (Yang et al., 2023a) . Indeed, the detection of full-length 2-LTR circles of well characterized autonomous LTR-RTs such as EVD (ATCOPIA93) has been directly linked to their actual transposition [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] . As a complement to the assembly-based analysis, we thus followed a stringent approach to analyse our mobilome-seq data. We aligned reads to a library comprising artificial 3'LTR-5'LTR fusions of all full-length LTR-RTs annotated in the B. distachyon reference assembly [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . Because such reads are not expected to be present in genomic DNA, this allowed us to speciflangically detect intact 2-LTR circles of extrachromosomal LTR-RTs capable of integrating into the genome. To control for possible traces of undigested genomic DNA that may subsequently be amplified by the Phi29 enzyme during mobilome-seq [START_REF] Silander | Whole genome amplification with Phi29 DNA polymerase to enable genetic or genomic analysis of samples of low DNA yield[END_REF] , we also included publicly available genomic reads of all nine accessions in our analysis.

We found that several LTR-RTs formed eccDNA with 2-LTR junctions. Yet, most of them occurred sporadically and we did not observe a recurring stress-specific formation of 2-LTR circles for any of the 37 LTR-RT families with annotated full-length copies (Fig 3A). For instance, we found a very strong signal for RLC_BdisC031 that was solely detected in glyphosate-treated Bd21-3 plants and therefore not further considered in the analysis. In contrast, and in accordance with the assembly-based approach, the Angela element HOPPLA showed a recurrent formation of 2-LTR circles. However, this formation was not triggered by a specific stress and only occurred in the two independent pol IV mutants (Fig 3B). Finally, our attempt to transiently inhibit LTR-RT silencing using alpha-amanitin and Zebularine (a combination of inhibitors shown to increase the activity of LTR-RTs in A. thaliana and rice), did not result in the consistent activation of HOPPLA or other LTR-RTs in multiple accessions (Fig 3A). The presence of 2-LTR circles of HOPPLA in Bd nrpd1-2 (-/-) was confirmed by an inverse PCR on total DNA that was not subjected to a rolling circle amplification, with outward facing primers specific to the two LTRs (Fig 3C). Notably, we also detected a faint signal for the outcrossed line Bd NRPD1 (+/+) suggesting a weak activity of HOPPLA in wild-type plants.

Since the stress-or tissue-dependent activity of LTR-RTs is mediated by the specific binding 230 of transcription factor (TF), we screened the consensus sequence of HOPPLA for motifs of known TF binding sites. First, we validated this approach by analyzing one of most active copies (AT1G11265) 232 of the heat-responsive ONSEN family of A. thaliana. As expected, a GO term analysis indicated a strong enrichment of heat-responsive TFs for this element (S2 Fig) . In contrast to the well-known, stress-responsive ONSEN LTR-RT, the GO terms of TFs that could bind to HOPPLA indicated that developmental processes and auxin-activated signaling pathways played a role in its activity, rather than specific stresses (Fig 3D).

Members of the HOPPLA family differ in activity

Because individual copies of the same LTR-RT family can differ in their activity [START_REF] Cavrak | How a Retrotransposon Exploits the Plant's Heat Stress Response for Its Activation[END_REF] , we also anad the relative abundance of 2-LTR-spanning reads for each annotated full-length copy of HOPPLA. This analysis revealed a great diversity of eccDNA formation among individual copies of the HOPPLA family and confirmed the strongest activity of HOPPLA in the two pol IV mutants (Fig 4A).

We also obtained a few reads spanning the 2-LTR-junction of the two most active HOPPLA copies (Bd3_22992889 and Bd4_25471847) in the Bd NRPD1 (+/+) control line (Fig 4A), which confirmed the weak but detectable band for the inverse 2-LTR PCR for Bd NRPD1 (+/+) (Fig 3C).

Using the meta information of individual HOPPLA copies described previously [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] , we further assessed which genomic factors (DNA methylation, CG content, distance to the closest gene or age of the copy) were linked to the 2-LTR circle formation for individual HOPPLA copies. While no clear pattern emerged from this analysis, the most active copies of HOPPLA tend to be rather young (Fig 4B).

HOPPLA is targeted by Pol IV-dependent 24 nt siRNAs in the wild type and transposes in pol IV

mutant plants

The pivotal role of Pol IV in producing TE-specific 24-nt siRNAs for RNA-directed DNA methylation has been demonstrated in many plant species including A. thaliana [START_REF] Ito | An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress[END_REF] , rice [START_REF] Xu | Regulation of Rice Tillering by RNA-Directed DNA Methylation at Miniature Inverted-Repeat Transposable Elements[END_REF] and for the Alesia family (RLC_BdisC010) in B. distachyon [START_REF] Böhrer | Integrated Genome-Scale Analysis and Northern Blot Detection of Retrotransposon siRNAs Across Plant Species[END_REF] . To confirm that the increased production of 2-LTR eccDNA circles of HOPPLA in mutants deficient in B. distachyon NRPD1 is correlated with a depletion of 24 nt siRNAs, we performed a small RNA blot including To complete their life cycle, reverse transcribed extrachromosomal copies of LTR-RTs have to 270 integrate into the host genome [START_REF] Schulman | Retrotransposon replication in plants[END_REF] . Because all our analyses congruently pointed to the activity of HOPPLA in the Bd nrpd1-2 (-/-) mutant, we sequenced the genome of seven 272

Genome-wide association studies for pCNVs do not recover known components of RdDM

To decipher the genetic basis of HOPPLA accumulations in natural populations, we first performed a genome-wide association study (GWAS) using HOPPLA pCNVs in the diversity panel of 320 natural accessions (Fig 1B) as a phenotype. We identified only one region with two significant peaks (FDRadjusted p-value < 0.05, Bd5 6920000-6960000 and 7210000-7240000) obtained for the mobile HOPPLA family (Fig 6). Because inserted copies of HOPPLA may themselves lead to significantly associated regions in the GWAS as shown in A. thaliana [START_REF] Quadrana | The Arabidopsis thaliana mobilome and its impact at the species level[END_REF] , we first verified that there were neither TIPs [START_REF] Horvath | The evolution of transposable elements in Brachypodium distachyon is governed by purifying selection, while neutral and adaptive processes play a minor role[END_REF] nor annotated reference insertions of HOPPLA in that region (Fig 6). As described above, our data suggested that the loss of 24-nt siRNAs in the Pol IV mutants was sufficient to mobilize HOPPLA. We therefore further tested whether any of the genes encoding subunits of Pol IV or Pol V would be localized in or near this region (window size 50 kb) (S3 Table ). We did not detect any known Pol IV or Pol V-related genes, but instead found 296

Bradi5g05225, an ortholog of the A. thaliana ROS1-associated methyl-DNA binding protein 1 (RMB1, AT1g63240) (Liu et al., 2021a) , to be co-localized with the peak. 298

To test whether genomic regions might be recurrently associated with their pCNVs, we finally extended the GWAS analyses to the four other most recently active families (RLC_BdisC022, 300 RLC_BdisC030, RLC _BdisC209 and RLG_BdisC004) (S6 Fig) . We extracted the candidate genes for andMaterial andMethods). Apart from the two closely related families RLC_BdisC030 and RLC_BdisC209 that shared the majority of their GWAS candidates, and RLC_BdisC022 and RLC_BdisC209 that shared three genes, we found no overlap of annotated loci potentially contributing to the pCNVs of recently active families (Fig 6).

Discussion

Understanding the dynamics of TEs and their role in adaptation is currently one of the major challenges in the field of evolutionary genomics. The fact that mobile TEs are a source of epi/genetic diversity and potential drivers of evolution has been demonstrated in many organisms including fungi [START_REF] Muszewska | Transposable elements contribute to fungal genes and impact fungal lifestyle[END_REF] , insects [START_REF] Gilbert | Transposable Elements and the Evolution of Insects[END_REF] , mammals [START_REF] Senft | Transposable elements shape the evolution of mammalian development[END_REF] and plants [START_REF] Lisch | How important are transposons for plant evolution?[END_REF] . However, while there are a number of examples showing that certain TE insertions facilitated the adaptation to changing environments (for review [START_REF] Bourgeois | On the Population Dynamics of Junk: A Review on the Population Genomics of Transposable Elements[END_REF] , TEs are generally harmful [START_REF] Bourgeois | Disentangling the determinants of transposable elements dynamics in vertebrate genomes using empirical evidences and simulations[END_REF][START_REF] Horvath | The evolution of transposable elements in Brachypodium distachyon is governed by purifying selection, while neutral and adaptive processes play a minor role[END_REF][START_REF] Langmüller | The genomic distribution of transposable elements is driven by spatially variable purifying selection[END_REF] and therefore controlled by complex silencing mechanisms. To foster our understanding of TE activity, we investigated the environmental conditions and genetic factors associated with the accumulation and mobility of LTR-RTs in plant genomes. By measuring LTR-RT pCNVs in a panel of 320 B. distachyon natural accessions, we show that the intra-specific variations of pCNVs of RLC elements, but not the pCNVs of the generally older and more abundant RLG elements [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] , separate accessions according to their genetic cluster of origin. This is even more striking for members of the Angela (RLC_BdisC022, RLC_BdisC024, RLC_BdisC030, RLC_BdisC209) and the Tekay (RLG_BdisC004) family, which are the youngest families in B. distachyon [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . Highly polymorphic among natural accessions of B. distachyon [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF][START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] , they are expectedly the main drivers of lineage-specific expansions of pCNVs in our study system. Hence, we do not only confirm that LTR-RT families in B. distachyon globally differ in size [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] but also demonstrate that the accumulation of genomic sequences derived from specific families varies significantly among natural accessions.

The transcriptional activity of LTR-RTs can be triggered by specific environmental stresses [START_REF] Grandbastien | LTR retrotransposons, handy hitchhikers of plant regulation and stress response[END_REF][START_REF] Negi | Moving through the stressed genome: Emerging regulatory roles for transposons in plant stress response[END_REF] . Given that B. distachyon occurs in a wide range of different habitats in the Mediterranean area [START_REF] Minadakis | The demographic history of the wild crop relative Brachypodium distachyon is shaped by distinct past and present ecological niches[END_REF] , this characteristic feature of LTR-RTs provides a potential explanation for the pCNVs we observed across natural accessions [START_REF] Quadrana | The Arabidopsis thaliana mobilome and its impact at the species level[END_REF][START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF] . Yet, for all LTR-RT families, pCNVs correlate only moderately with environmental factors. Consequently, our genomic data do not support a large effect of the environment on LTR-RT activity in B. distachyon. While this result could seem startling, it is not completely surprising. Indeed, many LTR-RT families, and especially the old RLG elements, do not show signs of increased activity in the recent past in B. distachyon [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . Considering that their copy number expansions took place in a climate that has drastically changed following the last glacial maximum [START_REF] Minadakis | The demographic history of the wild crop relative Brachypodium distachyon is shaped by distinct past and present ecological niches[END_REF] , a limited link between their activity and the current environmental conditions is actually expected for most families. In contrast, the lack of correlation between current bioclimatic variables and copy number variation for families with ongoing activity (RLC_BdisC022, RLC_BdisC024, RLC_BdisC030, RLC_BdisC209 and RLG_BdisC004; [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] , suggests a more complex mechanism than their pure dependence on specific stresses in certain environments. This hypothesis is supported by previous findings in A. thaliana. Indeed, a minor impact of the environment on transpositional activity was also found in this species, where the two most associated environmental variables ('seasonality of precipitation' and 'diurnal temperature range') only explained about 9 % of the observed variation [START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF] .

Genetic factors are well-known to be essential in regulating LTR-RT activity [START_REF] Miura | Mobilization of transposons by a mutation abolishing full DNA methylation in Arabidopsis[END_REF][START_REF] Tsukahara | Bursts of retrotransposition reproduced in Arabidopsis[END_REF][START_REF] Mirouze | Selective epigenetic control of retrotransposition in Arabidopsis[END_REF][START_REF] Bourque | Ten things you should know about transposable elements[END_REF] . As the loss of main players of the RdDM silencing pathway leads to increased TE activity [START_REF] Ito | An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress[END_REF][START_REF] Benoit | Environmental and epigenetic regulation of Rider retrotransposons in tomato[END_REF][START_REF] Xu | Regulation of Rice Tillering by RNA-Directed DNA Methylation at Miniature Inverted-Repeat Transposable Elements[END_REF][START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF] , the two B. distachyon Pol IV (NRPD1) mutants provided an ideal functional tool to experimentally validate our in silico analysis. Since transcriptionally active LTR-RTs are not necessarily able to transpose [START_REF] Bajus | A review of strategies used to identify transposition events in plant genomes[END_REF] , we used a mobilome-seq approach to detect TEderived eccDNAs and transpositionally active LTR-RT families. We deliberately followed a very stringent approach for analysing the data and by doing so, identified HOPPLA as the only highly active LTR-RT family in B. distachyon. Indeed, HOPPLA is the only family for which we further detect newly inserted copies in the Pol IV mutant.

The non-stress-specific activity of HOPPLA in the two independent Pol IV mutants supported our in silico approach and strengthened the idea that genetic, rather than environmental stresses, are major drivers of LTR-RT activity in B. distachyon. These results are also in line with our TF binding sites analysis. In contrast to the heat-responsive A. thaliana element ONSEN, for which we predominantly recovered TF-binding sites associated with heat response, HOPPLA seemed to be targeted by TFs involved in developmental processes and auxin signaling [START_REF] Leyser | Auxin Signaling[END_REF] . Follow-up studies should therefore address the question of whether the activity of HOPPLA or other families differs between tissues or developmental stages, as observed for the endosperm-specific mobility of PopRice in rice, for example [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] . Strikingly, despite the central role of Pol IV in the RdDM pathway, we did not observe bursts of multiple LTR-RT families but instead found that the loss of 24 nt siRNAs specifically activated individual copies of the HOPPLA family. Interestingly, we also detected a weak signal for 2-LTR eccDNAs in the Bd21-3 wt and the outcrossed line Bd NRPD1 (+/+) but not in other natural accessions. This suggests that the accession-specific composition of the mobilome, and hence the genetic background of the pol IV mutant line, plays an important role in LTR-RT activity. Related to this, we sporadically observed very strong signals for individual samples, which could indicate an accession-specific response of the mobilome to certain triggers.

Given that pCNVs vary greatly among genetic clades, assessing the effect of a genetic mutation of major components of the RdDM pathway in a set of genetically diverse natural accessions would be timely, yet labor-intensive as transformation works more efficiently in the Bd21-3 background than in the other accession tested. Our attempt to transiently reduce LTR-RT silencing in multiple accessions from different genetic clades using the chemical inhibition of Pol II and DNA methyltransferases [START_REF] Thieme | Inhibition of RNA polymerase II allows controlled mobilisation of retrotransposons for plant breeding[END_REF] did not result in an increased activity of HOPPLA or members of any other LTR-RT family. In addition, and despite the differences of activities observed among individual HOPPLA copies, we could not detect, in the present study, a clear link between their activity and GC contents or methylation states. Taken together, these findings suggest that the specific function of the canonical RdDM with Pol IV, rather than generic DNA methylation states are regulating HOPPLA activity. Yet, our GWAS failed to recover major components of the RdDM pathway. Instead, the diversity of activity within the HOPPLA family may suggests that the presence of single active copies could determine the fate of an entire family. In addition, pCNV is also dependent on the removal rate of LTR-RT families, which varies greatly in B. distachyon [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . This complexity of parameters affecting the dynamics of LTR-RT might explain why none of the genes known to be involved in silencing LTR-RTs are associated with pCNV for HOPPLA or any recently active family. Our candidate locus containing Bradi5g05225, a gene related to RMB1 whose loss of function has been shown to result in DNA hypermethylation (Liu et al., 2021b) , remains nonetheless a great candidate for functional validation.

Altogether, our work confirms that LTR-RTs in B. distachyon are 'well-behaved' [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF] and that the evolutionary consequences of their mobility are hard to study in real-time. Indeed, while mobilome-seq revealed a sporadic activity for other families, we only found recurring activity and new insertions of HOPPLA in the pol IV mutant. These results somewhat contrast with our population genomics analyses which clearly indicate an ongoing activity of several LTR-RT families in natural accessions. We propose that the activity of LTR-RTs is relatively low and might depend on a complex interaction between genetic factors, developmental stages and, more marginally, the punctual occurrence of stresses.

Material and methods

Estimation of LTR-RT pCNVs

We used publicly available genomic reads of 320 sequenced natural accessions of B. distachyon [START_REF] Gordon | Extensive gene content variation in the Brachypodium distachyon pan-genome correlates with population structure[END_REF][START_REF] Gordon | Gradual polyploid genome evolution revealed by pan-genomic analysis of Brachypodium hybridum and its diploid progenitors[END_REF][START_REF] Skalska | Genetic and methylome variation in Turkish brachypodium distachyon accessions differentiate two geographically distinct subpopulations[END_REF][START_REF] Stritt | Migration without interbreeding: Evolutionary history of a highly selfing Mediterranean grass inferred from whole genomes[END_REF]Minadakis et al., 2023a) to assess the natural variation of copy numbers of LTR-RTs. Because sequencing depth differed substantially between accessions (Minadakis et al., 2023a), we first downsampled all fastq files to the read number of the sample with the lowest number of reads (4.230.721 reads) using the reformat.sh function of BBtools (v 38.75, BBMap Bushnell B., sourceforge.net/projects/bbmap/). Downsampled reads were aligned to the TREP consensus sequences of LTR-RTs and the reference assembly of Bd21 (v 3.0) using BWA-MEM (v 0.7.17-r1188) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] with the -M and the -a options set, hence outputting all alignments found. Coverages of LTR-RTs and the reference assembly were assessed using bedtools (v 2.30.0) [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] genomecov with the -d and -split parameters set. For each of the LTR-RT families, a proxy for copy number was obtained by normalizing the coverage signals by the coverage of the entire reference assembly and by correcting for the length of the consensus sequences. We favoured this coverage-based approach over an analysis based on transposon insertion polymorphisms (TIPs) as estimates based on TIPs are reference-dependent and biased by the phylogeny in our study system. We have for instance previously shown that accessions from the B_East clade harbor significantly less TIPs than accessions from the A_East clade due to the fact that the reference genome Bd21 belongs the B_East clade [START_REF] Stritt | Recent Activity in Expanding Populations and Purifying Selection Have Shaped Transposable Element Landscapes across Natural Accessions of the Mediterranean Grass Brachypodium distachyon[END_REF] . In addition, whole-genome de novo assembly of 54 B. distachyon natural accessions and the subsequent pangenome analysis revealed that non-reference accessions display large genomic variations [START_REF] Gordon | Extensive gene content variation in the Brachypodium distachyon pan-genome correlates with population structure[END_REF] which may further bias the estimates of TIP abundances. pCNV raw data were processed using R (v 3.6.3 and 4.0.2) (R Core Team, 2020) in Rstudio (RStudio Team, 2016) .

Variation in pCNVs across the 320 natural accessions was visualized with a heatmap drawn with the heatmap() function natively provided in R version 4.0.2. We computed pairwise genetic distances between accessions with the R package pvclust v 2.2.0 [START_REF] Suzuki | Pvclust: an R package for assessing the uncertainty in hierarchical clustering[END_REF] . The resulting tree was used to order accessions phylogenetically on the heatmap. PCAs based on pCNVs were obtained with the R package ggbiplot v 0.55 [START_REF] Vu | ggbiplot: A ggplot2 based biplot[END_REF] .

To test for an association between pCNVs and environmental variables, we retrieved information about climatic variables at each local site from [START_REF] Minadakis | The demographic history of the wild crop relative Brachypodium distachyon is shaped by distinct past and present ecological niches[END_REF] . Linear mixed model analyses where pCNV per LTR-RT family was entered as the response variable, the bioclimatic variables entered separately as fixed factors and the clade of origin as random factors to account for population structure were ran with the R package lme4 [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF] . The part of the variance explained by the fixed-(marginal R 2 ) were computed following [START_REF] Nakagawa | A general and simple method for obtaining R2 from generalized linear mixed-effects models[END_REF] and visualized as bubble plot with the R package ggplot2 [START_REF] Wickham | ggpolt2 Elegant Graphics for Data Analysis[END_REF] . Classical linear models were run in base R.

Plant material, growth conditions and stresses for mobilome-seq

Brachypodium distachyon natural accessions used in this study comprised Bd21, Bd21-3, Cm18, Cb23, ABR2, Bd29-1 BdTR13c, RON2 and Arn1. Because Pol IV is known to play an important role in LTR-RT silencing in plants [START_REF] Stonaker | Diversity of Pol IV Function Is Defined by Mutations at the Maize rmr7 Locus[END_REF][START_REF] Ito | An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress[END_REF][START_REF] Ferrafiat | The NRPD1 N-terminus contains a Pol IV-specific motif that is critical for genome surveillance in Arabidopsis[END_REF] , we also included the sodium azide mutagenized pol IV mutant line NaN74 (Bd nrpd1-1) [START_REF] Dalmais | A TILLING Platform for Functional Genomics in Brachypodium distachyon[END_REF][START_REF] Böhrer | Integrated Genome-Scale Analysis and Northern Blot Detection of Retrotransposon siRNAs Across Plant Species[END_REF] , the T-DNA insertion pol IV mutant line JJJ18557 Nr31 [START_REF] Bragg | Generation and Characterization of the Western Regional Research Center Brachypodium T-DNA Insertional Mutant Collection[END_REF] Bd nrpd1-2 (-/-) and a corresponding sibling, outcrossed control line Bd NRPD1 (+/+) in the background of the natural accession Bd21-3. For in vitro experiments, seeds were soaked for 4 h in tap water and, without damaging the embryo, the lemma was carefully peeled off. Seeds were then surface-sterilized for 30 seconds in 100% ethanol and immediately rinsed three times with sterile tap water. Surface-sterilized seeds were placed with the embryo facing down and at an angle of about 30° towards the side, onto solid ½ MS-medium (2.15 g/L MS basal salt without vitamins (Duchefa Biochemie, Haarlem, NL), 0.5 g/L MES-Monohydrate, 10 g/L sucrose , pH 5.8 (KOH), 0.25 % Phytagel (Sigma-Aldrich, St. Louis, USA) in 'De Wit' culture tubes (Duchefa Biochemie, Haarlem, NL).

Plants were grown at 24 °C (day) / 22 °C (night), 16 h light under controlled conditions in an Aralab 600 growth chamber (Rio de Mouro, PT) for 25 to 29 days until the onset of stresses. For salt stress, seedlings were transplanted to solid ½ MS-medium supplied with 300 mM NaCl and grown for five days at 24/22 °C, 16h light. A solution of sterile-filtrated Glyphosate (Sintagro AG, Härkingen,CH) (20 mM, diluted in water) was applied to leaves using a piece of soaked sterile filter paper and plants were incubated for four days at 24/22 °C, 16h light. Drought stress was induced by uprooting plants from the medium and incubating them for 2:15 h at 24 °C in the light. Before sampling, plants were allowed to recover for two hours on fresh ½ MS-medium. For the infection with Magnaporthe oryzae (rice blast) six isolates (FR13,IK81, with spore concentrations between 130'000-200'000 K spores per isolate per mL sterile water, supplied with 0.2 % Tween 20 were mixed and applied with a cotton swab to plant leaves. Plants were incubated for 24 h in the dark (24/22 °C) and then grown for another three days at 24/22 °C, 16h light. For heat stress, plants were incubated for 8 h at 42 °C. Before sampling, heat-stressed plants were allowed to recover for 16 h at 24/22 °C. Cold stress was induced by incubating plants for 24 h at 2 °C on ice at 16 h light. Prior to sampling, plants were allowed to recover for two hours at 24 °C in the light. For submergence stress, two small holes were drilled just above the growth medium and at the top through the wall of the culture tubes. Tubes were then inverted and submerged upside down for 48 h at 24/22 °C, 16h light using a custom rack in a plastic beaker filled with 2.5 liters of 24 °C tap water.

In this way, it is possible to submerge plant leaves without the medium coming into contact with the water. Chemical de-methylation of DNA was conducted according to [START_REF] Thieme | Inhibition of RNA polymerase II allows controlled mobilisation of retrotransposons for plant breeding[END_REF] by germinating and growing plants for 28 days on ½ MS-medium supplied with a mixture of Zebularine (Sigma-Aldrich, St. Louis, USA) and alpha-amanitin (Sigma-Aldrich, St. Louis, USA). Because the drug treatment severely affected the growth of seedlings, we omitted a treatment of mutant plants and used reduced concentrations of 20 uM (Zebularine) and 2.5 mg/ml (alpha-amanitin), respectively for all natural accessions.

Mobilome sequencing and validation of eccDNAs

DNA was extracted using the DNeasy plant kit (Qiagen, Venlo, Netherlands) according to the protocol of the manufacturer. DNA concentration was measured using the Qubit high sensitivity kit (Invitrogen, Waltham, USA). Mobilome sequencing was performed according to [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] using pooled DNA of two biological replicates per sample. For this, 50 ng of DNA from both biological replicates were pooled and diluted to a volume of 58 μL. To enrich eccDNA, DNA was first purified using the GENECLEAN kit (MP Biomedicals, Santa Ana, USA) according to manufactures recommendations using 5 μL glass milk with an elution volume of 35 μL. Thirty μL of the eluate were digested using the Plasmid-Safe ATP-dependent DNase (Biosearch Technologies, Hoddesdon, UK) for 17 h at 37 °C. The digestion product was then subjected to an ethanolic precipitation and the precipitated eccDNA amplified using the illustraTempliPhi Amplification Kit (Cytiva, Marlborough, USA) according to [START_REF] Lanciano | Sequencing the extrachromosomal circular mobilome reveals retrotransposon activity in plants[END_REF] with an extended incubation time of 65 h at 28 °C. The templiphi product was diluted 1:10, quantified using the Qubit high sensitivity kit and 120 ng per sample were used for library preparation. Sequencing libraries were prepared using the Nextera DNA Flex Library Prep and the Nextera DNA CD Indexes (Illumina, San Diego, USA). Quality of libraries were assessed using the Tape Station (Agilent Technologies, Santa Clara, USA) with High Sensitivity D1000 screen tapes and concentrations were measured using the Qubit high sensitivity kit. Up to 12 indexed libraries were pooled and sequenced with an Illumina MiSeq sequencer using the MiSeq reagent kit v3 (600 cycles). Raw reads have been uploaded to ENA (accession number PRJEB58186).

The presence of extrachromosomal circular copies of HOPPLA (RLC_BdisC024) was validated by an inverse PCR using 7 ng/μL total DNA. Input quantities of DNA were controlled using primers specific to the S-adenosylmethionine decarboxylase (SamDC) gene [START_REF] Hong | Exploring valid reference genes for gene expression studies in Brachypodium distachyon by real-time PCR[END_REF] . Sequences of primers are listed in S5 Table.

Analysis of Mobilome-seq

Reads were trimmed using the BBDuk tool of BBtools (BBMap (v 38.75,Bushnell B., sourceforge.net/projects/bbmap/) with the parameters qtrim = rl and trimq = 20. Reads originating from organelles were removed by aligning reads to the chloroplast genome (NC_011032.1) [START_REF] Bortiri | The complete chloroplast genome sequence of Brachypodium distachyon: Sequence comparison and phylogenetic analysis of eight grass plastomes[END_REF] and the mitochondrion genome (v 1.0.0) of B. distachyon using BWA-MEM (v 0.7.17-r1188) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] with the -M parameter set. Unmapped reads were isolated using samtools (v 1.13) [START_REF] Danecek | Twelve years of SAMtools and BCFtools[END_REF] view -b -f 4 and bedtools (v 2.30.0) [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] bamtofastq.

Organelle-filtered mobilome reads were assembled using the SPAdes genome assembler (v 3.13.0) [START_REF] Prjibelski | Using SPAdes De Novo Assembler[END_REF] . From each assembly, the top ten contigs were extracted and jointly aligned to the reference assembly of Bd21 (v 3.0) using BWA-MEM (v 0.7.17-r1188) with the -M parameter set. Bam files were converted into bed files using bedtools (v 2.30.0) bamtobed with the -split option set and overlapping contigs were merged using bedtools merge with the -o distinct, count, count_distinct and -c 4 parameters set. Assembled, circle-forming regions were annotated with bedtools intersect using the version 3.1 annotation of the reference assembly and the annotation of all full-length LTR-RTs [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] of the reference assembly. Annotated regions were extracted with bedtools getfasta and all sequences longer than 2 kb were isolated using SeqKit seq (v 0.11.0) [START_REF] Shen | SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation[END_REF] . Circle-forming regions that occurred in less than three samples were not included in the analysis.

To specifically detect mobilized LTR-RTs, we first extracted all annotated full-length LTR-RTs of the Bd21 reference assembly [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF]. Using a custom python script, we then merged the last 300 bp of the 3' to the first 300 bp of the 5' LTR to obtain a 'tail-to-head' library containing all annotated full-length LTR-RT copies annotated in Bd21. We then aligned organelle-filtered mobilome reads to the tail-to-head library of LTR-RTs and used bedtools (v 2.30.0) intersect to extract aligned reads that were spanning the 2-LTR junction and that aligned to at least 5 bp of both LTRs.

The coverage of the junction-spanning reads was calculated using deeptools (v 3.5.1) [START_REF] Ramírez | deepTools2: a next generation web server for deep-sequencing data analysis[END_REF] with the parameters -bs 1, --ignoreDuplicates --outRawCounts set. To account for differences in sequencing depth, the obtained coverage for 2-LTR-junction spanning reads was normalized with the total coverage obtained with bedtools (v 2.30.0) genomecov with the -d and -split parameters set, from the alignments of filtered reads to the reference assembly of Bd21 (v 3.0) obtained from Phytozome 12 (International Brachypodium Initiative, 2010) generated by BWA-MEM (v 0.7.17-r1188) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] . To plot the overall activity per family, normalized signals were summed up for every individual TE family.

mRNA-sequencing and small RNA northern blotting

Leaves of 4-week-old B. distachyon plants were ground in liquid nitrogen and 500 μL of this powder was subjected to TRIzol extraction following the supplier instructions (Invitrogen, CA, USA). 20 µg of total RNA was treated with DNase I for 30 min., then repurified via phenol-chloroform extraction and ethanol precipitation. DNase-treated total RNA samples were sent to Fasteris/Genesupport (Plan-les-Ouates, Switzerland), subjected to poly(A)-tail selection, and then aliquoted for library construction via the Illumina TruSeq Stranded mRNA Library Prep kit. Resulting stranded polyA+ RNA-seq (mRNA-seq) libraries were sequenced on an Illumina NovaSeq 6000. The raw paired-end read data were deposited at the NCBI Gene Expression Omnibus (GEO accession: GSE243693).

For the small RNA blot analysis, 200 µg of each total RNA were size-fractionated using the RNeasy Midi Kit (QIAGEN), as described previously [START_REF] Böhrer | Integrated Genome-Scale Analysis and Northern Blot Detection of Retrotransposon siRNAs Across Plant Species[END_REF] . Low molecular weight (LMW, <200 nt) RNAs are not bound by the silica membrane of the columns and were isolated from the collected flow-through and wash aliquots. LMW RNAs were precipitated overnight using isopropanol.

Following a centrifugation step (45 min. at 24000 x g, 4°C) and the removal of the supernatant, the pellet was washed with 75% ethanol, centrifuged (15 min. at 24000 x g, 4°C), dried at RT for 20 min. then at 65°C for 5 min., and resuspended in 41 µL of DEPC-treated MilliQ water. LMW RNAs were quantified using a Nanodrop device and 12.3 µg of LMW RNAs from each sample were loaded into the 16% polyacrylamide gel [START_REF] Böhrer | Integrated Genome-Scale Analysis and Northern Blot Detection of Retrotransposon siRNAs Across Plant Species[END_REF] . After running, transfer and UV crosslinking, membrane was prehybridized in PerfectHyb Plus buffer (Merck, Darmstadt, Germany) at 35°C and then and hybridized at 35°C with the Klenow internally-labeled probe (HOPPLA), or with the 5'-end labeled probe (miR160) [START_REF] Böhrer | Integrated Genome-Scale Analysis and Northern Blot Detection of Retrotransposon siRNAs Across Plant Species[END_REF] . After overnight hybridization, washing was performed at 37°C. Signal detection requires 5-7 days exposure for HOPPLA and 1-2 days for miR160.

Oligonucleotide sequences for the probes are listed in S5 Table.

LTR-RT expression analysis

RNA-seq raw reads of Bd nrpd1-2 (-/-) and Bd nrpd1-2 (+/+) were trimmed for adapters using fastp (v 0.23.2) [START_REF] Chen | Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp[END_REF] with the following options: --qualified_quality_phred 15 --unqualified_percent_limit 4 --n_base_limit 20 --low_complexity_filter --overrepresentation_analysis --correction --detect_adapter_for_pe. Cleaned reads were then analysed using SalmonTE (v 0.4) [START_REF] Jeong | An ultra-fast and scalable quantification pipeline for transposable elements from next generation sequencing data[END_REF] to measure global expression of LTR-RTs. LTR-RT consensus sequences of B. distachyon obtained from the TRansposable Elements Platform (TREP, https://trep-db.uzh.ch/) were used to generate the custom library for SalmonTE. Default options of SalmonTE quant and test function were used to quantify expression and to perform statistical analysis. Expression data were plotted using R (v 3.6.3) in RStudio (v 7d165dcf).

Motif analysis

The consensus sequence of HOPPLA was screened for known transcription factor binding sites obtained from the PlantTFDB [START_REF] Jin | PlantTFDB 4.0: Toward a central hub for transcription factors and regulatory interactions in plants[END_REF] using FIMO (v 5.1.1) [START_REF] Grant | FIMO: scanning for occurrences of a given motif[END_REF] . To functionally annotate transcription factors that could bind to HOPPLA, we used GO-terms of the Gramene (release 50) database [START_REF] Tello-Ruiz | Gramene 2021: Harnessing the power of comparative genomics and pathways for plant research[END_REF] downloaded from the platform agriGO (v 2.0) [START_REF] Tian | AgriGO v2.0: A GO analysis toolkit for the agricultural community, 2017 update[END_REF] . Generic, TF specific GO terms (GO:0003700, GO:0006355, GO:0005634, GO:0003677, GO:0043565, GO:0046983, GO:0003682 GO:0045893) such as 'positive regulation of DNA-templated transcription' were removed from the list of GO terms as they would interfere with the downstream analysis. The remaining GO terms of transcription factors potentially binding to HOPPLA were visualized with REVIGO [START_REF] Supek | REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms (C Gibas[END_REF] using the 'SimRel' semantic similarity measure, the option 'small' and the GO terms of the Oryza sativa Japonica Group. The total number of occurrences of individual GO terms was taken into account with the option 'higher value is better'.

GO terms occurring more than five times were labelled in the plots. As a proof of concept, we followed the exact same approach using the sequence of one of the most active ONSEN copies (AT1G11265) and the GO terms of Arabidopsis thaliana.

Detection of novel HOPPLA insertions in Bd nrpd1-2 (-/-)

DNA of adult plants was extracted using the DNeasy plant kit (Qiagen, Venlo, Netherlands) according to the protocol of the manufacturer subjected to whole genome sequencing. Reads were trimmed using the BBDuk tool of BBtools (BBMap (v 38.75,Bushnell B.,sourceforge.net/projects/bbmap/)

with the parameters qtrim = rl and trimq = 20. Trimmed reads were aligned to the reference assembly of Bd21 (v 3.0) using BWA-MEM (v 0.7.17-r1188) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] with the -M parameter set. Samtools (v 1.13) [START_REF] Danecek | Twelve years of SAMtools and BCFtools[END_REF] was used to obtain sorted and indexed bam files. TIPs were detected with detettore (v 2.0.3) (https://github.com/cstritt/detettore) with the options -require_split, -q 30 and using the consensus sequences of LTR-RTs of B. distachyon (TREPdatabase) and the annotation of all full-length LTR-RTs of Bd21 [START_REF] Stritt | Diversity, dynamics and effects of long terminal repeat retrotransposons in the model grass Brachypodium distachyon[END_REF] . Because both Bd nrpd1-2 lines were in the Bd21-3 background we were able to exclude all Bd21-3 specific TIPs by removing those insertions that were detected in multiple individuals with more than one genetic background. Remaining TIPs were manually curated using the genome browser IGV (v 2.15.4.12) [START_REF] Robinson | Integrative genomics viewer[END_REF] . HOPPLA TIPs were visualized with JBrowse 2 (v 2.6.1) [START_REF] Diesh | JBrowse 2: a modular genome browser with views of synteny and structural variation[END_REF] .

Raw genomic reads of the re-sequencing of Bd nrpd1-2 (-/-), Bd NRPD1 (+/+) and Bd21-3 have been uploaded to ENA (accession number PRJEB64053).

GWAS for pCNV

GEMMA 0.98.5 [START_REF] Zhou | Genome-wide efficient mixed-model analysis for association studies[END_REF] was used to test for associations between SNPs (Minadakis et al., 2023a) and the LTR-RT families pCVNs, while correcting for population structure [START_REF] Stritt | Migration without interbreeding: Evolutionary history of a highly selfing Mediterranean grass inferred from whole genomes[END_REF]Minadakis et al., 2023a) . Α centered relatedness matrix was first created with the option -gk 1 and association tests were performed using the option -maf 0.05 to exclude rare alleles, and the default SNP missingness threshold applied by GEMMA that excludes SNPs with missing data in more than 5% of the accessions. We selected 20 kb genomic regions with a 10 kb overlap that contained at least two SNPs above the False Discovery Rate of 0.05 or Bonferroni correction threshold as candidate region using the R package rehh (v 3.2.2) [START_REF] Gautier | rehh: an R package to detect footprints of selection in genome-wide SNP data from haplotype structure[END_REF] . Genes overlapping with candidate regions were selected with the BEDTOOLS (v 2.26.0) [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] intersect command using the version 3.1 of the B.

distachyon annotation file (https://phytozome-next.jgi.doe.gov) that are contained in the significant regions. The UpSetR [START_REF] Conway | UpSetR: an R package for the visualization of intersecting sets and their properties[END_REF] R package was used to visualize the intersections of significant genes between the variables. Protein constituents of the Pol IV and Pol V enzymes (see S4 Table ) were downloaded from the plant RNA polymerase database http://rna.polymerase.eu/.

(

  RLC) elements (Wilcoxon test, W = 6517781, p-value < 2.2e-16; Fig 1B, C). Furthermore, a PCA based on RLG pCNVs did not allow us to discriminate accessions based on their phylogenetic relationship, while a PCA performed with RLC elements separated the samples by genetic lineage (Fig 1D). The strongest result was found with a PCA performed with the five youngest and putatively most recently active LTR-RT families found in the pangenome of B. distachyon (the Angela families RLC_BdisC022, RLC_BdisC024, RLC_BdisC030, RLC_BdisC209 and the Tekay family RLG_BdisC004), with the first two axes together explaining more than 77.8% of the variance. Finally, with the exception of samples from the most recently diverged clades A_East and A_Italia (13 kya), they further allowed us to 148 discriminate samples based on the genetic clade of origin (Fig 1 B,D).

Fig 1 .

 1 Fig 1. Natural diversity of proximal copy number variation (pCNV) of LTR-RTs in B. distachyon. (A) Origin of the 320 natural accessions included in this study. Accessions that were used for the mobilome-seq are labelled in the map and a PCA with the bioclimatic variables of their place of origin is depicted. Colors of points correspond to the genetic clades whose estimated split is shown in the phylogenetic tree. Black points indicate that the accession cannot be clearly assigned to one genetic clade. (B) Heatmap with read counts of TE-derived sequences (proxy for the copy number variations, pCNVs) of all 40 annotated LTR-RTs in 320 natural accessions of B. distachyon. Accessions are sorted according to their phylogeny and names of recently active TE-families are highlighted in red. (C) Overall estimates of the copy numbers of Copia and Gypsy-type LTR-RTs in 320 natural accessions. (D) PCAs of pCNVs of all Gypsy, Copia and all recently active LTR-RT families, highlighted in (B) belonging to the Angela & Tekay families (RLG_Bdis004, RLC_BdisC030, RLC_BdisC209, RLC_BdisC024 and RLC_BdisC022). Colors of points indicate the genetic clade of accessions. (E) Output of the LMM analyses between pCNVs of LTR-RTs and bioclimatic variables at the accessions` origin. Bubbles indicate a significant association (P-value <0.05). Colors and sizes of bubbles show the part of the variance (marginal R 2 ) explained by the bioclimatic variables in %.

  and running classical linear model analyses (S1 Fig). With the exception of RLC_BdisC010 and RLG_BdisC265, for which more than 40% of the variance in pCNVs was explained by environmental factors, the LTR-RT families showed non-significant to mild associations with environmental variables (S1 Fig).

Fig 2

 2 Fig 2 Assessment of LTR-RT mobility in B. distachyon. Stress (A) and genotype (B) specificity of the formation of eccDNA as determined by the alignment of assembled mobilome-seq reads. The color represents the degree of specificity and numbers indicate the count of samples from which one of their ten longest contigs aligns to each of the circle-forming regions. Annotations of regions are indicated on the y-axis. Multiple annotations in the same circle-forming region were concatenated. The two pol IV mutants and the controls Bd21-3 and the outcrossed line BdNRPD1 (+/+) are summarized as RdDM and Bd21-3, respectively

Fig 3

 3 Fig 3 HOPPLA forms 2-LTR circles in the pol IV mutants. Rela:ve abundance of 2-LTR-junc:on spanning reads

Fig 3

 3 Fig 3 HOPPLA forms 2-LTR circles in the pol IV mutants. Relative abundance of 2-LTR-junction spanning reads depending on the stress (A) and the genotype (B) of individual mobilome-seq samples. (C) Inverse PCR using total DNA not subjected to a rolling circle amplification for the confirmation of an increased amount of extrachromosomal 2-LTR circles of HOPPLA in the Bd nrpd1-2 (-/-) mutant compared to the Bd NRPD1 (+/+) outcrossed line. Three biological replicates are shown. (D) GO enrichment analysis of transcription factors for which binding sites have been detected in the consensus sequence of HOPPLA. Colors indicate number of TF-binding sites found. GO-terms that occur at least six times are highlighted in the plot. All GO-terms and their number of occurrences are listed in S1 Table.

Fig 4

 4 Fig 4 Members of the HOPPLA family differ in activity. (A) Relative abundance of 2-LTR-junction spanning reads of individual full-length copies of the HOPPLA family. (B) Age (age_a), closest distance to gene (d_closest), GC content (GC), and methylation levels in CG, CHG and CHH contexts of all individual genomic full-length copies of HOPPLA in percent. The color indicates relative abundance of 2-LTR-junction spanning reads from the mobilome-seq in (A).

  two pol IV mutants and their respective wild-type controls. Using a hybridization probe specific to the HOPPLA LTRs, 24 nt siRNAs were detected in the control lines Bd21-3 and 262 Bd NRPD1 (+/+), but not in either of the pol IV mutant lines (Fig 5A). This finding strongly suggests that HOPPLA is under control of the Pol IV-RdDM pathway, and that the absence of 24 nt siRNAs 264 results in the upregulation and increased production of 2-LTR eccDNAs from HOPPLA. Furthermore, RNA-seq data from part of the same mutant panel shows that HOPPLA is the most upregulated LTR-266 RT family in the Bd nrpd1-2 (-/-) background compared to the Bd NRPD1 (+/+) control line, indicating that the reduction of 24 nt siRNAs is likely associated with an increased expression and subsequent 268 formation of HOPPLA eccDNAs in both pol IV mutants (Fig 5B).

Fig 5

 5 Fig 5 Loss of 24-nt siRNAs leads to an increased activity of HOPPLA in the Pol IV mutants (A) Northern plot for the detection of 24-nt siRNAs specific to the 3` LTR of HOPPLA in the pol IV mutants Bd nrpd1-1 and Bd nrpd1-2 (-/-) and their control lines Bd21-3 and the outcrossed line Bd NRPD1 (+/+). (B) SalmonTE analysis of the expression of LTR-RTs in Bd nrpd1-2 (-/-) relative to the outcrossed control line Bd NRPD1 (+/+). LTR-RTs with a log2 fold change of at least two are labeled, three biological replicates were analysed.

Fig 5

 5 Fig 5 Loss of 24-nt siRNAs leads to an increased acvity of HOPPLA in the Pol IV mutants (A) Northern plot for the

Fig 6

 6 Fig 6 Genetic regions associated with pCNVs of recently active LTR-RTs (A) Manhattan plot depicting the GWAS results of pCNV of HOPPLA in 320 accessions of B. distachyon. Colored points indicate SNPs linked to TIPs (blue), TAPs (red) or known components of the Pol IV and V holoenzymes (asterisks). Threshold of significance (false discovery rate adjusted p-value <0.05) is marked with dashed lines. A significant region containing the candidate gene Bradi5g05225 (window size 50 kb) is highlighted. (B) UpSet plot of genes in 20 kb windows surrounding significant regions with at least two SNPs above the threshold of significance (FDR-adjusted p-value <0.05 for HOPPLA, RLC_BdisC209 and RLC_BdisC022 and Bonferroni correction for RLC_BdisC030 and RLG_BdisC004) of the five most recently active LTR-RT families in B. distachyon. To visualize potential overlaps, a list of the componenlangts of the Pol IV and Pol V holoenzymes is included in the UpSet plot.

AFig 6

 6 Fig 6 Genec regions associated with pCNVs of recently acve LTR-RTs (A) Manha&an plot depic:ng the GWAS results of pCNV of HOPPLA in 320 accessions of B. distachyon. Colored points indicate SNPs linked to TIPs (blue), TAPs (red) or known components of the Pol IV and V holoenzymes (asterisks). Threshold of signiHcance (false discovery rate

  five families (see S3 Table

S2Fig.

  TFs binding to ONSEN are heat-stress inducible. GO-enrichment analysis of transcription factors for which binding sites have been detected in AT1G11265, a member of the heat-responsive ONSEN (ATCOPIA78) LTR-RT family in A. thaliana. Colors indicate number of TF-binding sites found. GO terms that occur at least six times are highlighted in the plot. All GO-terms and their number of occurrences is listed in S1 Table.

  Figs. TIPs of HOPPLA in one of the resequenced Bd nrpd1-2 (-/-) plants. JBrowse screenshot of the insertion site in Bd nrpd1-2 (-/-) (bottom) compared to the Bd21-3 wt (top). The target side duplication (TSD) is annotated and soft clipped parts of reads are coloured.

Acknowledgements and author contributions

We thank the Genetic Diversity Center Zürich for providing the infrastructure for the mobilome-seq. M.T. and A.C.R. conceived the study and planned experiments. M.T., B.K., K.R., C.H., B.R. and M.B. generated data. M.T., A.C.R., N.M. and W.X. analysed data. D.L., J.V. and R.S. provided mutagenized germplasm and online TILLING/T-DNA browsing tools. C.H., B.R. and T.B. identified and screened for the pol IV null mutants. M.T and A.C.R wrote the paper with contributions from T.B., C.S. All authors read the manuscript. The authors declare that they have no conflict of interest. Professional editing performed by Emanuelle Botté.

Funding

University of Zurich Research Priority Programs (URPP) Evolution in Action to M.T. and A.C.R.; Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung. Grant Number: 31003A_182785 to A.C.R., W.X., N.M. and B.K.; Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung. Grant Number: PZ00P3_154724 to C.S.; T.B. was supported by the Interdisciplinary Thematic Institute IMCBio (ITI 2021-2028 program), including funds from IdEx Unistra (ANR-10-IDEX-0002), SFRI-STRAT'US (ANR 20-SFRI-0012) and EUR IMCBio (ANR-17-EURE-0023) in the framework of the French Investments for the Future Program.