To obtain putative TER loci from species for which whole genome sequence was available, a bioinformatics approach was employed. BLAST searches were queried with the portion of the A. thaliana RAD52 cDNA that overlaps the TER1 locus and the first 100 base pairs (bp) of the TAD3 predicted cDNA adjacent to the TER2 locus. Nucleotide BLAST searches against the A. lyrata whole genome sequence available through the JGI web portal were performed with default settings to identify orthologs of both genes. Using GenBank, we queried BAC sequences in the nucleotide collection specifying B. rapa as the target organism. When the Schrenkiella parvula (Eutrema parvulum) genome sequence was completed and deposited on GenBank, nucleotide BLAST was performed against its genome using the same strategy. Reciprocal nucleotide BLASTs were performed for all returned sequences to determine orthology with the query sequences.

Nucleotide BLAST searches against Brassica rapa, Capsella rubella, and Eutrema salsugineum (formerly Thellungiella halophila) whole genome sequences employed the standard settings in Phytozome. Searches were performed as described above, using RAD52 or TAD3 cDNA as query, and reciprocal BLASTs were performed on all returned sequences. In addition to these publicly available data, additional BLAST results were provided from the emerging genome sequence of Aethionema arabicum pre-publication (http://vegi.cs.mcgill.ca/).

To obtain putative TER loci from species lacking whole genome data, a degenerate PCR approach was employed. Primer design took advantage of the fact that coding regions flanking the two TER loci in A. thaliana are relatively well conserved. First we identified the RAD52 and TAD3 orthologs from A. lyrata and B. rapa using available genomic resources and aligned each gene with its respective A. thaliana copy. Using highly conserved coding regions within these genes, a series of modestly degenerate primers were designed for hiTAIL PCR in PrimaClade (Gadberry et al., 2005), and subsequently amplified regions adjacent to each gene for 10 species lacking genomic data: A. arenosa, A. neglecta, Boechera platysperma, Cardamine hirsuta, Cardamine rhomboidea, Crucihimalaya lasiocarpa, Dimorphocarpa wislizenii, Lepidium draba, Olimarabidopsis pumila, Turritis glabra. In each case, all robustly amplified bands >150 bp produced during hiTAIL PCR were gel extracted, cloned (TOPO-TA, Invitrogen Inc.), and sequenced. A minimum of two clones were sequenced for each PCR fragment to increase the chances of finding alternative loci. For each species resulting sequence data were used in BLAST searches against the A. thaliana genome to determine orthology of the coding portions of the amplified fragments. hiTAIL PCR reactions amplified a single locus in all species (see results) regardless of the band sequenced. Sequences were aligned, and degenerate primers were used to amplify the entire locus using standard PCR methods. Degenerate primers for standard PCR were placed in the second exon of RAD52 and the third exon of TAD3.

The PCR and BLAST screens were designed to return sequences orthologous to RAD52 and TAD3, along with adjacent sequence, from all sampled species. The sequences were aligned to determine whether adjacent regions contained putative TER loci. First, the RAD52 and TAD3 orthologous coding regions were aligned by eye from all sampled species with the A. thaliana genes using the visual alignment interface of the computer software MacClade (Maddison and Maddison, 2005). With these genes as anchors, the regions adjacent to them were aligned. The TER telomere template domain served as a landmark to aid alignment. For species lacking a template, the alignment program MUSCLE (Edgar, 2004) was used with default parameter settings to test if sequence similarity not evident under manual alignment was present.

To determine the evolutionary history of the TER duplication, we inferred a maximum likelihood phylogeny of all the recovered loci and TER1 and TER2 from Arabidopsis. Likelihood searches were performed in RAxML version 7.2.3 (Stamatakis, 2006) using the GTRGAMMA option from an alignment of 330 nt. Support for nodes in the tree was assessed using 100 bootstrap replicates and implemented in RAxML. The gene tree we recovered agreed with the accepted organismal phylogeny of the family (Beilstein et al., 2006, 2008), and thus we integrated the duplication history inferred from this gene tree with the organismal phylogeny. To test scenarios of evolution, we used the alignment of all 16 loci in the program MacClade, which reconstructs the evolution of each site in the alignment under a parsimony framework by inferring ancestral character states throughout a given tree. For these analyses, we chose the well-resolved organismal tree, which did not differ significantly from the gene tree.

Arabidopsis thaliana TER1 and TER2 are encoded on different chromosomes. Thus, our initial bioinformatics screens were designed to identify RAD52 orthologs in A. lyrata and B. rapa. For both species, a single locus with sequence similarity to A. thaliana RAD52 was retrieved; reciprocal BLASTs returned A. thaliana RAD52. Next A. thaliana TAD3 was queried via BLAST and a single locus with similarity to A. thaliana TAD3 was recovered. Surprisingly, these BLAST results also revealed that the two genes are located on the same scaffold in A. lyrata and the same BAC in B. rapa, although they are encoded on opposite strands. In addition, the start codons are separated by a relatively short intergenic region of 1,711 bp (A. lyrata) and 1,120 bp (B. rapa ; Figure 3).

RAD52 and TAD3 are located on the same assembled genomic scaffold in Capsella rubella, Eutrema salsugineum, Schrenkiella parvula, and Aethionema arabicum, whose genome sequences were released during the course of the project or were obtained from collaborators. Sequence data from A. arabicum is particularly informative because the genus Aethionema is sister to all other Brassicaceae. As for A. lyrata and B. rapa, BLAST results from these species returned a single locus with sequence similarity to either RAD52 or TAD3. In each case, an intergenic region of varying length separated the two genes: 607 bp (A. arabicum), 1567 bp (C. rubella), 895 bp (E. salsugineum), and 635 bp (S. parvula; Figure 3).

Both hiTAIL and standard PCR screens indicated that RAD52 and TAD3 are adjacent in Arabidopsis arenosa, A. neglecta, Camelina hispida, Cardamine hirsuta, C. rhomboidea, Crucihimalaya lasiocarpa, Dimorphocarpa wislizenii, and Lepidium draba. Initial sequence data were generated using hiTAIL PCR, which amplified multiple bands from degenerate primers. Fragments greater than 150 bp were gel-extracted and cloned. However for each species, a single locus with sequence similarity to the target region (either TAD3 or RAD52) was obtained from cloned PCR fragments, regardless of size. In C. hispida and C. lasiocarpa, hiTAIL PCR products were longer (~2000 bp) than in other species we sampled. Moreover, sequences generated from the cloned C. hispida and C. lasiocarpa PCR fragments included both the RAD52 and TAD3 orthologs (Figure 3). Subsequent PCR amplifications used standard PCR protocols with primers placed in the RAD52 and TAD3 genes. For all species we were able to generate PCR products between these two primers, indicating that the two genes are in the same genomic region (Figure 3). Sequence data from these cloned PCR products confirmed the region amplified contained RAD52 and TAD3. The length of the intergenic space for sequenced species was: 587 bp (A. arenosa), 538 bp (A. neglecta), 597 bp (C. hirsuta), 640 bp (C. rhomboidea), 564 bp (Crucihimalaya lasiocarpa), and 653 bp (D. wislizenii). For C. hispida, we did not sequence the entire length of the intergenic region, but the sequenced fragment was greater than 550 bp and contained a putative template domain.

In addition to identifying the RAD52 and TAD3 orthologs from the sampled species, we analyzed the sequence of each intergenic region to determine whether it shares greater identity with the A. thaliana TER1 or TER2 locus. In nine of the sampled species the majority of the intergenic space shares identity of at least 69% with the 545 bp domain upstream of the TER2 start position (Figure 3). In five species, less sequence similarity was detected. In A. lyrata, which has the longest intergenic region sampled, 445 bp of the 1711 bp intergenic region shares 77% sequence identity with the TER2 locus, while the 550 bp region upstream of the RAD52 start codon shares 74% sequence identity to a DNA tract located on A. thaliana chromosome 2. Roughly 25% (212 bp) of the E. salsugineum intergenic region shares 71% identity with TER2. For C. hispida, the region 86 bp upstream of the template domain shares 79 % identity with TER2. The C. rubella intergenic space shares little identity with TER2, rather a 134 bp segment of its intergenic space shares sequence identity with a DNA tract located several megabases away on A. thaliana chromosome 5. Finally, the A. arabicum intergenic space shared less than 25% sequence identity with either TER1 or TER2. Thus, the intergenic region in most of the sampled species is more similar to the region 5′ of A. thaliana TER2, rather than TER1.

In all species we sampled, RAD52 and TAD3 are located within a segment of 538 to 1711 bp and are encoded on alternate strands. In striking contrast, A. thaliana RAD52 and TAD3 are encoded on different chromosomes, and each is adjacent to one of the two TER genes. For the majority of species sampled, the intergenic region is highly similar to the stretch of DNA between the TAD3 start codon and the TER2 start, increasing confidence in the orthology of the region. Integration of this finding with the organismal phylogeny indicates that the event leading to the TER duplication occurred on the branch leading to A. thaliana (Figure 3), since even other species of the genus have only a single locus with sequence similarity to the two loci in A. thaliana. Altogether, these data indicate that the genomic region encoding the two A. thaliana TER occurs as a single locus in all other sampled Brassicaceae.

We aligned the TAD3 to RAD52 intergenic regions of the 14 sampled species. Using the alignment, we searched for retention of an A. thaliana like template domain upstream of the RAD52 start codon (Figure 1). An exact match to the A. thaliana template domain was recovered in six species: A. arenosa, A. neglecta, C. hispida, D. wislizenii, E. salsugineum, and S. parvula (Figures 1 and 3). Variations on the template domain occurred in five species. Nucleotide changes include the addition of a single extra adenine in B. rapa, increasing the stretch of A residues in the putative template domain from three to four nt (Figure 1). The same insertion is evident in C. hirsuta, C. rhomboidea, and A. arabicum. The putative template domains of C. rhomboidea and C. hirsuta are also altered by the insertion of an extra thymine at the 3′ end of the template region. C. hirsuta has a substitution of the terminal cytosine to thymine. Aethionema arabicum exhibits even more elaborate deterioration of the template domain with a two nt substitution at both the 5′ and 3′ ends. Only the core AAACCC repeat remains intact. Crucihimalaya lasiocarpa lacks even this core repeat, with a substitution at one of the central cytosines to thymine. Strikingly, we could find no evidence of a template or template-like domain in the entire length of the intergenic region from three species: A. lyrata, C. rubella, and L. draba. Moreover, no other loci with sequence similarity to the two A. thaliana TERs were recovered by BLAST or PCR approaches. We conclude that the template domain within TER is either evolving rapidly, or that alternative TER loci exist for several of the members of the Brassicaceae.

Phylogenetic analysis of the TER loci recovered a single likelihood tree (-ln 2564.5971; Figure 2). Model parameters from the most likely tree in RAxML included a shape parameter for the γ distribution (α = 1.78325) and a GTR transition matrix (a–c = 0.719862, a–g = 0.697629, a–t = 1.035392, c–g = 0.576520, c–t = 1.201706, g–t = 1.000000). Bootstrap support was generally low, however the Arabidopsis clade was resolved with 100% support, and A. arenosa, A. lyrata, and A. neglecta formed a well-supported monophyletic group within the Arabidopsis clade (92% bootstrap), indicating that the duplication that produced TER1 and TER2 is restricted to A. thaliana. Moreover, the recovered tree is consistent with the accepted organismal phylogeny of the family, and thus we used the organismal tree to infer evolutionary events. Integration of the putative TER loci with the organismal tree revealed that lineage-specific events affecting the locus occurred several times throughout the history of the family. The most parsimonious reconstruction of the evolution of the template indicates that the 5′-CTAAACCCTA-3′ encoding domain is ancestral in Brassicaceae (Figure 3). Species in both lineages I and II of the family have the canonical template, indicating they descended from a common ancestor with the template. Data for lineage III are currently missing, despite our efforts to obtain sequences for the group. In A. arabicum, there are numerous substitutions in the template domain. However, our data suggest that this is more likely a lineage-specific degradation of the template domain rather than the ancestral condition, especially since such modifications are relatively frequent, occurring in five of the 14 sampled species.

Substitutions in the template domain are independent evolutionary events. For example, insertion of adenine in the A stretch of the template occurred independently at least three times in the family: once each along the branch to A. arabicum and B. rapa, and once along the branch uniting C. hirsuta and C. rhomboidea. These latter two species share a thymine insertion at the 3′ end of the template, but C. hirsuta also has an independent thymine substitution at the 5′ end. Additional independent lineage-specific substitutions and insertions occurred in A. arabicum and C. lasiocarpa. Loss of the template domain occurs as an independent lineage-specific event in A. lyrata, C. rubella, and L. draba, since the closest sampled relative of each of these species has the canonical template. Interestingly, loss of the template domain does not appear to be correlated with a general loss of sequence identity in the intergenic region. For example, the L. draba intergenic region contains a 287 bp segment with 71% sequence identity to the TER2 locus, and in A. lyrata there is a 449 bp segment with 77% identity to TER2. Only in C. rubella is there both a loss of the template and a corresponding loss of sequence similarity. Moreover, losses are relatively common, occurring in three of the fourteen species (Figure 3). Losses are restricted to lineage I, however our sampling within this lineage is denser than in other groups of the family, suggesting denser sampling might reveal additional losses in other lineages. Finally, duplications of the TER locus are rare; the duplication producing two TER loci occurred along the branch to A. thaliana (Figure 3).

Interestingly, we found that all of the putative TER loci contain CR1 and CR2 overlapping with the RAD52 gene (Figure 3). This configuration is TER1-like in other Brassicaceae, indicating that the 529 nt non-conserved stretch that splits these conserved regions in TER2 is of recent origin. All other sampled species of Arabidopsis lack the 529 nt non-conserved stretch, hence the most recent common ancestor of the group is inferred to lack this stretch. The coalescence of A. thaliana and other Arabidopsis species occurred approximately 13 million years ago (Beilstein et al., 2010), and thus the insertion of this internal non-conserved sequence occurred sometime after this date. Similarly, the 235 nt non-conserved region 5′ of the template in TER1 is highly divergent from both TER2 and other Brassicaceae, supporting the conclusion that this region is of recent origin. In A. lyrata, the 550 bp stretch 5′ of the RAD52 start is also highly divergent from TER2 and other Brassicaceae. However, this stretch in A. lyrata is highly similar to a segment of A. thaliana chromosome 2, but not TER1. In A. arenosa and A. neglecta, the region 5′ of the template domain is highly similar to TER2. Hence, the A. thaliana TER1 non-conserved region and the A. lyrata non-homologous stretch are of independent and recent origin. A well-documented genome rearrangement and chromosome number reduction occurred on the branch leading to A. thaliana (Lysak et al., 2009), which may explain the duplication of the TER locus, as well as the origin of the non-conserved regions of both TER1 and TER2.