Lectures Notes Strategic Management Transposable elements

Transposable elements (TEs) are ubiquitous in eukaryotic genomes and their

mobility impacts genome structure and function in myriad ways. Because of their

abundance, activity, and repetitive nature, the characterization and analysis of TEs

remains challenging, particularly from short-read sequencing projects. To overcome this

difficulty, we have developed a method that estimates TE copy number from short-read

sequences. To test the accuracy of our method, we first performed an in silico analysis of

the reference Sorghum bicolor genome, using both reference-based and de novo

approaches. The resulting TE copy number estimates were strikingly similar to the

annotated numbers. We then tested our method on real short read data by estimating TE

copy numbers in several accessions of S. bicolor and its close relative S. propinquum.

Both methods effectively identify and rank similar TE families from highest to lowest

abundance. We found that de novo characterization was effective at capturing qualitative

variation, but underestimated the abundance of some TE families, specifically families of

more ancient origin. In addition, interspecific reference-based mapping of S. propinquum

reads to the S. bicolor database failed to fully describe TE content in S. propinquum,

indicative of recent TE activity leading to changes in the respective repetitive landscapes

over very short evolutionary timescales. We conclude that reference-based analyses are

best suited for within-species comparisons, while de novo approaches are more reliable

for evolutionarily distant comparisons.

In plant genomes, the transposable elements (TE) community is both

quantitatively and qualitatively dynamic. TEs occupy a variable proportion of the

genome, and many studies have shown that their differential accumulation and deletion

strongly correlates with genome size (Bennetzen 2000; Bennetzen 2002; Tenaillon et al.

2010; Michael 2014). In some cases, rapid rates of amplification and decay can lead to

variable TE composition among closely related plant lineages. For instance, although

similar numbers of TE families occupy both the rice and maize genomes, transpositional

bursts of a few Long Terminal Repeat (LTR) retrotransposon families have inflated the

maize genome to six times that of rice (Baucom et al. 2009; Baucom et al. 2009).

Variation is also common among more closely related taxa, such as those that belong to

the same subfamily (SanMiguel and Bennetzen. 1998), and even among members of a

single genus. Interspecific comparisons performed in Gossypium and Oryza revealed that

proliferation of a small subset of TE families were responsible for the observed genome

size variation among species (Hawkins et al. 2006; Piegu et al. 2006). TE copy number

variation is not always due to the activity of just a few families, however, as has been

demonstrated by the proliferation of several different TE families in Zea and Arabidopsis

(Hollister et al. 2011; Hu et al. 2011; Tenaillon et al. 2011). Given that ongoing and

punctuated TE activity has been widely documented in plants, it is important to

understand the contribution of TEs to lineage-specific novelty, as insertional

polymorphisms may ultimately contribute to species diversification.

Although TEs are predominantly silenced, protecting the genome from rampant

insertional mutagenesis, the host suppression system can be circumvented in some

situations. TE activation has been reported following both hybridization and

polyploidization (Kashkush et al. 2002; Beaulieu et al. 2009; Parisod et al. 2009; Xu et al.

2009). An example is that of hybrid Sunflower species, in which the hybrid genomes are

50% larger (~1,130 Mb additional DNA) than either of the parental genomes (Ungerer et

al. 2006). This increase was attributed to recent proliferation (0.5 - 1 mya) of Ty3/Gypsy

elements (Ungerer et al. 2009). In addition, TEs are often activated in response to various

biotic and abiotic stresses such as infection, temperature, wounding and salinity (Mhiri et

al. 1997; Grandbastien et al. 1998; Takeda 1998; Pecinka 2010; Tittel-Elmer et al. 2010;

Fujino et al. 2011; Cavrak et al. 2014; Makarevitch et al. 2015; Finatto et al. 2015). Such

TE activity can generate numerous insertional polymorphisms, each with the potential for

functional consequences when inserted into genes or gene regions (Kashkush et al. 2003;

Kashkush and Khasdan 2007; Chu et al. 2011; Ito et al. 2011; Studer et al. 2011; Butelli

et al. 2012). Although TE activation might negatively affect genome integrity, it has been

suggested that species or populations that are prone to strong diversifying selection would

benefit from TE driven genome variability (Huang et al. 2008; Martin et al. 2009;

Fernandez et al. 2010; Lockton and Gaut 2010).

To date, analyses of TE copy number that could inform studies of activity and

accumulation have been hindered by the inability to analyze the repetitive fraction from

short-read sequencing projects. This difficulty is due to the inability to accurately

assemble short sequence reads that belong to repeats, as these reads often create assembly

gaps, incorrectly collapse onto a single chromosomal position, and/or map to multiple

locations in the genome, resulting in misassembled arrangements. Most modern

assemblers attempt to resolve these issues by employing alignment strategies that either

discard multiply mapping reads, report all possible mapping locations, or randomly map

reads to the position of best alignment; however, the random placement or all-together

removal of multiply mapping reads clearly prevents detailed analyses of repeat regions.

Additionally, although there are assemblers that are efficient enough to report all possible

mapping locations for repetitive sequences, they work best for high-coverage datasets

(>20x) or with longer sequence read lengths (Phillippy et al. 2008, Treangen and Salzberg

et al. 2012). These limitations pose a serious challenge to the characterization of repeat

content and to the determination of repeat copy numbers from NGS datasets.

Here, we demonstrate a method to evaluate TE composition and accurately estimate copy

number using Illumina short-read sequence data. We first tested the accuracy of our

method by performing an in silico analysis of the reference BTx623 Sorghum bicolor

genome, which resulted in copy number estimates that are strikingly mathematically

similar to the annotated TE copy numbers. Following this simulated analysis, we tested

our methods on real short read datasets by estimating copy numbers in several accessions

of S. bicolor and its close wild relative, S. propinquum. Copy number estimations were

performed using both reference-based and de novo methods for comparison. Both

methods rank the families in similar order from highest to lowest abundance; however,

the estimated copy numbers via de novo analysis differ from that of the reference-based

approach, and these differences correlate with the relative insertional timing of individual

TE families. Specifically, we find that de novo approaches are more effective in

estimating copy numbers for young TE families but tend to underestimate the abundance

of older families. In addition, interspecific reference-based mapping failed to fully

describe TE content in S. propinquum due to inefficient mapping of reads to the S. bicolor

database. We conclude that reference-based methods are best for estimating within

species variation whereas de novo approaches are more reliable for evolutionarily distant

comparisons.

Materials and methods

In silico development of method for copy number estimation

As the availability of a high-quality reference genome was required to develop our

approach, we focused on the genus Sorghum. The Sorghum bicolor reference genome

(BTx623, v1.0 Paterson et al. 2009) was downloaded from Phytozome v9.0, and full-

length Long Terminal Repeat (LTR)-retrotransposons were identified with LTRharvest

(Ellinghaus et al. 2008) using the default settings, except for the following: a motif for 5’

and 3’ LTRs as each LTR should begin and end with TG and CA nucleotides, minimum

and maximum length LTRs of 100 – 5,000 bp, and seed length set to 60 bp. Sequences

that were incorrectly identified as LTR-retrotransposon (false positives) were removed

from the output by performing a nucleotide BLAST against the Plant Genome and

Systems Biology (PGSB, formerly MIPS) Poaceae repeat database (Nussbaumer et al.

2013) with an e-value cutoff of 1e-10 and sequence identity of 80%. Sequences that did

not match to grass-specific LTR-retrotransposons were removed. An all-by-all BLAST (e-

value cutoff 1e-10 and at least 80% sequence similarity) was performed with the 5’LTRs

and the result was clustered into families using RepMiner (J. Estill, code available at

http://repminer.sourceforge.net, Baucom et al. 2009), and resulting clusters were

visualized as a network using the imaging program Cytoscape 3.0 (Shannon et al. 2003).

In addition, the percent sequence divergence between the 5’ and 3’ LTR of each element

was extracted from the LTRharvest output and used to determine the insertional timing of

each element in the BTx623 reference genome using a grassspecific substitution rate of

1.3 x 10-8 per site per year (Ma and Bennetzen, 2004).

RepeatMasker (employing default settings) was used to identify additional LTR

sequences that are not present as full-length sequences (solo LTRs) using the LTRharvest

output sequences as a database (Smit, Hubley and Green RepeatMasker Open-4.0.

20132015). All of the full-length elements identified by LTRharvest were first masked

using the feature maskfasta from Bedtools 2.17.0 package (Quinlan and Hall 2010), so as

not to count these LTRs twice. The RepeatMasker output was filtered for hits that were

less than 5% divergent from the database sequences and were at least 150 bp in length.

LTR sequences from both the LTRharvest and Repeatmasker analyses were combined

into a single dataset for further analysis. In addition, exemplars that represent the entire

population of sequences were selected from the extracted 5’ LTRs using affinity

propagation clustering, which were then used for de novo analysis (see below) (Frey and

Dueck Science 2007; Bodenhofer et al. 2011).

To devise an accurate method to estimate TE copy number from the short read data, we

tested the equation from Hawkins et al. (2006) on an in silico dataset generated from the

BTx623 reference using the short read simulator, DWGSIM v.0.1.11. The program was

run on ‘illumina’ mode to generate 7.5 million reads of 100 bp in length. We used only

single-end reads for our analysis in order to treat each read as a mathematically

independent sample. The simulated reads were mapped back to the reference genome

sequence using Bowtie2 under default settings (Langmead and Salzberg 2012). With the

known chromosomal positions for each identified element, the total number of reads that

are strictly and uniquely aligned (i.e., entire 100 bp read) within the first and last

nucleotide position of each identified 5’ LTR was extracted from the BAM alignment file

using the intersectBed and coverageBed tools of Bedtools 2.17.0.

The copy number of each element (n) was estimated using the following equation

(Hawkins et al. 2006):

n =  XobsN *(G −e)* (lt −21m + e)

Where Xobs is the total number of reads aligned to a given element (5' LTR), N is the

total number of reads used during mapping, lt is the length of target sequence (5’ LTR), m

is the overlap required to count a match to the target region (in this case, the entire length,

or 100 bp), e is the length of the sequence read, and G is genome size. Simply put, this

equation estimates how frequently a sequence of a particular length (the LTR) must be

present in a genome of a given size, if a proportion of random samples (reads that

match/total random samples) from that genome match the particular LTR. The greater the

proportion of reads that map to a particular LTR of a given size, therefore, the higher the

copy number estimate for that LTR in the genome. To evaluate possible sampling effects,

we repeated the in silico subsampling and statistical estimates for 100 independently

simulated datasets. Confidence intervals (95%) were calculated for all copy number

estimates as described in Hawkins et al (2006).

We also tested our method using a de novo assembly approach, implemented in

RepeatExplorer, for the same simulated short read sequence data used in the

referencebased approach described above (Novak et al. 2013). RepeatExplorer employs a

graphbased clustering method by quantifying the similarities between reads. The program

begins by filtering reads that pass a specific threshold (> 90% sequence similarity over

55% of the read length). Using these similarities, the program constructs a graph and

creates clusters from frequently connected reads that represent individual repetitive

families. Reads within each cluster are then assembled into contigs using CAP3 with an

overlap length cutoff of 50 bp (for 100 bp reads) and sequence identity of at least 80%

(Novak et al. 2010). To classify clusters into specific LTR-retrotransposon families,

BLAST was performed between the RepeatExplorer contigs from the largest clusters and

the 5’ LTR exemplars from the reference genome. After classifying the clusters into

families, copy numbers were estimated from the contigs present in each cluster

employing the same probability equation used in the reference-based quantification.

Copy numbers were estimated for the ten largest TE families for each accession.