Highly syntenic regions in the genomes of soybean, Medicago truncatula, and Arabidopsis thaliana

Identifying and mapping homologous contigs

There is a gap in all three species between synteny blocks 1a and 1b which we were unable to span. In soybean, these two sequence assemblies are genetically linked and located 2 cM apart on LG-G. In addition, synteny block 1b contained two gaps in M. truncatula (Figure 1), one a 90 kbp gap toward the bottom M. truncatula sequence assembly. This gap and its surrounding M. truncatula sequence (totaling about 175 kbp) correspond to an insertion/deletion in soybean just 25 kbp in size and containing two gene models without hits to Genbank's nonredundant (nr) database along with one repetitive sequence. There is also a gap of unknown size between the top and bottom sequence assemblies in synteny block 1b that could not be spanned with end sequenced BACs. In M. truncatula, the top sequence assembly maps to M. truncatula chromosome 4, while the bottom maps to chromosome 3 (Table 1). These two M. truncatula assemblies therefore appear to be unlinked, even though they show substantial synteny and are apparently both orthologous (see below) to a contiguous region in soybean.

M. truncatula sequences in synteny block 2 also map to chromosomes 3 and 4. One of the M. truncatula homoeologs in synteny block 2, Mt_2ii, maps 5–8 cM below the M. truncatula bottom sequence assembly in synteny block 1b (Figure 1, 2, Table 1). The other M. truncatula duplicate in synteny block 2, Mt_2i, maps to chromosome 4, more than 25 cM from the top sequence assembly in synteny block 1b (Figure 1, 2, Table 1).

Synteny

Perspectives on cyst nematode resistance genes

Despite the absence of an rhg1 homolog in the syntenic M. truncatula region examined here (synteny block 1b, Mt, top sequence assembly in Figure 1), it is clear that there is an rhg1 homolog elsewhere in M. truncatula. Though not present in any of the full-length BAC sequences that currently comprise >40% of M. truncatula's genespace [33], a homolog does exist on a M. truncatula BAC-end sequence (mth2-60m7), showing 78% amino acid identity to rhg1 over 279 amino acids (e-value = e-118). This percent identity is higher than that of the syntenic M. truncatula homologs compared to soybean's Rhg4 locus in synteny block 2 (~70%). BAC mth2-60m7 and the next three BAC-end hits all belong to the same region of M. truncatula FPC contig 949 [33], which tentatively maps to chromosome 5. This indicates that either the M. truncatula rhg1 homolog was translocated out of the remainder of the syntenic region on chromosome 4 or that contig 949 represents a paralogous region with the rhg1 homolog lost from the orthologous region that we examined here. Given the extensive synteny that exists in the region surrounding rhg1, it is surprising that the M. truncatula homolog of this gene, in particular, has undergone such rearrangement and/or loss. Notably, a homolog to rhg1 does exist among the network of A. thaliana duplicates (Figure 1, Synteny block 1b, At3_1iii, gene 5).

On the other hand, homologs of Rhg4 were found in the corresponding regions of M. truncatula and A. thaliana, including both M. truncatula homoeologs (Figure 2, Mt_2i gene 8, Mt_2ii gene 7, and At3_2ii gene 11). The M. truncatula Rhg4 homologs both show approximately 70% identity with Rhg4 and the surrounding region is greatly conserved as well. Mt_2i shares 75% and Mt_2ii 47%, of confirmed genes with soybean in the region immediately surrounding the Rhg4 homologs. Mt_2i even shows conservation of 100% of confirmed genes in a region over 300 kbp away. The high conservation of genome context surrounding Rhg4 indicates that, had the M. truncatula sequence been available before Rhg4 was cloned, it would have greatly facilitated cross genomic chromosome walking and cloning of the Rhg4 gene.

Tandem duplications

Genes have undergone tandem duplication in all species (soybean: 7 genes, M. truncatula: 9 genes, A. thaliana: 6 genes). In four cases, homologous soybean/M. truncatula genes are both duplicated. In no cases are homologous A. thaliana/legume genes both duplicated.

Tandemly duplicated genes with the highest copy numbers occur in a highly rearranged region in the middle of synteny block 2 (Figure 2). The rearranged region in soybean contains 11 copies of chalcone synthase genes in three separate groups of four, four, and three genes (genes 39–42, 50–51, 53–54, 57–59 in Figure 2). The latter group appears to have originated from a 25 kbp segmental duplication of the top CHS group and surrounding genes. While soybean has 11 copies of the CHS gene in this region, including CHS1, CHS2, CHS3, CHS4, and CHS5, the Mt_2ii region has only one CHS cluster with two genes, CHS1A and CHS1B (genes 52–53 in Figure 2). In addition, Mt_2ii contains a group of 10 genes with similarity to A. thaliana auxin-induced proteins 6B and X10A that are absent in soybean (genes 24–33 in Figure 2). It was not possible to analyze the corresponding region in Mt_2i, as this genome segment has not yet been sequenced.

Tandem duplications occur in other regions as well (Figures 1, 2). There are examples of tandemly duplicated genes whose homolog(s) are not duplicated, as well as cases in which two or more homologs have duplicated. For example, soybean and both M. truncatula duplicates in synteny block 2 have three copies of a glucosyltransferase (Mt_2i genes 2–4, Gm genes 18–20, and Mt_2ii genes 4–6 in Figure 2). Cases of differential tandem duplication may have resulted from duplication in only one species or loss of duplicates from one species.

Phylogenetic analysis

Phylogenetic trees were successfully generated for 21 gene families with members in synteny block 1 and for 23 in synteny block 2, many of which included homologs from both M. truncatula duplicates (Figures 1, 2) [see Additional Files 1, 2]. These phylogenies were examined to determine whether soybean and M. truncatula homologs within synteny blocks were more closely related to each other than to homologs elsewhere in the genomes, as represented by expressed sequences. For all 16 phylogenies in synteny block 1b, M. truncatula/soybean homologs were more closely related to each other than to homologs in other genomic regions, strongly suggesting orthology (Figure 1) [see Additional File 1]. Synteny block 1a contained a mix of tentative orthologs (two comparisons) and paralogs (three comparisons) (Figure 1) [see Additional File 1]. In synteny block 2, soybean genes showed orthologous relationships with their homologs in M. truncatula block Mt_2i every time (20 of 20 comparisons) and paralogous relationships in M. truncatula homoeolog Mt_2ii (11 of 11 comparisons) (Figure 2) [see Additional File 2].

Nucleotide substitution levels were determined to measure the evolutionary distance between soybean and M. truncatula (synonymous substitution levels) and to identify differences, if any, in selection pressure (nonsynonymous substitution levels). In synteny block 1, estimates of synonymous and nonsynonymous substitution levels were obtained for 34 sets of M. truncatula and soybean homologs (Figure 3a), six in synteny block 1a and 28 in synteny block 1b. In comparing these two blocks, we observed no difference in the number of synonymous substitutions per site (Table 3; 1a: 0.71, 1b: 0.71, p = 0.96), suggesting similar times of divergence between soybean and M. truncatula in both regions. This result is somewhat surprising given the fact that block 1a is composed of a mixture of apparent orthologous and paralogous relationships, while block 1b exhibits exclusively orthologous relationships.

In synteny block 2, the two M. truncatula homoeologs shared the same eight genes with soybean. We therefore focused on paired comparisons using these eight genes. The extent of synonymous substitutions between soybean and Mt_2i (0.87), an orthologous relationship based on phylogenetic tree analysis, was significantly lower than the extent between soybean and Mt_2ii (1.21), a paralogous relationship (p = 0.008) (Table 3; Figure 3b). All eight paired comparisons show higher levels of synonymous substitutions in the paralogous comparison. Not surprisingly, therefore, the paralogous region has evolved farther from soybean in evolutionary time than the orthologous region, implying that a duplication spanning the entire synteny block 2 preceded speciation between M. truncatula and soybean. The number of nonsynonymous substitution levels per site were comparable between the orthologous and paralogous M. truncatula/soybean relationships, with no significant difference (p = 0.69) (Table 3).

Comparisons of the distance between the two M. truncatula homoeologs in synteny block 2 revealed levels of synonymous substitutions (0.82) comparable to those of the orthologous Mt_2i/soybean comparison in the same block (0.87) (Table 3; Figure 3b; p = 0.85), suggesting that the duplication may have occurred close in time with speciation. It is surprising that the synonymous distance between M. truncatula homoeologs (0.82) is not closer to that of the paralogous M. truncatula/soybean comparison (1.21), which should be comparable given a duplication event followed by speciation. However, the difference between them was not significant (p = 0.22). There were no significant differences when comparing homoeologs between synteny blocks (data not shown) for either synonymous or nonsynonymous substitution levels. No differences were observed between tandemly duplicated and single copy genes for nonsynonymous or synonymous substitution levels (data not shown).

Estimates of synonymous substitution distance between orthologous soybean and M. truncatula regions allowed us to estimate the time of the Medicago/Glycine speciation event, as did the synonymous distance between the two M. truncatula duplicates in synteny block 2 in timing the underlying genome duplication. Orthologous regions between soybean and M. truncatula in both synteny blocks 1b (Mt/Gm) and 2 (Mt_2i/Gm) (Figures 1, 2) give similar estimates of divergence since speciation. Synteny block 1b has a median synonymous substitution level of 0.61 per site (Table 3), suggesting 50 My since the divergence through speciation, using an estimate of 6.1 × 10^-9 substitutions per synonymous site per year [35]. Synteny block 2 has a median synonymous substitution level of 0.59 per site (Table 3) when comparing all orthologs, suggesting 48 Mya since speciation. By contrast, the duplication event in M. truncatula evident in synteny block 2 (Figure 2) appears to have predated speciation, with a median of 0.79 synonymous substitutions per site and an inferred divergence between duplicates of 64 Mya.

Synteny

Gene density

Although soybean's genome size is more than double that of M. truncatula [51], gene density is comparable in these regions (Synteny block 1: M. truncatula = 1 gene/7.2 kb, soybean: 1 gene/6.7 kb; Synteny block 2: M. truncatula = 1 gene/6.2 kb, soybean: 1 gene/5.8 kb). These values are not far from those of A. thaliana overall (1 gene/5 kb) [52]. In this region, A. thaliana gene density is approximately half to two-thirds that of the legumes (Synteny block 1: A. thaliana = 1 gene/3.4 kb; Synteny block 2: A. thaliana = 1 gene/4.1 kb), although soybean and M. truncatula genome sizes are 7X and 3X that of A. thaliana, respectively [51]. Higher than expected gene densities in soybean and M. truncatula suggest the possibility of gene clustering. Indeed, gene clustering has identified in M. truncatula [53] and forms the basis of targeting M. truncatula sequencing to the gene-rich regions of the genome [14].

Genome duplication and speciation

The networks of synteny we have identified reflect duplication and gene loss between species, including M. truncatula regions with both orthologous and paralogous relationships to soybean. In synteny block 2, for example, there are clear-cut examples of regions with orthologous and paralogous relationships. Mt_2i shows only orthologous relationships with soybean, while Mt_2ii shows only paralogous relationships (Figures 1, 2). M. truncatula regions in synteny block 1b also unambiguously display orthology to soybean.

M. truncatula duplications in synteny block 2 allowed us to systematically examine corresponding orthologous and paralogous relationships to soybean. The percentage of conserved genes between soybean and Mt_2i (orthologous region) was twice as high as with Mt_2ii (paralogous region) (Figure 2, Table 2). Given that orthology indicates the most closely related regions evolutionarily (reflected in the phylogenetic trees) [see Additional File 2], it is not surprising that fewer genes have been deleted/inserted or experienced substitutions in orthologous comparisons. The fact that all the M. truncatula genes in the orthologous region (Mt_2i) are more closely related to soybean as evidenced by phylogenetic trees, synonymous substitution levels, percent identity, and extent of synteny, than either is to the M. truncatula genes in the paralogous region (Mt_2ii) suggests that the duplication seen in M. truncatula occurred before the speciation event splitting Medicago and Glycine lineages. Presumably, soybean also has (or had) a duplicate region as well, a possibility with some phylogenetic support [see Additional File 2].

We date the duplication event, possibly as a part of a genome duplication, in M. truncatula at 64 Mya, preceding a speciation event approximately 48–50 Mya. Indeed, the possibility of a genome duplication event predating the split between M. truncatula and soybean has been suggested previously [16, 54, 55]. Median synonymous substitution levels between the two M. truncatula duplicates in synteny block 2 (0.79 synonymous substitutions per site) fall within [55] or near [54] synonymous distance peaks, which were interpreted by the authors as a genome duplication event in M. truncatula. Schleuter et al. [55] estimates that this event occurred 58 million years ago, while Blanc and Wolfe [54] inferred a more recent event based on a substantially different molecular clock [56]. Likewise, we estimate the speciation event between Medicago and Glycine at 48 – 50 Mya, while Blanc and Wolfe [54] inferred a much more recent date of 13.3–15 million years ago, though again, the differences are primarily due to the use of differing molecular clocks [56].

Comparatively long (~500 kbp) and contiguous sets of homologous segments from different species with known phylogenetic relationships and nucleotide substitution levels bring power to the study of molecular evolution. Though median synonymous substitution levels of duplication and speciation events correspond well to published values [54, 55] (see above), the extent of synonymous substitutions varies significantly between neighboring genes despite a common genomic context (Figure 3). Estimates comparing the two M. truncatula segments created by a duplication event range from 0.62–1.12 while those comparing soybean and M. truncatula orthologs (speciation event) range from 0.42–2.68. Since the duplicates in each one of these cases presumably diverged at the same moment, one must postulate different evolutionary trajectories for the different gene lineages. Knowing that all the genes on a contiguous genomic block duplicated (and later speciated) together removes an important unknown from evolution analyses in contrast to comparable EST-based studies [54, 55].

Sequences

Glycine max sequences [GenBank:AX196294.1, GenBank:AX196295.1, GenBank:AX196297.1, and GenBank:AX197417.1] [32] were obtained from Genbank. All soybean sequences were derived from cultivar 'A3244'. [GenBank:AX196295.1] (Figure 1, gene 29) includes the susceptible allele of soybean rhg1 gene on molecular linkage group G. [GenBank:AX197297] (Figure 2, gene 21) includes the susceptible allele of soybean Rhg4 gene on molecular linkage group A2.

M. truncatula BACs were sequenced as part of an international effort to sequence the genespace of this model legume [14]. Two additional M. truncatula BACs were sequenced and examined before the international genome sequencing had begun [57]. Putative homologs of soybean sequences in M. truncatula and A. thaliana were identified by searching the soybean sequences against all sequenced M. truncatula BACs and the A. thaliana proteome using BLAST [58] (The Institute for Genomic Research, Arabidopsis Proteome version 5). After identifying genes (see below), protein/protein comparisons (BLASTP) were performed in order to confirm that BACs were syntenic and to identify syntenic genes (see below). Genbank accessions for soybean and M. truncatula sequences, A. thaliana gene numbers, and mapping information are shown in Table 1.

Sequence assemblies

Sequences were aligned and merged in regions of sequence overlap on the basis of 99% identity or better. End-sequenced BAC clones that tentatively spanned gaps in the sequence were identified based on strong hits (e-value = 0, ≥99% identity) to sequenced BACs on either side. Gap sizes were estimated by removing overlap from the estimated size of end-sequenced BAC(s).

Nomenclature

Throughout the manuscript, the following nomenclature is used. Regions surrounding and syntenic to the SCN resistance rhg1 locus are collectively referred to as synteny block 1 (Figure 1). Regions surrounding and syntenic to SCN Rhg4 gene are collectively referred to as synteny block 2 (Figure 2). Synteny block 1 is divided into blocks 1a and 1b, which are separated by gaps in all three species (Figure 1). Within each synteny block, species are labeled as Gm (soybean), Mt (M. truncatula), or At (A. thaliana). The chromosome number follows the "At" abbreviation for A. thaliana. If more than one homoeolog is present, the species abbreviation is appended with an underscore followed by the synteny block and lower case roman numerals (i.e. Mt_2i, Mt_2ii, At4_2i, and At3_2ii in Figure 2) (Figures 1, 2). Sequence assemblies separated by physical gaps are labeled as sequence assemblies "top" and "bottom" in arbitrary order (Figure 1).

Gene prediction and identification of synteny

Genes were predicted in G. max and M. truncatula genomic sequences using the dicot (Arabidopsis) matrix of FGENESH [59, 60]http://www.softberry.com. BLASTP was used to compare predicted proteins between databases containing these G. max or M. truncatula predicted genes and all A. thaliana proteins with an e-value cutoff of e-8 and percent identity cutoff of 40% for the top high scoring segment pair for soybean and M. truncatula comparisons and an e-value cutoff of e-8 for comparisons to A. thaliana [58]. These cutoff values generally identified homologs in syntenic positions while rejecting related genes in nonsyntenic positions.

In this study, we defined synteny to include both conservation of gene content and order between species. In estimating syntenic density (the percentage of genes conserved between two species), repetitive sequences (genes with similarity to transposable elements, including retroelements) were not included and tandemly duplicated genes were counted as one. Synteny between two species was estimated from the first to the last pair of conserved genes in the available sequence for both species.

Phylogenetic analysis

To distinguish between orthologous and paralogous regions, we constructed phylogenetic trees as follows. BLASTP or TBLASTN, as appropriate, were used to compare all G. max genes with the following sequences: all G. max and M. truncatula proteins in the corresponding genomic regions of this analysis; the nonredundant A. thaliana proteome; soybean and M. truncatula EST unigene sets [61] (GMGI v.11 and MTGI v.7; The Institute for Genomic Research. Rockville, MD). The top 25 hits ≥100 amino acids with e-values ≤ e-10 were included in the analysis. Tandem duplications and highly related genes in the same gene family were grouped for analysis.

Initial alignments were calculated using T-COFFEE [62] with manual evaluations and edits in Jalview [63] for poorly aligning sequences. For subsequent phylogenetic analysis, an HMM calculated for each alignment using hmmer [64] was used to realign sequences and to identify and remove indel regions and sequences with fewer than 60% matches to the model. Parameters for hmmbuild were: archpri = 0.7, gapmax = 0.3.

Parsimony trees were calculated using the protpars of Phylip [65], with maximum likelihood branch lengths calculated using TREE-PUZZLE [66]. Parameters for protpars were: randomize input order; use ordinary parsimony; search for best tree; select one best tree for further analysis in TreePuzzle. Parameters for TreePuzzle were: user defined tree (from parsimony search); approximate parameter estimates; Whelan-Goldman substitution model [67] estimate amino acid frequencies from data set; allow rate heterogeneity with eight gamma-distributed rates.

Nucleotide substitutions

Codon-aligned nucleic acid sequences were created with TranslateAlign.pl (courtesy Dan Kortschak, University of Adelaide, Adelaide, Australia). Nucleotide substitutions levels were calculated using these alignments with SNAP (Synonymous/Non-synonymous Analysis Program) [68, 69]. In this program, the levels of synonymous and nonsynonymous substitutions per site are approximated using methods developed by Nei and Gojobori [70], incorporating Ota and Nei's statistic [71]. Median synonymous substitution levels were converted into estimates of time since divergence using an estimate of 6.1 × 10^-9 substitutions per synonymous site per year [35].