Recurrent evolution of heat-responsiveness in Brassicaceae COPIA elements
© The Author(s). 2016
Received: 5 April 2016
Accepted: 23 September 2016
Published: 11 October 2016
The mobilization of transposable elements (TEs) is suppressed by host genome defense mechanisms. Recent studies showed that the cis-regulatory region of Arabidopsis thaliana COPIA78/ONSEN retrotransposons contains heat-responsive elements (HREs), which cause their activation during heat stress. However, it remains unknown whether this is a common and potentially conserved trait and how it has evolved.
We show that ONSEN, COPIA37, TERESTRA, and ROMANIAT5 are the major families of heat-responsive TEs in A. lyrata and A. thaliana. Heat-responsiveness of COPIA families is correlated with the presence of putative high affinity heat shock factor binding HREs within their long terminal repeats in seven Brassicaceae species. The strong HRE of ONSEN is conserved over millions of years and has evolved by duplication of a proto-HRE sequence, which was already present early in the evolution of the Brassicaceae. However, HREs of most families are species-specific, and in Boechera stricta, the ONSEN HRE accumulated mutations and lost heat-responsiveness.
Gain of HREs does not always provide an ultimate selective advantage for TEs, but may increase the probability of their long-term survival during the co-evolution of hosts and genomic parasites.
KeywordsBrassicaceae COPIA Evolution Heat stress ONSEN
Transposable elements (TEs) are ubiquitous components of eukaryotic genomes. Their functions and roles range from DNA parasites, through regulators of gene transcription to facilitators of genome evolution (reviewed in [1, 2]). Together with other types of repeats, TEs comprise 10–80 % of plant genome content and specific families of long terminal repeat (LTR) retrotransposons can reach thousands of copies per genome [1, 3]. Plants evolved several layers of sophisticated epigenetic silencing mechanisms in order to suppress TE activity. Their transcripts are degraded by the post-transcriptional gene silencing (PTGS) pathway, which greatly reduces possible transposition events . In parallel, transcriptional gene silencing (TGS) stably silences TEs by deposition of DNA methylation via RNA directed DNA methylation (RdDM) mechanism (reviewed in [5, 6]). The repressed state is further stabilized by accumulation of specific histone modifications and faithfully transmitted in a DNA replication-dependent manner to the next generations. External or internal factors [7, 8] can lead to transient loss of silencing, but the epigenetic control will be re-established through tissue-specific RdDM activity . In addition to the nimble epigenetic silencing system, entire TEs can be physically removed from the host genome by deletion-biased homologous recombination processes .
In spite of the multi-layer amplification barriers, many TE families show signs of recent transpositions [11–13], suggesting that TEs occasionally escape epigenetic surveillance. There is increasing evidence that stress treatments affect chromatin structure and may lead to transposon activation (reviewed in [1, 14]). A possible mechanism was proposed based on the analysis of stress-induced TEs. LTRs of medicago cold-inducible repetitive element (MCIRE) retrotransposon contain a putative cold-responsive element (CRE) in alfalfa (Medicago sativa) . The CRE is specified by a conserved 5-bp core sequence (CCGAC) typical for C-repeat (CRT)/dehydration-responsive elements (DRE) that are recognized by cold-specific transcription factors (TFs) . LTRs of heat-responsive COPIA78/ONSEN (used as synonyms in this study) retrotransposon in Arabidopsis thaliana [7, 8, 17], contain a cluster of four nGAAn motifs forming a heat-responsive element (HRE) . During heat stress (HS), the ONSEN HRE is bound by heat shock factor A 2 (HSFA2), which triggers its transcriptional activity. This regulation is very specific and greatly independent of TGS control as the loss of decreased DNA methylation 1 (DDM1) in mutant plants did not trigger ONSEN transcriptional activation , in contrast to other typical LTR retrotransposons .
Presence of HRE and CRT/DRE motifs in ONSEN and MCIRE, respectively, suggested that the TEs’ response to stresses may be mediated by specific TF binding motifs. HREs were previously classified into four types based on their structure and, most likely, also activity . The strongest 4P HRE contains at least four adjacent nGAAn motifs and is bound by two HSFA2 trimers. The 3P HRE is bound by a single HSFA2 trimer and represents a moderately responsive HRE. In contrast, gap and step HREs with irregularly and more distantly spaced nGAAn motifs have on average lower HRE activity. Therefore, the HRE composition needs to be considered in order to define the strength of transcriptional response.
Here we identified multiple heat-responsive COPIA families in Arabidopsis lyrata and A. thaliana, two closely related species, using RNA sequencing (RNA-seq). Subsequently, we extended our analysis to five other Brassicaceae species and reconstructed putative HREs, their evolutionary history, and validated our predictions by transcriptional analysis after HS treatment.
Identification of heat-responsive TE families in A. thaliana and A. lyrata
To assess the extent of plant responses to HS, we monitored transcript levels from 32,793 A. lyrata and 32,678 A. thaliana protein-coding genes. This revealed significant upregulation (adjusted P <0.05; DESeq) of 21.8 % A. lyrata genes (n = 7156) and 18.9 % A. thaliana genes (n = 6165) after 6 h HS (Fig. 1c; Additional files 1 and 2). After recovery, we found only 2.9 % (n = 980) of genes still upregulated in A. lyrata and 0.6 % (n = 192) in A. thaliana. A. lyrata showed 21.3 % (n = 6992) downregulated genes after HS and 1.5 % (n = 491) after recovery (Fig. 1d). There were 17.3 % (n = 5650) significantly downregulated genes after HS and only 0.3 % (n = 89) after recovery in A. thaliana. Hence, HS treatment induced a similar degree of transient transcriptional changes in both species.
Because there is no publicly available A. lyrata TE annotation, we prepared custom-made catalogues of 53,089 A. lyrata and 17,009 A. thaliana repetitive elements (Additional files 3 and 4, respectively). Although the two species differed threefold in their TEs numbers, their spectra of TE families were similar (Additional file 5: Figure S1). The multi-copy nature of many TEs hinders RNA-seq analysis using standard protocols. Therefore, we developed the COMparative EXpression of TEs (COMEX) method, which allows quantification of transcripts derived from individual TE copies and effective removal of the RNA-seq reads mapping across TE families (see “Methods;” Additional files 6 and 7). We found 197 and 132 significantly (adjusted P <0.05; DESeq) upregulated TEs, representing 90 and 60 families (26 in common), after 6 h HS in A. lyrata and A. thaliana, respectively (Fig. 1e, f; Additional files 8 and 9). Comparing the major upregulated TE groups versus those in the whole genome revealed general under-representation of DNA transposons and HELITRONs and A. lyrata-specific under-representation of SINEs. In contrast, we found an over-representation of heat-responsive SADHU and LINE retrotransposons in A. lyrata, GYPSY elements in A. thaliana, and COPIA TEs in both species (Fig. 1g). Heat-responsive AlCOPIAs (n = 60; 100 %) comprised six families with at least two heat-inducible elements (Fig. 1h): AlCOPIA31 (n = 3; 3 %), AlCOPIA79 (n = 2; 3 %), AlCOPIA37 (n = 5; 11 %), AlCOPIA20 (n = 9; 14 %), AlONSEN (n = 19; 37 %), and a so far unknown family which we named TEMPERATURE RESPONSIVE TRANSPOSON (TERESTRA, n = 6; 10 %), as well as a bulk of single copies from different families (n = 16; 22 %). A. thaliana heat-responsive COPIAs (n = 34) were represented by six families with more than one heat-responsive TE. However, only AtONSEN (n = 8; 29 %) and AtCOPIA37 (n = 4; 12 %) were common between both species (Fig. 1i). A prominent A. thaliana-specific family was ROMANIAT5, comprising 12 % (n = 4) of all heat-responsive AtCOPIAs. After recovery, all TEs were re-silenced in A. thaliana and only five (AlCOPIA37, AlRE1, SADHU6-1, AlATN9_1, and AlLINE1_3A) showed increased transcript amounts in A. lyrata (Fig. 1e; Additional file 8). Surprisingly, ONSEN was fully silenced after two days of recovery, most likely owing to a shorter HS applied here compared to the previous study . The families representing at least 10 % of heat-responsive COPIA elements in each species were considered for further analysis (Fig. 1h).
Next, we tested whether heat-responsive COPIA families represent a particular COPIA clade. We reconstructed phylogeny of HS-responsive COPIA37, ONSEN, TERESTRA, ROMANIAT5, and seven HS-non-responsive COPIA families (19, 23, 26, 35, 45, 66, 75) based on their RT sequences (Fig. 1j; Additional file 5: Figure S2). The coding sequence was preferred over LTRs for the similarity analysis because this is strongly influenced by length of the input sequences, which may vary drastically in case of LTRs from different families. All heat-responsive families formed distinct and early separated branches, suggesting multiple independent origins of COPIA heat-responsiveness.
The structure and evolution of ONSEN heat-responsiveness
Copy numbers of elements within analyzed COPIA families in Brassicaceae species
We found 55 COPIA78 TEs in A. lyrata. Ten are full-length elements and 45 are fragments, either solo LTRs or incomplete according to the gaps in the genome assembly (Table 1; Additional file 5: Figure S3 and Table S2). In total, 15 copies contain at least one putative HRE with three or more adjacent (≤5 bp) nGAAn motifs. Remarkably, a high number of AlONSENs carry HREs identical to A. thaliana copies (Fig. 2a; Additional file 5: Figure S3). AlONSEN 2 and 8 have A. thaliana-like gap HREs; the 4P type is present in AlONSEN 10 and both co-occur in AlONSEN 6, 7, 9, 11, 14, 15, 17. In addition, we observed putative low efficiency gap/step HREs substituting the 4P HRE in AlONSEN 1, 4, 5, 12, and 13. All AlONSENs with predicted 4P HREs were upregulated after 6 h HS (Additional file 8). AlONSEN 3 and 16 were also found upregulated although they did not contain putative HREs. This was most likely caused by ambiguity in RNA-seq analysis, as 100 % of the reads mapping to these elements were multiply mapping to other ONSENs. Hence, there is a high correlation between the predicted HREs and RNA-seq results.
Conservation of the most frequent HRE haplotype between the Arabidopsis species raised the question about the evolutionary history of ONSEN heat-responsiveness in the Brassicaceae. Therefore, we searched for COPIA78 elements in whole-genome assemblies of Boechera stricta v1.2, Brassica rapa FPsc v1.3 (both JGI; Phytozome), Capsella rubella , Eutrema salsugineum , and low coverage assembly of Ballantinia antipoda (Vu, Finke, and Pecinka; unpublished data) using genome-wide BLAST searches. We confirmed the absence of COPIA78 in Capsella , but found at least one ONSEN copy in all other species (Table 1). RT nucleotide sequence identity was >80 % (Additional file 5: Figure S4), fitting previously proposed criteria for a single TE family . The LTR identity was lower (typically <70 %) due to the presence of insertions and deletions and decreased with phylogenetic distance (Additional file 5: Figure S5). Nevertheless, it allowed us reconstructing putative HREs. None of the other species contained the A. thaliana-like gap HRE (Fig. 2a; Additional file 5: Figure S3). However, there was a perfectly conserved 4P HREs in two out of three ONSENs in B. antipoda. ONSENs of other species either did not contain any HREs (B. stricta) or they represented only lower efficiency types and were non-homologous to the Arabidopsis HREs (Fig. 2a). To challenge the predicted HREs, we grew all species in vitro and quantified ONSEN transcript levels after 6 h and 12 h of HS (Fig. 2b). In agreement with RNA-seq results, we found massive 884–976-fold upregulation in A. lyrata and A. thaliana. There was also high (185-fold) upregulation in B. antipoda containing the putative 4P HRE, but lacking an additional gap HRE (Fig. 2a, b). In contrast, B. stricta, B. rapa, and E. salsugineum predicted to have no or only low efficiency HREs did not show strongly increased ONSEN transcript amounts (Fig. 2b).
To test whether HREs represent a major cis-regulatory element in ONSEN LTRs, we performed phylogenetic shadowing of the LTR consensus sequences (Fig. 2c). Although the 4P HRE region was partially conserved, there are several other similarly conserved regions. The longest stretch of conserved LTR sequence comprises approximately the first 25–30 bp (Fig. 2c), which may be required for TE RT.
By anchoring the structural information on the Brassicaceae chalcone synthase-based phylogeny, we reconstructed the evolutionary trajectory of ONSEN HREs (Fig. 2d). The nTTCnnGAAn motif, which can be considered as the non-functional sequence preceding the 4P HRE (proto-HRE), is present in B. rapa and E. salsugineum (Fig. 2a; Additional file 5: Figure S3), suggesting that it existed already at the onset of Brassicaceae evolution. Later, proto-HRE became duplicated and instantly created the high affinity 4P HRE. Molecular dating of the split of the B. antipoda lineage  suggests that this motif was maintained over 6–9 million years of evolution. However, the 4P HRE was occasionally lost due to accumulation of the point mutations (B. stricta) or deletion of whole elements (Capsella).
Species-specific gain of HREs in COPIA37 and the novel family TERESTRA
We also identified TERESTRA as a new retrotransposon heat-responsive family. The A. lyrata genome contains six TERESTRA copies sharing 97 % similarity (Fig. 3d). BLAST searches using TERESTRA sequences revealed only local similarities to ONSEN LTRs and COPIA46 GAG and POL domains and no other significant hits. Therefore, we performed de novo TERESTRA analysis. Based on the order of GAG and POL, TERESTRA was unambiguously identified as Ty1/COPIA LTR-retrotransposon (Fig. 3c). Based on only 70 % similarity in an alignment of TERESTRA to COPIA46 and ONSEN elements (Fig. 3d), we defined TERESTRA as a novel COPIA family. The consensus length of the complete AlTERESTRA element was 5116 bp and the 5’ and the 3’ LTR were 529 and 536 bp long, respectively (Fig. 3c; Additional file 5: Figure S7). Sequence analyses of AlTERESTRAs revealed that all copies are full length, contain a tRNA primer binding site and a polypurine tract, suggesting their autonomy (Fig. 3c). TERESTRA LTRs are relatively A-T-rich (69 %) and the consensus sequence contained only a small number of cytosines in symmetrical contexts (CG = 5, CHG = 0; H = A, T or C), which resembles LTR nucleotide composition of ONSEN . TERESTRA was missing in the A. thaliana Col-0 genome. Therefore, we extended our search to 50 A. thaliana accessions by genotyping them with TERESTRA-specific primers (Additional file 5: Table S3). This screen also gave negative results and suggested the absence of TERESTRA in A. thaliana. However, we found TERESTRA TEs in Arabidopsis cebennensis (95 % identity; Additional file 5: Figure S8) and Arabidopsis halleri (91 % identity; Additional file 5: Figure S9) using the NCBI sequence database. Furthermore, there were TERESTRAs in B. stricta (n = 14), B. rapa (n = 2), and E. salsugineum (n = 6; Table 1), but not outside of the Brassicaceae.
All six AlTERESTRAs were heat-responsive (Fig. 1h; Additional files 8). Screening of AlTERESTRA LTRs for possible HREs revealed a cluster of six nGAAn motifs, which can assemble either two partially overlapping gap HREs or a 4P HRE (Fig. 3e; Additional file 5: Figure S10). Based on the high AlTERESTRA transcriptional heat response (Fig. 3f), we favor the latter possibility. Another species with high TERESTRA transcriptional activation after 6 h and 12 h HS was B. stricta (Fig. 3f). By BLAST we found 14 TERESTRA copies in the B. stricta genome (Table 1). Eleven copies among them contain a complex cluster of up to five adjacent nGAAn motifs in their 5’ LTRs (Fig. 3d; Additional file 5: Figure S10). According to a conservative approach, three nGAAn motifs within this cluster can putatively form a low affinity gap HRE, but high TERESTRA activation in B. stricta suggests that all five motifs can establish 4P HREs as compatible with a more relaxed prediction (Fig. 3e, f). Importantly, all predicted B. stricta HREs are at positions different from those in A. lyrata HREs (Fig. 3e), highlighting their species-specific evolution (Fig. 3g). In B. rapa, one TERESTRA copy carries a putative gap HRE and another one a step HRE (Fig. 3e). Out of six TERESTRAs in E. salsugineum, three had putative step and one also an additional gap HRE. However, none of the HREs found in B. rapa and E. salsugineum was homologous to Arabidopsis or Boechera HREs and their predicted low HSF binding efficiency was congruent with the absence of heat-responsiveness in the RT-qPCR experiments (Fig. 3f).
Decreased AlCOPIA37 and AlTERESTRA transcript amounts after 12 h versus 6 h HS (Fig. 3b, f) contrasted with continuous transcript accumulation for AlONSEN (Fig. 2b). We hypothesized that the failure to maintain high transcript level could be caused by the TGS. Due to a lack of A. lyrata TGS mutants, we used a pharmacological approach to interfere with TE silencing . We treated 14-day-old A. lyrata plants with 10 μM 3-deazaneplanocin A (DZNep) and 40 μM zebularine, including control plants without treatment. DZNep is an S-adenosylhomocysteine synthesis inhibitor, which blocks the production of SAM, the methyl group donor required for DNA and histone methylation. Zebularine is a cytidine analog leading to DNA de-methylation and loss of silencing from specific transposons [27–29]. After two days of drug treatment, plants were heat-stressed for 0, 6, and 12 h and the amount of transcript analyzed by RT-qPCR (Fig. 3h–j). Control DZNep and zebularine treatment increased COPIA37 and ONSEN transcript tenfold and fivefold, respectively (Additional file 5: Figure S11), suggesting that both TEs can be weakly activated by TGS interference also without HS treatment. TERESTRA was not activated by the drug treatment. A combination of HS with drug treatments had strong additive effects in all cases, except for zebularine and HS-treated COPIA37 (Fig. 3h–j). Both ONSEN and TERESTRA transcripts accumulated at much higher levels that were not decreasing at 12 h HS (Fig. 3h, i). The effect was generally stronger for DZNep and weaker for zebularine. This suggests that the heat-induced TE transcript accumulation is rapidly suppressed by epigenetic means, in particular for TEs carrying lower affinity binding HREs.
AtROMANIAT5 contributes to transcriptional regulation of APUM9 under HS
We found ROMANIAT5 elements in genomes of all species except for B. antipoda (Table 1). Putative HREs were present in at least some copies of ROMANIAT5 in all species except for C. rubella that contained only solo LTRs. There were step and gap HREs in A. lyrata, B. rapa, and B. stricta, 3P/gap HREs in E. salsugineum, and 3P HREs in A. thaliana (Fig. 4i; Additional file 5: Figure S13). The predicted HSF binding affinity of individual HREs correlated well with the amount of ROMANIAT5 transcripts found after 6 and 12 h of HS (Fig. 4j). The only exception was B. rapa, which showed 42-fold upregulation after 12 h HS but the analyzed copies carried at most only low affinity step HRE. This could be due to the presence of heat-responsive ROMANIAT5 copies in the part of the B. rapa genome that is not yet assembled.
Transpositions and insertions of TEs may lead to loss of gene functionality [32, 33]. Therefore, TEs activity and mobility are tightly controlled by epigenetic means throughout the entire plant development [5, 6]. On the other hand, new insertions contribute to genome evolution and regulation of gene transcription . Therefore, it was already suggested in the early days of transposon research that, under conditions when diversity of regulatory patterns in a population may provide a better basis for selection, limited transposon activation could be beneficial . However, how occasional TE expression is provoked and how control is regained later is still a matter of debate. There is a rapidly increasing number of reports showing transient TE activation under various stress conditions reviewed in . Therefore, it was hypothesized that stresses may open a window for transpositions. Here, we introduced the ONSEN (COPIA78) family as a model for understanding TE control and behavior under HS. ONSEN shows massive transcriptional upregulation upon HS in A. thaliana and new insertions in progenies of heat-stressed Pol IV mutant [7, 8, 17]. The molecular basis of ONSEN heat-responsiveness was puzzling until recently, when a typical HSFA2 TF binding HRE was identified in its cis-regulatory region . Presence of canonical TF binding motifs in TE promoters was described for D. melanogaster and M. truncatula [15, 35]. However, the frequency of such activation strategy among TEs was unknown.
We analyzed LTRs of A. thaliana and A. lyrata heat-responsive COPIA TEs ONSEN, COPIA37, TERESTRA, and ROMANIAT5 for putative HREs. A minimum of three adjacent (<5 bp) nGAAn motifs can form a basal HRE, whose activity will depend on their distance and the total number . Heat-responsive COPIAs featured the whole spectrum of HREs ranging from the 4P types in ONSEN and TERESTRA, through 3P types in COPIA37 and ROMANIAT5 to a dozen of variable gap and step HREs in all these families. By comparing predicted HREs with transcriptional data, we conclude that gap and step HREs are mostly not sufficient to trigger HS-induced TE upregulation. This is congruent with their proposed low HSF binding efficiency . Predicted 3P HREs correlated with up to a hundred-fold (COPIA37, ROMANIAT5) and 4P HREs with up to a thousand-fold (ONSEN, TERESTRA) transcript accumulation upon HS. This suggests a strong correlation between putative HREs and the transcriptional response of the TEs.
Previously it was shown that the TGS machinery antagonizes the TE activation [7, 17]. We found that the speed of re-silencing during or after HS depends on the HRE type. While ONSEN, with the strong 4P HRE, accumulated transcript during entire HS exposure, TEs carrying lower affinity HREs typically showed a maximum transcript amount at 6 h HS and lower levels at 12 h HS. This silencing can be reduced by treatment with DNA methylation inhibitors. Hence, stressed plants take active measures to prevent TE transpositions already during ongoing HS treatment. However, HS-induced TE activation must not always aim at transposition, but can be part of the plant regulome . In A. thaliana, we found that heat-responsive AtROMANIAT5-2 controlled transcription of the APUM9 gene located downstream of the element. As we did not observe any evidence for high amount of a read-through transcript from ROMANIAT5-2 towards APUM9, we hypothesize that this transcriptional activation may be mediated rather by a specific three-dimensional chromatin organization at this locus. APUM9 gene was previously shown to be under control by HDA6 and synergistically by MOM1 and RdDM pathways, but not DDM1 and MET1 [30, 31]. Therefore, AtROMANIAT5-2 may represent a domesticated transposon with fine-tuned HS-regulated activation, contributing to transcriptional control of APUM9.
To challenge the hypothesis that HREs could be beneficial for TE amplification (but not necessarily for the host genome stability), we reconstructed evolutionary trajectories for HREs of ONSEN, COPIA37, TERESTRA, and ROMANIAT5 in the Brassicaceae. ONSEN was not heat-responsive in the early separated lineages represented by B. rapa and E. salsugineum, because its LTRs contained only one half of the 4P HRE (proto-HRE), which does not constitute a functional HRE. The proto-HRE became duplicated approximately 6–9 millions of years ago  and directly formed the present days 4P HRE found in the genus Arabidopsis and in the Australian species B. antipoda. Hence, ONSEN 4P HRE represents an evolutionary conserved cis-regulatory element. However, it should be noted that there are several other similarly or even more conserved regions within the ONSEN LTR. Whether they represent other TF binding sites and/or enhancers remains currently unknown. Furthermore, the ONSEN example shows that even high affinity HREs do not allow a TE to overrule the host genome defense, because their heat-responsiveness was lost in B. stricta, and the whole family became vanished from the C. rubella genome. In TERESTRA, high affinity 4P HREs evolved independently at two different LTR regions in the closely related species A. lyrata and B. stricta, while 3P HREs of COPIA37 emerged multiple times from a common nTTCn-rich LTR region. In contrast to ONSEN, HREs of these families are evolutionary young and species-specific. Whether they will be evolutionary successful, is an open question, but we speculate this to be the case for A. lyrata TERESTRA, where all genomic copies are full length, carry strong HRE, and respond to heat.
At present it is unknown whether higher temperatures in southern latitudes lead to greater amplification of heat-responsive TEs in subtropical relative to temperate zones. Although this is possible, there are also several factors that may act against such correlation. First, southern populations may reduce effects of HS by adaptation and growth at favorable microclimatic and/or temporal conditions . Second, the genomes are subject to purification mechanisms and the higher transposition rate may be opposed by a greater frequency of TE removal . Indeed, HS was shown to increase frequency of DNA sequence removal by a single strand annealing type of homologous recombination in transgenic constructs structurally resembling a LTR retrotransposon [38, 39]. Therefore, the final number of stress responsive TEs per genome may be the result of multiple effects acting in a complex network.
TEs evolve cis-regulatory elements, such as HREs, rapidly and independently in many groups. This may represent a strategy to produce new copies, constantly challenging the host defense system by searching for potential weak points. Successful regulatory elements may become evolutionary conserved and spread by new TE insertions in a self-reinforcing loop. However, these copies will be silenced and frequently removed from the genome. Hence, stress-mediated TE activation is likely not an unequivocal and straightforward winning principle, but rather a necessary strategy to survive under the pressure of the host defense systems. It is also likely that the host genome can benefit to some extent, and in specific cases, from cis-regulatory elements spread by TEs.
Plant materials and growth conditions
We used: Arabidopis thaliana Col-0 and Silex , Arabidopsis lyrata subsp. lyrata MN47, Ballantinia antipoda, Boechera stricta ES9, Brassica rapa FPSc, Capsella rubella, and Eutrema salsugineum. Before standard HS experiments, A. thaliana and A. lyrata seeds were placed on wet soil, stratified for one week at 4 °C, and then grown in a growth chamber (Percival) at 21 °C during the day and 16 °C during the night (16 h light/8 h dark) until plants reached approximately the five-leaves stage. Subsequently, a part of the plants was transferred to 37 °C for 6 h. RNA samples for sequencing were collected from some of the stressed plants and the controls directly after stress. The remaining stressed plants were allowed to recover at control conditions and collected after 48 h. Later, HS and drug-treatment experiments were performed with in vitro grown plants. First, the seeds were surface-sterilized with 8 % sodium hypochlorite for 6 to 12 min, washed with copious amounts of sterile water, dried under sterile conditions, and spread on sterile ½ Murashige-Skoog medium. After one week of stratification at 4 °C, the Petri dishes with seeds were transferred to growth chamber with a long day regime (16 h light/8 h dark) and constant temperature of 21 °C. Plates with rosettes at the pre-bolting stage were then placed in another chamber with 37 °C for 6 h. For combined drug and heat treatments, A. lyrata plants were grown as described above, then transferred to plates with no inhibitor, 10 μM DZNep, or 40 μM zebularine (both Sigma-Aldrich) for 48 h and then exposed to 0, 6, or 12 h at 37 °C HS. Aerial plant tissues were harvested immediately after the stress, flash frozen in liquid nitrogen, and stored at –80 °C.
Seeds of the Silex reporter line were sown directly on potting soil and stratified at 4 °C for 48 h. The pots were then placed in a Percival CU-22 L chamber at 21 °C with 12 h light (140 mmol m−2s−1) and 12 h dark. When the plants turned 14 days old, the pots were placed at 6 °C under the same light conditions for 24 h. At this time, control plants were moved back to the 21 °C chamber while HS plants underwent 24 h HS at 37 °C with light conditions as before. Immediately after the HS treatment, all pots were placed again at 21 °C. Fluorescence pictures of control and HS plants were taken at 0 and after 1, 2, and 5 days of recovery. Fluorescence imaging was performed using an Aequoria dark box with a mounted ORCAII CCD camera (Hamamatsu, Japan).
Nucleic acids extraction, cDNA synthesis, and RT-qPCR
Genomic DNA was isolated using the Phytopure gDNA Kit (GE Healthcare). Total RNA was isolated with the RNeasy Plant Mini Kit (Qiagen) with an on column DNaseI (Roche) digestion or by the standard Trizol method with additional DNaseI (Thermo Scientific) digestion. cDNA was synthesized from 1 μg total RNA per sample using the Revert Aid H-Minus First Strand cDNA synthesis kit with the oligo-d(T) primer (all Thermo Scientific). For strand-specific RT, total RNA of 6 h HS A. thaliana plants was divided into five aliquots which were converted into cDNA using (1) oligo-d(T) primer, APUM9 (2) sense and (3) antisense transcript primer, and ROMANIAT5-2 (4) sense and (5) antisense transcript primers. RT-qPCR analysis was performed on three biological replicates with at least two technical replicates in a CFX384 instrument (BIO-RAD) using the SensiMix Plus SyBr Kit (PEQLAB). Expression values were calculated relative to control-treated samples using the standard curve method  and normalized using the glyceraldehyde-3-phosphate-dehydrogenase C2 (GAPC-2) or the UBC28 gene with a stable expression under mock, HS, and recovery conditions. Primers used in this study are listed in Additional file 5: Table S4.
One μg total RNA per sample with RIN >8.0 (Agilent Bioanalyzer 2100) was used to construct strand non-specific sequencing libraries with the Illumina TruSeq RNA Library Kit v2 according to the manufacturer’s instructions. Library quality was tested on a Bioanalyzer and high-quality libraries were subsequently sequenced in the 100 bp single-end read mode using a HiSeq 2500 sequencer (Illumina). Adaptor sequences and low quality bases were trimmed and low quality reads were filtered out with the FAST-X toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) using custom-made scripts. Subsequently, reads were mapped to the corresponding reference genome (TAIR 10 genome assembly or A. lyrata genome assembly v1.0) using tophat2  with default settings. On average, >15 million sequencing reads per library passed trimming and quality filtering. The numbers of reads mapping to specific genomic positions were retrieved using Qualimap and the latest A. thaliana genome annotation TAIR10 and A. lyrata genome annotation v2  for genes and custom-made repeat annotations for TEs. The TE data were further processed with COMEX (see below) and data for genes were analyzed directly using the DESeq package in R software [43, 44].
COMparative EXpression of transposable elements (COMEX)
Accurate quantification of TEs expression using short read sequences is hampered by high similarity of potentially many genomic copies. We developed a simple protocol called COMEX (https://github.com/bpietzenuk/COMEX) that partially overcomes this problem and allows analysis of TE transcription from RNA-seq data. Out of >10 million reads per average sequencing library, 0.12 % and 0.73 % high-quality mappable reads corresponded to TEs within our custom made A. thaliana and A. lyrata, respectively, TE annotations. This suggests that TE expression analysis using RNA-seq can be made more sensitive by high sequencing depth. The reads were processed via a shell-script that merges the pipeline as follows. First, the binary mapping.bam file is converted into a readable .sam file. Subsequently, ends are printed (ToPrint_end1.py) to the .sam file and mapping errors are removed (Selectnonrepeated1.py). In the following step (Selectmultiplymapped1.py), the output files for the uniquely mapping and the multiply mapping reads are created. The high-quality uniquely mapping (UM) TE reads were directly accepted for expression analysis. High-quality multiply mapping TE reads were analyzed to identify those providing usable information. We classified multi-mapping reads into two categories: (1) informative reads mapping to multiple members of the same TE family (Specifically Multiply Mapped – SMM); and (2) non-informative reads mapping across TE families (Non-specifically Multiply-Mapped – NMM) using the TE annotation gff-file. Reads of the second category were discarded (new_cases1.py). Afterwards, UM and SMM are merged into a single .sam-file and converted into a binary .bam-file. Subsequently, the output file of the COMEX2.0-pipeline (filename.output.final.bam) containing the number of SMM and UM reads from the same TE family was retrieved using a strand non-specific protocol in Qualimap. To avoid a bias by repeated counting of SMM reads, we used the proportional read count method that divides the power of a read by the number of mapped positions. This provided the number of reads per individual TE families and TEs, which were subjected to statistical analysis using the DESeq package in R software . To avoid considering potentially large number TEs with minimal transcriptional changes, which would be later difficult to validate experimentally, we considered only those which had at least 0.55 RPKM in one of the experimental time points.
In silico sequence analysis
Sequences of interest were extracted from corresponding TE annotation files using bedtools . LTR reconstruction was carried out in LTR-Finder  or manually by pairwise and multiple alignments of the 3'end to the 5'ends of TE annotated regions using MUSCLE or multalin with the DNA 5–0 comparison table option. Structural analysis and annotation of TERESTRA was performed using LTR Finder and blastx using NCBI non-redundant protein sequences database. LTR_Finder was used in both analyses with the threshold option set to 2.0 using the tRNA database of A. thaliana to predict PBS. The LTR length range was set from 100–3500 and the minimum LTR distance was set to 1000. Other parameters were left at default settings. Search for ONSEN sequences within genomes of various Brassicaceae was done using BLASTN within Phytozome 10 [47, 48]. Hits with a sequence identity of >70 % were extracted and manually investigated. Positive hits with a query coverage <70 % were analyzed manually for sequence similarity with Multialign using the DNA 5–0 comparison table option. The input ONSEN RT and LTR sequences are provided in Additional files 10 and 11, respectively.
To analyze the evolutionary distance of the Ty1/COPIA LTR-retroelements, multiple sequence alignments of the RT domains were performed using the genomic nucleotide sequences in MUSCLE . RT protein sequences used for construction of the network and the tree (Fig. 1j and Additional file 5: Figure S2, respectively) are provided in Additional file 12. The evolutionary history was inferred using the Neighbor-Joining method (Kimura-2-Parameter method) with 1000 bootstrap replicates. Positions containing missing data and gaps were removed (pairwise deletion option) leading to a total of 862 position in the final dataset. The tree was visualized as an unrooted tree. Phylogenetic network of genomic RT domain sequences from Ty1/COPIA LTR-retroelements was constructed using Neighbor-Net  within the splitstree 4.0 package [51, 52]. The phylogenetic distances were calculated by LogDet-pairwise genetic distances using LDDIST  with imputed missing matrix entries. Multiple sequence alignments of CHS genomic sequences were performed using MUSCLE . The CHS phylogeny was inferred using the Maximum Likelihood tree based on the Kimura-2-parameter model with 1000 bootstrap replicates. CHS sequences were retrieved from . All positions containing gaps and missing data were eliminated. There were a total of 1267 positions in the final dataset. All phylogenetic trees were constructed within MEGA 7 . Phylogenetic shadowing and analysis of motif conservation was performed with mVISTA [55, 56] using LTR consensus sequences of different species prepared in BioEdit , allowing fasta ambiguity codes for low conserved positions. Sequences were aligned using AVID . The cutoff was defined as ≥70 % conservation over a 20 bp sliding window with the minimal consensus of 7 bp relative to A. lyrata 5’ LTR sequence.
Short sequence reads were deposited in the NCBI GEO archive under accession number GSE69077.
We thank J. de Meaux, E. Schranz, S. Woody, and K.S. Schumaker for seed material. We are grateful to B. Eilts, R. Gentges, P. Pecinkova, and D. Hamacher for plant care and technical assistance, A. Abdelsamad for basic identification of TEs using RepeatMasker, J. Ali for preparation of the shell script, O. Mittelsten Scheid for critical reading of the manuscript, and M. Koornneef, V. Cavrak, K. Schneeberger, and J. Jimenez-Gomez for numerous discussions.
This work was supported by the German Research Foundation (DFG) grant AP1859-2 to A.P within the Science Priority Program 1529 Adaptomics and by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) fellowship BEX 10896/14-7 to C.M.
Availability of data and materials
RNA-sequencing data have been uploaded to the Gene Expression Omnibus and are available under accession numbers [GEO: GSE69077]. Lists of A. lyrata and A. thaliana TEs are provided in Additional files 3 and 4. Primers for PCR analyses are available in Additional file 5: Table S4. The COMEX method is described in detail and is available to download from https://github.com/bpietzenuk/COMEX; this software is licensed under GNU General Public License version 3. The seeds of plants used in this work can be requested from the corresponding author.
Experiments in this study were conceived by AP, EB, and AMJ. Experiments were performed by BP, CM, HG, NB, and AA. RNA sequencing was analyzed by BP and NB. Phylogenetic analyses were performed by BP. Transcript quantifications were done by CM, BP, and HG. Silex experiments were performed by HG. AP wrote the manuscript with contribution of all authors. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Grandbastien M-A. LTR retrotransposons, handy hitchhikers of plant regulation and stress response. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms. 2015;1849:403–16.View ArticleGoogle Scholar
- Lisch D. How important are transposons for plant evolution? Nat Rev Genet. 2013;14:49–61.View ArticlePubMedGoogle Scholar
- Vu GTH, Schmutzer T, Bull F, Cao HX, Fuchs J, Tran TD, et al. Comparative genome analysis reveals divergent genome size evolution in a carnivorous plant genus. Plant Genome. 2015;8:0021. doi:10.3835/plantgenome2015.04.0021.
- Mari-Ordonez A, Marchais A, Etcheverry M, Martin A, Colot V, Voinnet O. Reconstructing de novo silencing of an active plant retrotransposon. Nat Genet. 2013;45:1029–39.View ArticlePubMedGoogle Scholar
- Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11:204–20.View ArticlePubMedPubMed CentralGoogle Scholar
- Matzke MA, Mosher RA. Rna-directed DNA methylation: An epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15:394–408.View ArticlePubMedGoogle Scholar
- Pecinka A, Dinh HQ, Baubec T, Rosa M, Lettner N, Scheid OM. Epigenetic regulation of repetitive elements is attenuated by prolonged heat stress in Arabidopsis. Plant Cell. 2010;22:3118–29.View ArticlePubMedPubMed CentralGoogle Scholar
- Tittel-Elmer M, Bucher E, Broger L, Mathieu O, Paszkowski J, Vaillant I. Stress-induced activation of heterochromatic transcription. PLoS Genet. 2010;6:e1001175.View ArticlePubMedPubMed CentralGoogle Scholar
- Baubec T, Finke A, Mittelsten Scheid O, Pecinka A. Meristem-specific expression of epigenetic regulators safeguards transposon silencing in Arabidopsis. EMBO Rep. 2014;15:446–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Devos KM, Brown JKM, Bennetzen JL. Genome size reduction through illegitimate recombination counteracts genome expansion in Arabidopsis. Genome Res. 2002;12:1075–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Hu TT, Pattyn P, Bakker EG, Cao J, Cheng J-F, Clark RM, et al. The Arabidopsis lyrata genome sequence and the basis of rapid genome size change. Nat Genet. 2011;43:476–81.View ArticlePubMedPubMed CentralGoogle Scholar
- Seymour DK, Koenig D, Hagmann J, Becker C, Weigel D. Evolution of DNA methylation patterns in the Brassicaceae is driven by differences in genome organization. PLoS Genet. 2014;10:e1004785.View ArticlePubMedPubMed CentralGoogle Scholar
- Willing E-M, Rawat V, Mandáková T, Maumus F, James GV, Nordström KJV, et al. Genome expansion of Arabis alpina linked with retrotransposition and reduced symmetric DNA methylation. Nature Plants. 2015;1:14023.View ArticlePubMedGoogle Scholar
- Gutzat R, Mittelsten Scheid O. Epigenetic responses to stress: Triple defense? Current Opin Plant Biol. 2012;15:568–73.View ArticleGoogle Scholar
- Ivashuta S, Naumkina M, Gau M, Uchiyama K, Isobe S, Mizukami Y, et al. Genotype-dependent transcriptional activation of novel repetitive elements during cold acclimation of alfalfa (Medicago sativa). Plant J. 2002;31:615–27.View ArticlePubMedGoogle Scholar
- Nakashima K, Ito Y, Yamaguchi-Shinozaki K. Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 2009;149:88–95.View ArticlePubMedPubMed CentralGoogle Scholar
- Ito H, Gaubert H, Bucher E, Mirouze M, Vaillant I, Paszkowski J. An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress. Nature. 2011;472:115–9.View ArticlePubMedGoogle Scholar
- Cavrak VV, Lettner N, Jamge S, Kosarewicz A, Bayer LM, Mittelsten Scheid O. How a retrotransposon exploits the plant’s heat stress response for its activation. PLoS Genet. 2014;10:e1004115.View ArticlePubMedPubMed CentralGoogle Scholar
- Zemach A, Kim MY, Hsieh P-H, Coleman-Derr D, Eshed-Williams L, Thao K, et al. The Arabidopsis nucleosome remodeler DDM1 allows DNA methyltransferases to access H1-containing heterochromatin. Cell. 2013;153:193–205.View ArticlePubMedPubMed CentralGoogle Scholar
- Sakurai H, Enoki Y. Novel aspects of heat shock factors: DNA recognition, chromatin modulation and gene expression. FEBS J. 2010;277:4140–9.View ArticlePubMedGoogle Scholar
- Slotte T, Hazzouri KM, Agren JA, Koenig D, Maumus F, Guo Y-L, et al. The Capsella rubella genome and the genomic consequences of rapid mating system evolution. Nat Genet. 2013;45:831–5.View ArticlePubMedGoogle Scholar
- Yang R, Jarvis DJ, Chen H, Beilstein M, Grimwood J, Jenkins J, et al. The reference genome of the halophytic plant Eutrema salsugineum. Frontiers Plant Sci. 2013;4:46.Google Scholar
- Ito H, Yoshida T, Tsukahara S, Kawabe A. Evolution of the ONSEN retrotransposon family activated upon heat stress in Brassicaceae. Gene. 2013;518:256–61.View ArticlePubMedGoogle Scholar
- Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, et al. A unified classification system for eukaryotic transposable elements. Nat Rev Genet. 2007;8:973–82.View ArticlePubMedGoogle Scholar
- Mandakova T, Joly S, Krzywinski M, Mummenhoff K, Lysak MA. Fast diploidization in close mesopolyploid relatives of Arabidopsis. Plant Cell. 2010;22:2277–90.View ArticlePubMedPubMed CentralGoogle Scholar
- Pecinka A, Liu C-H. Drugs for plant chromosome and chromatin research. Cytogenet Genome Res. 2014;143:51–9.View ArticlePubMedGoogle Scholar
- Baubec T, Pecinka A, Rozhon W, Mittelsten SO. Effective, homogeneous and transient interference with cytosine methylation in plant genomic DNA by zebularine. Plant J. 2009;57:542–54.View ArticlePubMedGoogle Scholar
- Foerster AM, Dinh HQ, Sedman L, Wohlrab B, Mittelsten SO. Genetic rearrangements can modify chromatin features at epialleles. PLoS Genet. 2011;7:e1002331.View ArticlePubMedPubMed CentralGoogle Scholar
- Liu C-H, Finke A, Diaz M, Rozhon W, Poppenberger B, Baubec T, et al. Repair of DNA damage induced by the cytidine analog zebularine requires ATR and ATM in Arabidopsis. Plant Cell. 2015;27:1788–800.View ArticlePubMedPubMed CentralGoogle Scholar
- Yokthongwattana C, Bucher E, Čaikovski M, Vaillant I, Nicolet J, Scheid OM, et al. MOM1 and POL-IV/V interactions regulate the intensity and specificity of transcriptional gene silencing. EMBO J. 2010;29:340–51.View ArticlePubMedGoogle Scholar
- Hristova E, Fal K, Klemme L, Windels D, Bucher E. HISTONE DEACETYLASE 6 controls gene expression patterning and DNA methylation-independent euchromatic silencing. Plant Physiol. 2015;168:1298–308.View ArticlePubMedPubMed CentralGoogle Scholar
- Mirouze M, Reinders J, Bucher E, Nishimura T, Schneeberger K, Ossowski S, et al. Selective epigenetic control of retrotransposition in Arabidopsis. Nature. 2009;461:427–30.View ArticlePubMedGoogle Scholar
- Miura A, Yonebayashi S, Watanabe K, Toyama T, Shimada H, Kakutani T. Mobilization of transposons by a mutation abolishing full DNA methylation in Arabidopsis. Nature. 2001;411:212–4.View ArticlePubMedGoogle Scholar
- McClintock B. The significance of responses of the genome to challenge. Science. 1984;226:792–801.View ArticlePubMedGoogle Scholar
- Guerreiro MPG. What makes transposable elements move in the Drosophila genome? Heredity. 2012;108:461–8.View ArticleGoogle Scholar
- Enoki Y, Sakurai H. Diversity in DNA recognition by heat shock transcription factors (HSFS) from model organisms. FEBS Lett. 2011;585:1293–8.View ArticlePubMedGoogle Scholar
- Wilczek AM, Cooper MD, Korves TM, Schmitt J. Lagging adaptation to warming climate in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2014;111:7906–13.View ArticlePubMedPubMed CentralGoogle Scholar
- Boyko A, Filkowski J, Kovalchuk I. Homologous recombination in plants is temperature and day-length dependent. Mut Res. 2005;572:73–83.View ArticleGoogle Scholar
- Pecinka A, Rosa M, Schikora A, Berlinger M, Hirt H, Luschnig C, et al. Transgenerational stress memory is not a general response in Arabidopsis. PLoS One. 2009;4:e5202.View ArticlePubMedPubMed CentralGoogle Scholar
- Larionov A, Krause A, Miller W. A standard curve based method for relative real time pcr data processing. BMC Bioinformatics. 2005;6:62.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg S. Tophat2: Accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36.View ArticlePubMedPubMed CentralGoogle Scholar
- Rawat V, Abdelsamad A, Pietzenuk B, Seymour DK, Koenig D, Weigel D, et al. Improving the annotation of Arabidopsis lyrata using RNA-seq data. PLoS One. 2015;10:e0137391.View ArticlePubMedPubMed CentralGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.View ArticlePubMedPubMed CentralGoogle Scholar
- R Core Team. R: A language and environment for statistical computing. Vienna: Foundation for Statistical Computing; 2013. http://www.R-project.org/.Google Scholar
- Quinlan AR, Hall IM. Bedtools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.View ArticlePubMedPubMed CentralGoogle Scholar
- Xu Z, Wang H. Ltr_finder: An efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 2007;35:W265–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.View ArticlePubMedGoogle Scholar
- Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, et al. Phytozome: A comparative platform for green plant genomics. Nucleic Acids Res. 2012;40:D1178–86.View ArticlePubMedGoogle Scholar
- Edgar RC. Muscle: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Bryant D, Moulton V. Neighbor-net: An agglomerative method for the construction of phylogenetic networks. Mol Biol Evol. 2004;21:255–65.View ArticlePubMedGoogle Scholar
- Huson DH. Splitstree: Analyzing and visualizing evolutionary data. Bioinformatics. 1998;14:68–73.View ArticlePubMedGoogle Scholar
- Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.View ArticlePubMedGoogle Scholar
- Thollesson M. Lddist: A perl module for calculating logdet pair-wise distances for protein and nucleotide sequences. Bioinformatics. 2004;20:416–8.View ArticlePubMedGoogle Scholar
- Kumar S, Stecher G, Tamura K. Mega7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4.View ArticlePubMedGoogle Scholar
- Mayor C, Brudno M, Schwartz JR, Poliakov A, Rubin EM, Frazer KA, et al. Vista: Visualizing global DNA sequence alignments of arbitrary length. Bioinformatics. 2000;16:1046–7.View ArticlePubMedGoogle Scholar
- Frazer KA, Pachter L, Poliakov A, Rubin EM, Dubchak I. Vista: Computational tools for comparative genomics. Nucleic Acids Res. 2004;32:W273–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Hall TA. Bioedit: A user-friendly biological sequence alignment editor and analysis program for windows 95/98/nt. Nucleic Acids Symp Ser. 1999;41:95–8.Google Scholar
- Bray N, Dubchak I, Pachter L. Avid: A global alignment program. Genome Res. 2003;13:97–102.View ArticlePubMedPubMed CentralGoogle Scholar