Plant materials and growth conditions
The A. thaliana Columbia-0 (Col-0) accession of wild type and the ddm1-1 mutant30 were used throughout. The seeds were germinated on a plate of 0.5× or 1× Murashige and Skoog medium. After keeping in dark at 4 °C for 2–3 days, the plants were grown at 22 °C under long-day conditions of 16 h of light and 8 h of dark. Plant samples were harvested directly from the plate or transferred to soil for genetics.
TEd-seq library construction
A detailed protocol for TEd-seq50 can be found at https://www.protocols.io/view/ted-seq-c7seznbe. Before preparation of the TEd-seq library, 30 µM of custom oligonucleotide adaptor was prepared as follows: 90 µl of 100 µM P7_adapter_up (5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATC*T-3′, *phosphorothioate bond), 90 µl of 100 µM P7_adapter_bottom (5′-pGATCGGAAGAGCATC**−3′, p denotes phosphorylation, **dideoxy-C), 30 µl of T4 ligase buffer and 90 µl of H2O were mixed. The mixture was heated at 95 °C for 2 min and cooled to 25 °C over 45 min. The TEd-seq library was prepared as follows. In most samples (see Supplementary Table 3 for details), genomic DNA was extracted from roughly 30 10-day-old seedlings using Nucleon Phytopure DNA extraction kit (GE Healthcare). Next. 500 ng of genomic DNA was suspended in 100 µl of TE buffer and sheared into the length of 250–500 bp using Bioruptor Pico (Diagenode) for three cycles of 30 s on/30 s off. The sheared DNA was cleaned up using 0.9× SPRIselect (Beckman Coulter), resuspended in 25 µl of 0.1× TE buffer and subsequently END-repaired, A-tailed and ligated to custom oligonucleotide adaptors using NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB) with the reaction scale by half. Adaptor-ligated DNA was size selected using 0.13× SPRIselect for first and second bead addition, and eluted in 15 µl of 10 mM Tris-HCl pH 8.0. Then, nested PCR was performed to enrich fragments containing the terminal sequence of specific TE. For the first PCR, 7.5 µl of purified adaptor ligated DNA, 2.5 µl of adaptor specific primer, 2.5 µl of TE specific primer, 12.5 µl of NEBNext Ultra II Q5 Master Mix were mixed and PCR amplified with the following condition: 98 °C for 30 s, 20 cycles of (98 °C for 10 s, 61 °C for 75 s), 61 °C for 5 min. The primer sequences used in first PCR are shown in Supplementary Table 4. The first PCR product was cleaned up using 0.9× SPRIselect and eluted in 15 µl of 0.1× TE. The purified first PCR product was diluted ten times with 0.1× TE and used as a template in second PCR. In the second PCR, P7_primers with index and P5_TE_primers with index shown in Supplementary Table 4 were used for dual indexing. Then, 2.5 µl of first PCR product, 6.25 µl of NEBNext Ultra II Q5 Master Mix, 1.25 µl of P7_primer with index, 1.25 µl of P5_TE_primer with index, 1.25 µl of H2O were mixed and amplified under the following conditions: 98 °C for 30 s, two cycles of 98 °C for 10 s, 61 °C for 75 s, eight cycles of 98 °C 10 s, 72 °C for 75 s, 72 °C for 5 min. The second PCR product was cleaned up with 0.9× SPRIselect and diluted in 15 µl of 0.1× TE. The final library product should contain 41 and 61 bp of the 5′ LTR terminal sequences of Tal1 and Evade, respectively. Sequencing was performed at Macrogen Japan Corp., using the Illumina Hiseq X platform or Novaseq X with a read length of 151 bp. The detailed experimental conditions for library preparation are shown in Supplementary Table 3.
TEd-seq analysis
Source code for the TEd-seq analysis can be accessed at (https://github.com/LeanQ/TED_seq_Tsukahara_2024). Briefly, TEd-seq fastq files were mapped to the 5′ terminal sequence of TE (1–144 bp of Evade or 1–105 bp of Tal1 shown in the file of ‘target_TE_sequence_extremity.fa’) using Bowtie2 (v.2.5.3)51 with the parameter ‘–local –very-sensitive’. Pair-end reads in which only one mate mapped over 5′ terminal sequence of TE were extracted using Picard tools (v.2.27.5) (https://broadinstitute.github.io/picard/) with the function of ‘FilterSamReads’ and then extracted discordantly mapped reads. Those reads were then mapped to Col-CEN_v1.2 reference genome using Bowtie2 with the parameter ‘–local –very-sensitive’, and the mapped data was converted to BAM files using SAMtools (v.1.9) to generate ‘clip_disc-local.sorted.bam’ files. To detect de novo somatic insertions of TEs, clip_disc-local.sorted.bam files were converted to bedfiles using BEDTools (v.2.31.1)52 with the ‘bamtobed’ function. The sequence reads that contain de novo insertion site should be soft-clipped reads, in which the terminal sequence of TE is soft-clipped. Within the soft-clipped reads, the start position of alignment to the reference, which is the flanking nucleotide of soft-clipped region was regarded as an insertion site. The reads mapped within the upstream 1 kilobase (kb) and downstream 1 kb of the regions that are annotated as Evade (AT5TE20395) and its related copy (AT1TE41580) (shown in the file of ‘targeted_TE_sequences.bed’) were removed using the ‘intersect’ function of BEDTools, as they are not de novo insertions. The reads mapped to 1–30 kb of chromosome 2 were also removed as they were systematically found in all samples, including wild-type controls, and therefore do not correspond to bona fide de novo insertions. To obtain genome-wide integration landscapes, the number of insertion sites without duplicates were counted in non-overlapping 10 kilobase pair (kbp) windows genome wide using the ‘coverage’ function of BEDTools. The number of insertion sites were plotted using the packages of ggplot2 (v.3.4.4)53, readr (v.2.1.5)54 and dplyr (v.1.1.4)55 in R software (v.4.3.2)56.
DNA extraction for PacBio sequencing
Genomic DNA was extracted from 2 g of aerial parts of 25-day-old Tal1 transgenic plants (ddm1 mutant background), which were put in dark place for 2 days before sampling. Genomic DNA was extracted by the method of Carlson lysis buffer containing cetyl trimethyl ammonium bromide (CTAB)57 using Genomic-tip (Qiagen) in the following conditions. The frozen plant tissue was ground with a mortar and a pestle with liquid nitrogen into fine powder, and put it into Carlson lysis buffer (100 mM Tris-HCl pH 9.5, 20 mM EDTA, 1.4 M NaCl, 1% PEG 6000, 2% CTAB, 0.1% b-mercaptoethanol) heated at 74 °C in advance and incubated at 74 °C for 20 min inverting every 5 min. After incubation, when the sample was cooled down to room temperature, 20 ml of chloroform:isoamyalcohol (24:1) was added and mixed until homogenized and centrifuged at 3,000 rpm for 10 min. Then 20 ml of chloroform:isoamyalcohol (24:1) was added to the supernatant and mixed until homogenized. The sample was centrifuged at 3,000 rpm for 1 min and 20 ml of 2-propanol was added, mixed and incubated at 4 °C overnight. The sample was centrifuged at 3,500 rpm for 30 min and the supernatant was discarded. Next, 70% ethanol was added and centrifuged at 3,500 rpm for 10 min and the supernatant was discarded completely. The pellet was suspended with 1 ml of TE buffer on ice and 9 ml of Buffer G2 was added and mixed. Then, 18 µl of 100 mg ml−1 RNase A was added to the tube and incubated at 37 °C for 30 min. Next, 90 µl of proteinase K (Qiagen) was added to the sample and incubated for 50 °C for 1 h. After centrifugation at 10,000g, the supernatant was purified with Genomic-tip_100/G (Qiagen) following the manufacturer’s protocol. Then, 0.7× volume of 2-propanol was added to the eluted DNA, inverted several times and centrifuged at 10,000g for 20 min. Next, 70% ethanol was added to the pellet and centrifuged at 10,000g for 20 min. The pellet was dissolved in 100 µl of TE buffer. Extracted DNA was quantified with Qubit double-stranded DNA High Sensitivity Assay kit (Thermo Fisher Scientific) and Nanodrop 2000 (Thermo Fisher Scientific).
Library preparation and analysis of PacBio sequencing
Extracted DNA was sheared with g-tube (Covaris). The SMRT library was prepared by SMRTbell express template prep kit (Pacific Bioscience) and the libraries were size selected (greater than 30 kb) using BluePippin system (Saga Science). The libraries were sequenced by PacBio Sequel (Pacific Bioscience) using Binding Kit v.3.0 and Sequencing Kit v.3.0 with the continuous long-read sequencing mode. One SMRT cell (Sequel SMRT Cell 1 M v.3 LR) was used for each library and ran for 20 h. The SMRT sequencing data was generated at a coverage of 101×, N50 of 30,210 bp and a mean read length of 16,298 bp. The longest subreads were extracted from the sequencing data (682,486 subreads). The extracted subreads were mapped to 484 bp of the LTR sequence of Tal1 (1,024 reads). The reads with less than 200 nucleotides (nt) of mapped region were removed from the analysis. The rest of the reads were mapped to Col-CEN reference genome3 (ASM2311539v1 (GCA_023115395.1)) and 917 reads were mapped. Minimap2 (v.2.15-r905) was used for the mapping with the parameter ‘-c -x map-pb’. Among the 917 of mapped reads of Tal1(ddm1), the reads that had more than 3,560 bp (the length of 20 copies of CEN178) of the genomic sequence flanking Tal1 were extracted. Many Tal1 insertions are often detected in each PacBio read, but for the analyses here, we used only the best Tal1 insertion match from each read. The number of Tal1-inserted site per 100 kbp of genomic sequence were counted.
Construction of CENH3 overexpression line
CENH3 overexpression line was generated as follows. For CENH3-OX-1, 2 and 4, CENH3 (At1g01370) coding region including intron was amplified from genomic DNA of wild-type Col-0. The amplified CENH3 fragment was introduced into XmaI/BamHI-digested pPLV01-pRPS5a vector that was generated by inserting pRPS5a promoter sequence in the HpaI site of pPLV01 vector58. For CENH3-OX-3, upstream region and coding region of CENH3 was amplified from genomic DNA of wild-type Col-0 and was introduced into XhoI/BamHI-digested pPLV01 vector. The primers used for the constructions are shown in Supplementary Table 4. The amplified CENH3 fragment and the digested vector were assembled using NEBuilder HiFi assembly Master Mix (NEB). The assembled product was introduced into the Escherichia coli DH5a strain by the heat shock method. The extracted plasmid was introduced into Agrobacterium tumefaciens GV3101::pMP90 by electroporation. The agrobacterium with transgene was introduced into Arabidopsis wild-type Col-0 by the floral dip method59. Transgenic T1 plants were selected in the Murashige and Skoog medium with 50 µg ml−1 of Basta. Plants with a homozygous transgene were selected in T2 plants. T2 or T3 plants were used for chromatin immunoprecipitation with sequencing (ChIP–seq).
Western blotting
Here, 0.5 g of 2-week-old seedlings from non-transgenic line (wild-type Col-0) or CENH3 overexpressing lines were ground into fine powder with liquid nitrogen, and nuclei were isolated by the method described previously60. Proteins were separated by 15% SDS–PAGE and then transferred to a polyvinyl difluoride (Cytiva) membrane using Trans-Blot SD Semi-Dry Cell (Bio-Rad). Primary antibodies against H3 (0.2 µg ml−1; Abcam, ab1791), CENH3 C-terminal (0.2 µg ml−1; affinity-purified rabbit polyclonal antibody against the peptide CRKDFELARRLGGKGRPW), HTR12 (CENH3 N-terminal) (0.25 µg ml−1; affinity-purified rabbit polyclonal antibody against the peptide RTKHRVTRSQPRNQTDAC) and H4 (0.24 µg ml−1; affinity-purified rabbit polyclonal antibody against the peptide CKRQGRTLYGFGG), and peroxidase-linked secondary antibody against rabbit IgG (1:10,000 dilution; Cytiva, NA934) were used for western blotting. Signals were developed using Western BLoT Quant Horse Radish Peroxidase Substrate (Takara) and detected using iBright Imaging System (Thermo Fisher Scientific).
Expression and purification of recombinant Arabidopsis H3 and CENH3
The DNA fragment encoding AtCENH3 was inserted into the pET-15b vector (Novagen), in which the tobacco etch virus protease recognition site was introduced instead of removal of the thrombin recognition site. The expression and purification of A. thaliana H3.1 and CENH3 proteins were performed as described previously61. Purified recombinant AtH3 and AtCENH3 proteins (50 ng each) were used for the validation of antibodies used in this study.
CENH3 ChIP–seq
Here, 0.4–0.5 g of 2-week-old whole seedlings were frozen with liquid nitrogen, ground into fine powder with motor and pestle and lysed with fixing buffer (PBS with 1% formaldehyde, 0.3% Triton X-100, 1 mM Pefabloc SC (Roche) and cOmplete EDTA-free Protease Inhibitor Cocktail (Roche)). The suspension was rotated at room temperature for 10 min for crosslinking. Then 0.2 M of glycine was added for quenching and rotated for 5 min at room temperature. The sample was centrifuged at 5,000g, 4 °C for 5 min, the pellet was washed once with PBS and it was resuspended with low-salt ChIP buffer without Triton X-100 (50 mM HEPES-KOH, 150 mM NaCl, 1 mM EDTA, 0.1% sodium deoxycholate, 0.1% SDS) to make the volume 900 µl. The sample was divided equally into three tubes, and sheared into 200–600 bp sizes using the Picoruptor sonication device (Diagenode) with 12 cycles of 30 s on/30 s off. The sonicated samples were centrifuged at 20,000g at 4 °C for 10 min, and Triton X-100 was added to the supernatant (final concentration: 1%). The sonicated chromatin was incubated with 0.5 µg of anti-HTR12 (CENH3 N-terminal) antibody62 overnight at 4 °C, and then incubated with Dynabeads Protein G (Veritas) at 4 °C for 2 h. The incubated beads were washed once with 1 ml of low-salt ChIP buffer (50 mM HEPES-KOH, 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 0.1% SDS), twice with 1 ml of high-salt ChIP buffer (50 mM HEPES-KOH, 500 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 0.1% SDS), once with 1 ml of LiCl buffer (10 mM Tris-HCl, pH 7.8, 1 mM EDTA, 0.25 M LiCl, 1% IGEPAL CA-630, 1% sodium deoxycholate) and once with 1 ml of TE buffer (10 mM Tris-HCl, pH 7.8, 1 mM EDTA), while rotating for 10 min at 4 °C each time. DNA was eluted in 100 µl of elution buffer (10 mM Tris-HCl, pH 7.8, 0.3 M NaCl, 5 mM EDTA, 0.5% SDS), incubated overnight at 65 °C and then purified with the Monarch PCR & DNA Cleanup kit (New England Biolabs). The DNA was quantified with the Qubit dsDNA High Sensitivity Assay kit (Thermo Fisher Scientific), and 1–2 ng of DNA was used for preparing library with ThruPLEX DNA-seq kit (Clontech). The libraries were sequenced using HiSeq X Ten sequencer (Illumina) or Novaseq X plus sequencer at Macrogen Japan Corp.
ChIP–seq analysis
ChIP–seq data was processed as previously described63 with some modifications. The fastq file of ChIP–seq data was quality filtered using Trimmomatic (v.0.39)64 and mapped to Col-CEN_v1.2 reference genome. Mapping was conducted with Bowtie2 (v.2.4.4)51 with no option. The mapped data was converted to BAM files using SAMtools (v.1.6)65 and converted to BED files using BEDTools (v.2.26.0)52. The number of reads overlapped with 10 kb bin of chromosomes were counted using coverage function of BEDTools to make bedgraph files. Signals were plotted in units of reads per million mapped reads using ggplot2 (v.3.4.4)53 package in R software (v.4.3.2)56. Reads that aligned to the chloroplast or the mitochondrial genome were not included in normalization. In Figs. 2c–f and 3d, a few (less than seven in each panel) dots are out of the ranges shown. The range of TR region around each centromere follow previous analyses3. The range of pericentromeric (PC) region surrounding each centromeric TR is defined as continuous 10 kb block with average mCHG > 0.25.
DNA methylation profiling by ONT
DNA methylation profiling of genomes including the centromeric TR regions has been performed as described previously3. Three weeks old seeding were ground in liquid nitrogen and the DNA was extracted using NucleoBond HMW DNA kit (MACHEREY-NAGEL) as per the manufacturer’s instructions. A sequencing library was prepared using LSK-110 ligation kit, and sequenced using two MinION R9 flowcells. Methylation calling was performed as previously reported66 with the modification that R9 reads were filtered for length and accuracy using Filtlong (v.0.2.0) (–min_mean_q 90, –min_length 5,000) and methylation prediction for the CG, CHG and CHH contexts were called using DeepSignal-plant (v.1.6.1) using the model: model.dp2.CNN.arabnrice2-1_120m_R9.4plus_tem.bn13_sn16.both_bilstm.epoch6.ckpt.
Constructions of chimeric retrotransposons
To generate each construct of chimeric retrotransposons, DNA fragments were amplified by PCR using the plasmid Tal1_pRI909 (ref. 27) or EVADE_pRI909 as templates. To obtain Evade_pRI909 that contains a full length of EVADE (AT5TE20395, 5,329 bp), the 5′ half and 3′ half of AT5TE20395 were amplified by nested PCR from A. thaliana Col-0 genome using primers with the sites of restriction enzyme as shown in Supplementary Table 4. Second PCR products were digested with restriction enzymes: PstI and EcoRI for the 5′ half fragment, and EcoRI and BamHI for the 3′ half fragment. Both fragments and pRI909 binary vector (Takara) digested by PstI and BamHI were ligated with Mighty Mix (Takara). The plasmids and primers used in the construction of the chimeric retrotransposons are shown in Supplementary Table 4. PCR-amplified fragments and HpaI-linearized pPLV02 vector or pPLV03 (ref. 58) vector were separated by electrophoresis. The gel bands were extracted and purified by QIAquick Gel Extraction Kit (Qiagen) or Fastgene Gel/PCR Extraction Kit (Fastgene). Purified DNA fragments were assembled into a pPLV02 or pPLV03 vector using NEBuilder HiFi DNA assembly (New England Biolabs) and cloned in E. coli. Transgenic lines were generated by Agrobacterium-mediated gene transfer as described previously27. To confirm the activation of the introduced TE, the presence of extrachromosomal circular DNA of the TE in the transgenic plants was examined by PCR67. The primers and conditions used for the PCR are shown in Supplementary Table 4.
Analysis of nucleotide sequence bias around integration sites of EVD and Tal1
To examine local integration bias of EVD and Tal1, the region around integration site of each TE was extracted as follows. TE sequence in Read1 of TEd-seq reads (61 bp of 5′ terminals sequence for EVD or 41 bp of 5′ terminal sequence for Tal1) were trimmed with Cutadapt 4.4 (ref. 68) with the parameter ‘-gGCCCACTCTCTTGTAGTACATATCCAATACTAGGCCTTTCTTATTTGAGTCTTGATCAATA-m 30’ for Evade and ‘-gATGTACGGATGGGTGCTTCACTCTTCGTTTCTTGATCAATA-m 30’ for Tal1. The trimmed reads were mapped to the Col-CEN_v1.2 reference genome using Bowtie2 (v.2.5.1) with the parameter ‘–local –very-sensitive’, and the mapped data was converted to BAM files using SAMtools (v.1.17). The central position of 5 bp target site duplication shown in Extended Data Fig. 9 was regarded as position zero. Next, 11 bp of integration site (−5 to +5 nucleotides from the centre of target site duplication) was extracted from Col-CEN_v1.2 using the ‘getfasta’ function of BEDTools (v.2.31.0). The number and the ratio of each nucleotide at each position in the 11 bp of integrated region were calculated.
Annotation of Arabidopsis LTR elements
TEs were identified for the NT1 and MN47 natural accessions of A. lyrata and for the Col-CEN reference genome of A. thaliana using the the Extensive De novo TE Annotator pipeline (v.2.0.1)69 with parameters –anno 1 and –sensitive 1, and the Repbase70 Arabidopsis-specific TE library (athrep.ref) as –curatedlib. Intact Ty1/Copia (110 in Col-CEN, 699 in NT1 and 800 in MN47) and Ty3 (96, 830 and 1,029) LTR retrotransposons were further classified into lineages using TEsorter 1.3 (-db rexdb-plant -nolib)71. For all analysis, we removed elements that did not have consistent superfamily classification between Extensive De novo TE Annotator and TEsorter (for example, an element needed to be classified as Ty3 by both pipelines). Following this approach, we retrieved 55, 393 and 295 ALE in Col-CEN, NT1 and MN47, respectively, and 49, 340 and 287 ATHILA.
Phylogenetic analysis and characterization of ALE and ATHILA
For the phylogenetic analysis, we further kept the subset of ALE (47, 368, 276) and ATHILA (11, 112, 71) that contained in correct order the hidden Markov models of all five genes (gag, protease, integrase, reverse transcriptase and RNaseH) according to TEsorter. For ATHILA, we also kept elements that contained only the gag and protease genes (35, 207, 205), because they represent a non-autonomous deletion derivative that is found in high numbers in A. thaliana and A. lyrata5. For ALE, we aligned the concatenated integrase and reverse transcriptase hidden Markov model domains retrieved from TEsorter for ALE with MAFFT (v.7.453, –globalpair –maxiterate 1,000)72. For ATHILA, we aligned the full-length DNA sequences with MAFFT (v.7.453, –retree 2 –maxiterate 50), so that elements that lacked the integrase, reverse transcriptase and RNaseH could be included in the tree. We used FastTree (v.2.1.11)73 with default parameters to generate maximum-likelihood trees.
The age of LTR retrotransposons was estimated using the sequence divergence between the two LTRs of an intact element. A pairwise alignment was produced for each pair of LTR sequences using MUSCLE (v.3.8.1551)74. We used the recent dataset of ATHILA based on the annotation of 66 A. thaliana accessions5. To identify the R/K polymorphism in the ALE4 branch, we looked for open reading frames with getorf from the EMBOSS:6.6.0.0 (ref. 75) using the internal region of every element. As in the phylogenetic analysis, we only surveyed elements that included all five genes. We focused on the longest open reading frame to examine and catalogue the R or K presence at the C-terminal end of the integrase gene with a custom R script. Downstream analyses were performed and plotted in R language and environment56, using packages included in Tidyverse collection (v.2.0.0)76, circlize (v.0.4.15)77 and ggtree (v.3.10.0)78.
Characterization of ALE4-like elements in other species
Sequences related to EVD and Tal1 were retrieved from Barbarea vulgaris (GCA_963667165.1), Eutrema japonicum (GCA_030161315.1 and GCA_030161335.1)79 and Raphanus sativus (GCF_000801105.2)80. The reverse transcriptase core domain region from A. thaliana COPIA20 was used to query the homology search by tBLASTn. Sequences with values less than 1 × 10−50 (for B. vulgaris and R. sativus) or 1 × 10−80 (for E. japonicum) were retrieved from genomic sequence assembly for analyses. Sequences including both 10 kb regions from the BLAST hit regions were obtained from genome assemblies to check presence of satellite sequences. The reverse transcription core domain regions were aligned for constructing phylogenetic trees. Aligned sequences were checked manually to delete sequences with more than 100 bp ambiguous or missing sites. Phylogenetic trees for reverse transcription regions were constructed by the neighbour-joining method with p-distances. EVD, AT1TE54585 (COPIA93) and AT2TE13385 (COPIA20) from A. thaliana and consensus sequences from clusters 1 to 4 of copia93/20 sequences from A. lyrata with other A. thaliana COPIA families (COPIA92, COPIA67, COPIA5, COPIA49, COPIA33, COPIA69, COPIA63, COPIA45 and COPIA35) were also included. All phylogenetic analyses were done by MEGA v.7.0 (ref. 81).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
