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

Abstract

Long terminal repeat retrotransposons (LTR-RTs) are powerful mutagens regarded as a major source of genetic novelty and important drivers of evolution. Yet, the uncontrolled and potentially selfish proliferation of LTR-RTs can lead to deleterious mutations and genome instability, with large fitness costs for their host. While population genomics data suggest that an ongoing LTR-RT mobility is common in many species, the understanding of their dual role in evolution is limited. Here, we harness the genetic diversity of 320 sequenced natural accessions of the Mediterranean grass Brachypodium distachyon to characterize how genetic and environmental factors influence plant LTR-RT dynamics in the wild. When combining a coverage-based approach to estimate global LTR-RT copy number variations with mobilome-sequencing of nine accessions exposed to eight different stresses, we find little evidence for a major role of environmental factors in LTR-RT accumulations in B. distachyon natural accessions. Instead, we show that loss of RNA polymerase IV (Pol IV), which mediates RNA-directed DNA methylation in plants, results in high transcriptional and transpositional activities of RLC_BdisC024 (HOPPLA) LTR-RT family elements, and that these effects are not stress-specific. This work supports findings indicating an ongoing mobility in B. distachyon and reveals that host RNA-directed DNA methylation rather than environmental factors controls their mobility in this wild grass model.

Fig 1. Natural diversity of proximal copy number (pCN) variation 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. 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. Numbers in dots indicate how many sequenced natural accessions were sampled in the marked area. The map has been obtained from (https://www.naturalearthdata.com/downloads/50m-cross-blend-hypso/50m-cross-blended-hypso-with-shaded-relief-and-water/). (B) Heatmap with read counts of TE-derived sequences (proxy for the copy number (pCN) variation) of all 40 annotated LTR-RTs in 320 natural accessions of B. distachyon. TE-families are clustered by their pCNs (dendrogram above the heatmap) and accessions are sorted according to their phylogeny. Names of recently active TE-families are highlighted in red. (C) Overall estimates of the copy numbers of Ty1/Copia and Ty3-type LTR-RTs in 320 natural accessions. (D) PCAs of pCNs of all Ty3, Ty1/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 pCNs 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²) explained by the bioclimatic variables in %.

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. Loci containing HOPPLA (RLC_BdisC024) are highlighted in red. Multiple annotations in the same circle-forming region were concatenated. The two pol IV mutants and the controls Bd21-3 and the outcrossed line Bd NRPD1 (+/+) are summarized as RdDM and Bd21-3, respectively. The following stresses were applied: c (control conditions), cold (2°C on ice, 24h), drought (uprooting, 2:15 h), drug (chemical de-methylation with Zebularine (20 uM) and alpha-amanitin (2.5 mg/ml), 28 days), glyphosate (20 mM, four days), heat (42°C, 8h), rice blast (Magnaporthe oryzae infection, four days), salt (300 mM NaCl, five days) and submergence (48 h). See materials and methods for details.

Fig 3. HOPPLA forms 2-LTR circles in the pol IV mutants.

Normalized abundance of 2-LTR-junction spanning reads depending on the stress (A) and the genotype (B) of individual mobilome-seq samples. HOPPLA (RLC_BdisC024) is highlighted in red. The alignment of publicly available genomic reads (genomic) served as a control for the presence of genomic reads aligning to 2-LTR-junctions. See Fig 2 and methods section for details about the applied c (control) conditions and stresses. (C) Coverage of mobilome-seq reads of the HOPPLA consensus sequence of sample Bd nrpd1-2 (-/-) submergence stress. The structure of the HOPPLA element is depicted. (D) 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. Loaded PCR reactions with primers specific to the HOPPLA LTRs (top gel) and to the genomic control gene SamDC (Bradi5g14640) (bottom gel) are shown. Three biological replicates are depicted. Schematic representation of primer design for the inverse PCR is shown.

Fig 4. Members of the HOPPLA family differ in activity.

(A) Normalized abundance of 2-LTR-junction spanning reads of individual full-length copies of the HOPPLA family. Samples with a high signal are labelled with c (control conditions) or the respective stress applied. See Fig 2 and materials and methods for details. (B) Age, closest distance to gene, GC content, 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). (C) 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 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, HOPPLA is highlighted in red, three biological replicates were analysed.

Fig 6. Genetic regions associated with pCN variations of recently active LTR-RTs.

(A) Manhattan plot depicting the GWAS results of pCN variation of HOPPLA in 320 accessions of B. distachyon. Colored points indicate SNPs linked (+/- 10 kb) to known Pol IV and V subunits (dark grey), HOPPLA reference insertions and TIPs (orange) and the genomic region containing the candidate gene Bradi5g05225 (red). 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 components of the Pol IV and Pol V holoenzymes is included in the UpSet plot.