Recent and dynamic transposable elements contribute to genomic divergence under asexuality

Plant material and DNA extractions

We selected five natural accessions of the apomictic T. officinale lineage Macranthoides to document and compare their repetitive genomic fractions. Macranthoides is a geographically widespread apomictic dandelion lineage, also referred to as a microspecies, that can be identified by expert taxonomist and whose unique clonal nature was confirmed by identical multi-locus microsatellite genotypes [25]. The same five accessions were previously screened for heritable divergence in gene expression using RNA-Seq [26]. The accessions were identified by Taraxacum taxonomists and sampled in Germany and Czech Republic from 3 different field locations (maximum distance 300 km) during spring 2012 (Additional file 1). Seeds were propagated for two generations in a common greenhouse environment. For each generation, seeds were surface-sterilized with a 0.5 % Sodium Hypochlorite solution and were germinated on agar (0.8 %) in petri dishes for 10 days in a climate chamber (10 h dark/14 h light, 15/20 °C). Then, seedlings were transplanted to individual pots (80 % potting soil with 20 % pumice) and grown in the greenhouse until leaf collection. For subsequent DNA extractions, 10 discs from two young leaves of the same individual were punched, pooled in one tube and directly frozen in liquid nitrogen. Samples were stored at −80 °C until DNA extractions.

Genomic DNA was isolated using NucleoSpin Plant II (Macherey-Nagel) columns following the manufacturer’s protocol with modifications as follows: leaf tissue was finely ground using liquid nitrogen before adding PL1 cell lysis buffer and RNase A solution. Lysate was incubated at 65 °C for 30 min and centrifuged for 2 min at full speed. Only supernatant was loaded onto the column. Quality and concentration were checked on a NanoDrop 1000 spectrophotometer.

Illumina sequencing

Sequencing of the five accession total genomic DNA was conducted by the CAT-AgroFood sequencing facility (Wageningen University, NL). Libraries with average insert size of 500 bp were prepared and pooled for subsequent 125 bp paired-end sequencing using half a lane of an Illumina HiSeq 2500 system. Genomic reads were analysed and trimmed for adapters and low-quality sequences using FastQC quality control. Cleaned sequencing output for each accession led to 41.4, 32.8, 55.0, 59.4 and 51.8 million reads for accessions 8, 11, 3, 12 and 13 respectively. From each accession a random sample of 5 million reads (125 bp length) was drawn from the quality-filtered sequencing output for subsequent downstream TE analysis that is comparable between the five accessions, representing 24 % of T. officinale triploid genome.

De novo identification of repetitive sequence families in Taraxacum

The RepeatExplorer pipeline via a Galaxy interface was applied to our datasets in order to cluster and assemble consensus sequences (contigs) of highly repeated fractions of the genome [27]. The five accessions were analysed separately. Read clustering criteria included 55 % minimum overlap and 90 % minimum similarity in the overlapping region. All clusters containing at least 0.2 % of the input reads were examined manually to identify clusters that required merging following recommendations and parameter settings from Kelly et al. [28]. Clusters were annotated according to similarity searches against the RepeatMasker Viridiplantae Database [29], the Conserved Domain Database [30] and GenBank non-redundant genomic database [31]. Chloroplast sequences were also clustered and excluded from subsequent analyses. After merging same-annotation clusters, all clusters representing >0.2 % of the input reads were explored manually to check for consistency in repeat annotation. This subset of high-abundant repeats was used for all subsequent analyses. To calculate the different repeat fraction sizes, total read lengths of same-annotation clusters representing similar TE clades, sub-clades, superfamilies and when possible families, were summed up [1]. Genome proportions were calculated based on the average monoploid genome size i.e. 865 Mb as calculated by Záveský et al. [32]. Differences in relative abundance of TE classes between accessions, as estimated by read counts, were tested using the likelihood ratio Chi-square (SAS 9.2, The SAS Institute, Carey NC).

Intra-family heterogeneity of repeats

To determine the level of mutational divergence between repeat copies within each cluster, we conducted a BLASTn all against all reads following recommendations by Kelly et al. [28]. Self-hit blast results were removed from the output files. For each cluster, a histogram of the pairwise percentage identity scores (between 0.9 and 1.0) was constructed using R [33]. Clusters with a bias towards high read similarity reflect sequence conservation and possible recent proliferation; whereas clusters for which the majority of reads have lower similarity scores reflect older families that exhibit mutational divergence. To automatically analyse patterns in the similarity histograms, two regression models (linear: Y = a + bX; and quadratic: Y = a + bX + cX²) were fitted to the histogram of each cluster using a custom R script. Regression coefficients were extracted from the best-fitting model based on the Bayesian Information Criterion [34]. Based on these regression coefficients, we classified individual clusters into six categories: 1) positive linear regression (linear model has lowest BIC and b > 0.001), 2) absence of linear relation (linear model has lowest BIC and −0.001 < b < 0.001), 3) negative linear regression (linear model has lowest BIC and b < −0.001), 4) positive quadratic function (quadratic model has lowest BIC and curve opens upward, b² > 0), 5) negative quadratic function (quadratic model has lowest BIC, curve opens downward b² < 0 and optimum is below 99 % pairwise identity) and 6) negative quadratic function with high optimum (quadratic model has lowest BIC, curve opens downward b² < 0 and optimum is above 99 % pairwise identity). Based on this classification, we identified TE clusters showing high sequence similarity of copies as a proxy for recent activity and proliferation and we distinguished them from low similarity clusters representing older copies accumulating mutations [28, 35]. We interpreted categories 1, 4 and 6 as evidence for young and dynamic TEs, and 2, 3, 5 as older and well-silenced TEs. This regression approach is a simple heuristic method to distinguish between these categories and is not intended to detect statistically significant differences between clusters.

TE family genomic content and level of transcription

From a previous study, we developed Illumina transcriptomic resources for these same T. officinale Macranthoides accessions [26]. For each accession, contigs assembled through the RepeatExplorer pipeline were used as a new TE database to which genomic and transcriptomic Illumina reads were mapped. This analysis enabled exploration of the most abundant TE family transcription activity at the cluster level. Both genomic and transcriptomic single reads were aligned and read abundance estimated per cluster using the Trinity pipeline, including Bowtie2 and RSEM methods [36–38]. A quality score was given for each alignment and the best one was kept for abundance estimation. Log ratios of FPKM normalized values (as described by [39]) were used to quantify transcription for each cluster. For each accession, we used simple linear regression to test for a relationship across all clusters between genomic abundance and expression.

Methylation profiles of the T. officinale accessions

To assess differential methylation levels in TE sequences among accessions, we used a recently developed reduced representation bisulfite sequencing (RRBS) protocol that does not require a reference genome [40]. Individuals originating from the same mother plants as the ones used for genomic sequencing were grown in a common greenhouse environment. DNA was extracted from leaf tissue as described above. For each accession, following van Gurp et al. [40], reduced genomic libraries were generated by PstI digestion, size selection and bisulfite treatment to convert unmethylated (but not methylated) cytosines before being pooled and sequenced on ¼ lane of Illumina Hi-Seq 2500.

Bioinformatics analysis for methylation calling followed van Gurp et al. [40]. Briefly, in-line barcodes present in the reads were stripped and a sample-specific read group tag was added to the read name. Overlapping paired-end reads were merged using PEAR [41]. Merged and non-merged reads were then mapped using BWA-Meth [42] to the pseudo-reference sequence of the T. officinale genome constructed previously by van Gurp et al. [40]. The first four adapter-derived nucleotides were clipped from the merged and non-merged reads. Watson and Crick specific bam files representing both DNA strands were created, identifying methylation polymorphisms as positions where (a proportion of) cytosines were converted, but guanines on the opposite strand were unmodified. Both Watson and Crick bam files were processed through Freebayes to identify polymorphic variants (with following parameters -F 0 -E 1 -C 0 -G 0 --haplotype-length 1 --report-all-haplotype-alleles --report-monomorphic --report-genotype-likelihood-max --haplotype-length 1) [43]. Combining variant calls from both Watson and Crick strands allowed for distinguishing genetic from epigenetic variation [44]. Resulting variants were merged into two files: one for single methylation polymorphisms, the other for genetic polymorphisms using custom-made scripts freely available in github.

The methylation output file included information on location and level of methylation for each cytosine position called. Methylation proportions were calculated as number of reads with methylated position divided by total number of reads matching that same position. Following analyses were based on a subset of highly confident positions called in all 5 accessions with a minimum of 50x coverage. Out of the 40,456 cytosines in the RRBS data, 4733 were found in CG context. Those positions were occurring on 623 RRBS loci. Because methylation in CG context is highly bimodal (Additional file 2) with usually very high or very low methylation proportion [45] we collapsed CG quantitative methylation calls in three qualitative levels: low methylation (<25 %), intermediate methylation (25–75 %) and high methylation (>75 %). We defined a “true” methylation polymorphism (differentially methylated position or DMP) if the five accessions showed variation ranging from low to high methylation levels. However, a larger number of positions showed variation between accessions from low to intermediate methylation or from intermediate to high levels of methylation and are therefore included in the analysis. As multiple DMPs can be clustered on the same locus, we also calculated the average level of methylation across all CG positions within each RRBS locus, to identify Differentially Methylated Regions (DMRs). DMRs were called using the same qualitative thresholds as for DMPs. Using BLAST, homologous sequences corresponding to all filtered loci and the DMR subset were retrieved from the genomic TE database, reference transcriptome and subsequent RNA-Seq data [26]. We were thus able to annotate each loci as TE, functional transcripts or unknown. For each locus annotated as TE, transcript level was retrieved in all accessions.

The repetitive landscape of T. officinale apomictic genomes

In the five accessions, 96 to 118 clusters were identified, each representing ≥ 0.2 % of input nuclear reads (Table 1). Total lengths and genome proportions for each TE family in the different accessions are shown in Table 1 and Fig. 1. The genomic proportions of highly repetitive DNA varied among accessions, ranging from 27.4 to 38.6 %. Except for DIRS, Penelope retrotransposons and TIR transposons, all superfamilies of TE known in plants were retrieved in the datasets [1]. Most widely distributed were Ty1/Copia and Ty3/Gypsy retroelements with Maximus/SIRE and Chromovirus lineages respectively, hAT transposons, Tandem repeats and ribosomal DNA. Only few LINE elements and other transposon families were found among the five accessions.

Annotations	Macranthoides 11				Macranthoides 8				Macranthoides 12				Macranthoides 13				Macranthoides 3
	Nb clusters	Nb clusters in model 1-4-6	Size (Mbp)	% TE	Nb clusters	Nb clusters in model 1-4-6	Size (Mbp)	% TE	Nb clusters	Nb clusters in model 1-4-6	Size (Mbp)	% TE	Nb clusters	Nb clusters in model 1-4-6	Size (Mbp)	% TE	Nb clusters	Nb clusters in model 1-4-6	Size (Mbp)	% TE
Retrotransposons	76	23	193.42	22.36	74	20	185.16	21.41	88	40	272.43	31.50	78	33	231.15	26.72	88	37	243.70	28.17
LTR Elements Ty1/Copia	38	15	93.59	10.82	37	12	84.05	9.72	36	14	97.63	11.29	36	14	99.99	11.56	38	16	97.32	11.25
LTR Elements Ty3/Gypsy	31	6	85.86	9.93	28	8	78.05	9.02	46	22	152.89	17.67	32	14	104.24	12.05	36	14	103.73	11.99
Unclassified LTR	4	0	8.96	1.04	5	0	13.44	1.55	5	4	18.76	2.17	10	5	26.92	3.11	5	1	11.55	1.33
LINE	1	1	1.77	0.20	0	0	0.00	0.00	0	0	0.00	0.00	0	0	0.00	0.00	2	0	4.52	0.52
ParaRetrovirus (Caulimovirus)	1	1	1.13	0.13	0	0	0.00	0.00	0	0	0.00	0.00	0	0	0.00	0.00	0	0	0.00	0.00
Unclassified retrotransposons	1	0	2.10	0.24	4	0	9.62	1.11	1	0	3.15	0.36	0	0	0.00	0.00	7	6	26.58	3.07
DNA transposons	5	1	10.67	1.23	6	3	23.55	2.72	8	3	23.16	2.68	14	6	41.12	4.75	12	5	32.14	3.72
Unclassified repeats	12	4	26.93	3.11	14	8	31.51	3.64	8	5	16.50	1.91	8	5	32.87	3.80	10	7	24.25	2.80
Tandem repeats	1	1	1.13	0.13	0	0	0.00	0.00	1	1	1.65	0.19	2	2	2.84	0.33	1	1	1.45	0.17
Ribosomal DNA	2	2	4.46	0.52	2	2	3.66	0.42	9	9	19.98	2.31	9	9	18.45	2.13	7	7	14.65	1.69
TOTAL	96	33	236.60	27.35	96	33	243.87	28.19	113	59	333.72	38.58	111	55	326.44	37.74	118	57	316.19	36.55

Number of clusters, genomic abundance and proportions of most abundant transposable elements classified according to Wicker et al. [1] in the five accessions of T. officinale apomictic lineage. Number of clusters showing high intra-family homogeneity are also indicated

The T. officinale genome is dominated by Class I LTR retrotransposons ranging from 21.4 % (185.2 Mb) to 31.5 % (272.4 Mb) depending on the accession. Ty1/Copia content was more constant among accessions (ranging from 9.8 to 11.7 %) than Ty3/Gypsy retroelements (9.0–17.7 %). Compared with Class I LTR retrotransposons, non-LTR retrotransposons and DNA transposons were relatively infrequent. Class II transposons only represented 1.2 % (10.7 Mb) to 4.8 % (41.1 Mb) of the genome. In all accessions, total read lengths of the top 10 clusters did not represent more than 8 % of the genome.

Abundance of ribosomal DNA was variable among accessions ranging from 0.4 % (3.7 Mb) to 2.3 % (20.0 Mb). Tandem repeats including satellites were not found in the most abundant clusters of all accessions representing a low-fraction of the repetitive compartment. Proportions of unclassified repeats account for 1.9–3.8 % of the genome and could include additional repeat types.

From histograms of pairwise sequence similarity between reads, we can infer aspects of relative age and dynamics of TEs. Overall in the clonal lineage, between 34.4 and 51.3 % of clusters showed evidence of recent TE proliferation as indicated by high similarity between reads. All classes and most superfamilies of TEs exhibited to some extent, proliferation of copies (Table 1). Other clusters displayed evidence for short and recent amplification, rapidly followed by divergence of copies as demonstrated by the largest number of clusters in model 5 (Additional file 3). Peaks of similarity scores were reached around 95.3–95.5 % for accessions 8 and 11 and at 95.7 % for accessions 3, 12 and 13 (Additional file 3).

Genomic differences are driven mostly by few accession-specific repeat superfamilies

The abundance of the different repeat annotations, as estimated by read counts, differed between the five accessions at all levels of annotation (order, superfamily and family: likelihood ratio Chi-square test, p <0.0001, Additional file 4). The first clade including accessions 11 and 8, showed less overall TE content than the second clade comprising accessions 3, 12 and 13 (Fig. 1). This last clade mostly differed from the previous by showing a higher number of Ty3/Gypsy Chromovirus elements and to a lesser extent higher proportions of transposons and ribosomal DNA copies. This last clade was also more dynamic with 47.8–52.2 % (compared to 34.3 % in clade 1) of repeat showing evidence of ongoing TE proliferation (Table 1, Additional file 5).

Within the second clade (accessions 3, 12 and 13), differences observed at the superfamily level between accessions were due to a small number of specific TE copies. Accession 12 has more Ty3/Gypsy Chromovirus sequences (130.0 Mb accounting for 15.0 % of the genome) which were showing recent TE amplification (16 out of 36 clusters). In accessions 3 and 13, main differences were driven by DNA transposons increasing total repeat length (32.1–41.1 Mb compared to 23.2 Mb in accession 12, Table 1) and rate of ongoing proliferation (41.7–42.8 % of DNA transposon clusters compared to 37.5 % in accession 12, Table 1). Most dynamic annotated DNA transposons included hAT, Mutator and Helitron superfamilies (Additional file 5).

Relationship between TE content and transcriptional activity

To look at the relationship between TE genomic content and transcriptional activity, we used the newly made TE database to map previously generated RNA-Seq data [26]. Regression analyses showed weak but positive relationships between RNA and DNA content for several accessions which were statistically significant for accessions 8 and 11 (Fig. 2). Conversely, no relationship was observed in accessions 12 and 13. The absence of a positive trend in these two accessions is at least partially caused by a few TE families that are low-abundant at the DNA level but that are very highly transcribed (e.g. two LTR elements: Copia Ale II and Copia Tork and one DNA transposon: Mutator-like Transposase genes). However, TE families with highly homogeneous copies do not exhibit higher transcriptomic activity (Additional file 6).

The repetitive fraction of T. officinale is heterogeneous, young and dynamic

To survive, apomictic lineages need to control proliferation of most active TEs. A large part of TEs show peaks in sequence similarity scores between copies well below 100 %; as expected for TE families whose activities are under host genome control. In T. officinale, 50–70 % of the clusters are showing patterns of pairwise similarities that are consistent with relatively recent proliferation followed by rapid degeneration of copies. This specific pattern was observed previously in the genome of Fritillaria where copies had amplified simultaneously accompanied by low rates of deletion creating an ancient TE fraction that degenerated and diverged overtime [28].

On the other hand, over 1/3 to half of the TE families within the apomictic lineage appear to have highly conserved intra-family sequences reflecting their relatively young age and recent transposition activity. A previous study looking at a small number of retroelement sequences in four young asexual species including T. officinale came to the same conclusion [47]. Presence of highly conserved sequences can be used as a proxy to detect recent or still ongoing proliferation of TEs in the host genome as these TEs have not degenerated by mutation accumulation and have conserved the possibility of being active and inserted in a new location. In cotton, highly similar young LTR elements found only on chromosome 1 exhibit the same pattern and active transcription of these young LTRs generates siRNAs to maintain the silenced state of these TEs [35]. Due to such silencing, young and conserved Gossypium elements are not necessarily active and transposing. In the same way, Taraxacum TE families with highly homogeneous copies do not automatically exhibit higher transcriptional activity in leaf tissue. Activity and transposition of young TE could be restricted to earlier developmental stages and seed development; therefore, not detectable in the leaf transcriptome [48]. Still, in two of the analysed accessions, a few low-abundant TE families exhibit high transcriptomic levels. We can speculate that they may be escaping host genome control; however, even when these TEs also showed variable methylation between the accessions (see Table 2), there was no clear relation between methylation level and transcription level. Alternatively, transcription of repetitive sequences can be part of the silencing machinery targeting TEs via transcript degradation into small RNAs [35, 49].

TE differences within one clonal lineage are governed by few repeat families

To characterize the potential actors of intra-clonal genomic divergence within one apomictic lineage, we compared five accessions of T. officinale genotype Macranthoides at the genomic, transcriptomic and epigenetic levels. We detected differences in TE genomic content but also in transcription and methylation levels which could have associated effects on nearby gene regulation and expression.