Main

Optical mapping indicates that the eight pseudomolecules of assembly Mt3.5 span a physical distance of 375 million base pairs (Mb), and fluorescence in situ hybridization indicates they extend from pericentromeres almost to telomeric ends (Supplementary Figs 1 and 2). Altogether, Mt3.5 consists of 2,536 bacterial artificial chromosomes (BACs; Supplementary Tables 1 and 2) with 273 physical gaps (including centromeres, Supplementary Table 3) and 101 internal sequencing gaps. The pseudomolecules contain 246 Mb of non-redundant sequence (Supplementary Table 2) located entirely within the optical map (Supplementary Fig. 3). Another 146 unfinished BACs/BAC pools that cannot be placed on the optical map contribute 17.3 Mb. Regions not represented in pseudomolecules or unanchored BACs were captured through assembly of approximately 40× coverage Illumina sequencing, yielding 104.2 Mb of additional unique sequence. Although not directly tested, the Illumina sequence is expected to lie predominantly within the boundaries of pseudomolecules (see below). On the basis of expressed sequence tag alignments, the combined data sets capture 94% of expressed genes, providing a highly informative platform for analysing the euchromatin of M. truncatula, although still at the draft stage.

Altogether there are 62,388 gene loci in Mt3.5 (Supplementary Table 4 and Supplementary Fig. 4), with 14,322 gene predictions annotated as transposons. Pseudomolecules and unassigned BACs contain a total of 44,124 gene loci, 177,271 retroelement-related regions and 26,487 DNA transposons, and non-redundant Illumina assemblies contribute an additional 18,264 genes, 75,777 retrotransposon regions and 8,476 DNA transposons (Supplementary Tables 5–9) along with 1,418 organellar insertions (Supplementary Data 1). For pseudomolecules and unassigned BACs, this translates to 16.8 genes, 67.6 retrotransposons and 10.1 DNA transposons per 100 kilobases (kb). Within Illumina sequence assemblies, gene density (17.1 per 100 kb) and retrotransposon density (72.2 per 100 kb) are similar to pseudomolecules and unassigned BACs, whereas DNA transposon density is lower (8.2 per 100 kb). Similarities in gene and transposon densities between BAC and Illumina sequences support the assertion that the Illumina sequence is euchromatic, although the possibility that some Illumina assemblies come from low-copy regions within heterochromatin can not be excluded. Considering only the 47,845 genes with experimental or database support (Supplementary Table 4), the average M. truncatula gene is 2,211 bp in length, contains 4.0 exons, and has a coding sequence of 1,001 bp. These values are similar to those observed previously in Arabidopsis thaliana (2,174 bp), Oryza sativa (3,403 bp) and Populus trichocarpa (2,301 bp)4,5,6.

Recent analyses of plant genomes indicate a shared whole-genome hexaploidy (WGH) preceding the rosid–asterid split at 140–150 Myr ago7. Duplication patterns and genomic comparisons strongly suggest an additional WGD approximately 58 Myr ago in the papilionoids8,9. Near the time of this WGD, papilionoids radiated into several clades, the largest of which split quickly into two subclades, the Hologalegina (including M. truncatula and L. japonicus) and the milletioids (including G. max and other phaseoloids) at about 54 Myr ago2. We therefore compared M. truncatula pseudomolecules with other sequenced plant genomes to learn more about shared synteny and genome duplication history.

There is significant macrosynteny among M. truncatula, L. japonicus and G. max (Fig. 1 and Supplementary Fig. 5a, b). Conserved blocks, sometimes as large as chromosome arms, span most euchromatin in all three genomes. A given M. truncatula region is typically syntenic with one other M. truncatula region as a result of the approximately 58-Myr-ago WGD, usually in small blocks showing degraded synteny (Fig. 2 and Supplementary Fig. 6). A given M. truncatula region is most similar to two G. max regions via speciation at about 54 Myr ago and the Glycine WGD at <13 Myr ago10 and less similar to two other G. max regions resulting from the 58-Myr-ago and <13-Myr-ago WGD events. A M. truncatula region is likewise most similar to one L. japonicus region via speciation at about 50 Myr ago and less similar to a second L. japonicus region as a result of the 58-Myr-ago WGD. Finally, each M. truncatula region and its homeologue typically show similarity to three Vitis vinifera regions via the pre-rosid WGH. Exceptions to these patterns could be due to gene losses, gains, or rearrangements specific to the M. truncatula lineage, resulting in synteny being more evident between M. truncatula and other genomes than in self-comparisons. Indeed, self-comparisons within M. truncatula reveal few remnants of the legume-specific WGD (Fig. 2 and Supplementary Fig. 6). Whereas this seems paradoxical, it is probably explained by extensive gene fractionation between WGD-derived homeologues in M. truncatula. In Fig. 3, two short regions on Mt1 and Mt3 resulting from the 58-Myr-ago WGD are displayed beside microsyntenic regions of G. max and V. vinifera. As expected, many genes are microsyntenic between M. truncatula and G. max (ranging from 7/19 between Mt3 and Gm14 to 10/20 between Mt1 and Gm17). Between the two M. truncatula homeologues, however, only 6 out of 33 genes (or collapsed gene families) are microsyntenic, with a homeologue missing from one or the other duplicate (Supplementary Table 10). Apparently, there have been many more changes, large and small, in M. truncatula than in G. max since the legume WGD. This is borne out by the fact that synteny blocks in M. truncatula are one-third the length of those remaining from the papilionoid WGD in G. max (524 kb against 1,503 kb) with the average number of homologous gene pairs per block correspondingly lower (12.4 against 31.0).

Figure 1: Circos diagram illustrating syntenic relationships between Medicago, Glycine, Lotus and Vitis.
figure 1

Homologous gene pairs were identified for all pairwise comparisons between M. truncatula, G. max, L. japonicus and V. vinifera genomes. Syntenic regions associated with the ancestral WGD events were identified by visually inspection of corresponding dot-plots. The large Mt5–Mt8 synteny block (yellow) was found to have two syntenic regions in L. japonicus (red), four syntenic regions in G. max (blue) and three in V. vinifera (green).

PowerPoint slide

Figure 2: Circos diagram illustrating the Medicago WGD and selected gene families.
figure 2

The 963 WGD-derived paralogous gene pairs were examined for overlap with the nodule-enhanced gene list (Supplementary Data 2). Resulting gene pairs were joined and plotted as either blue triangles (only one of the duplicates is nodule-enhanced) or red (both nodule enhanced). Gene densities of NBS-LRRs, NCRs and other defensin-like proteins are plotted against chromosome position. Density was calculated using a sliding window (100-kb window with 50-kb steps).

PowerPoint slide

Figure 3: Microsynteny comparison between Medicago homeologues and corresponding regions of Glycine and Vitis.
figure 3

Microsyntenic genome segments are centred around Medtr3g104510/Medtr1g015890 (Supplementary Table 10), a duplicated region derived from the 58-Myr-ago WGD event noted in orange. The <13-Myr-ago G. max-specific WGD is coloured yellow. Orthologous/paralogous gene pairs are indicated through use of a common colour. White arrows represent genes with no syntenic homologue(s) in this genome region. Some of these genes may actually have a syntenic sequence in soybean but no corresponding model reported in the current annotation (http://www.phytozome.net/soybean).

PowerPoint slide

The M. truncatula genome also has undergone high rates of local gene duplication. The ratio of related genes within local clusters compared to all genes in families is 0.339 in M. truncatula, 3.1-fold higher than in G. max and 1.6-fold higher than in A. thaliana or P. trichocarpa. (‘Local clusters’ are defined as genes in a family all within 100 gene models of one another.) The excess of local gene duplications in M. truncatula is observed genome-wide and affects many families. There are 2.63 times as many gene families with local duplications in M. truncatula compared with G. max (2,980 against 1,131), an excess that also is seen in detailed comparisons of syntenic regions in M. truncatula and G. max. We examined 16.3 Mb of Mt05 showing synteny to two large regions of Gm01 plus homeologous blocks on Gm02, Gm09 and Gm11. In these regions, 25.8% of M. truncatula genes are locally duplicated compared with just 8.0% in G. max. Local gene duplications and losses have contributed both to synteny disruptions (Fig. 3 and Supplementary Fig. 7) and to high gene count (62,388) in M. truncatula—a value nearly as high as the 65,781 total gene models in G. max despite its additional (<13 Myr ago) WGD. Local gene duplications are evident in certain gene families, such as F-box genes, which have undergone pronounced expansions (Supplementary Fig. 8 and Supplementary Table 11). M. truncatula also has experienced higher rates of base substitution compared to other plant genomes (Supplementary Fig. 9). Assuming 58 Myr ago as the date of the legume WGD, then the rate of synonymous substitutions per site per year in M. truncatula is 1.08 × 10−8, 1.8 times faster than estimates in other vascular plants11. Higher rates of mutation and greater levels of rearrangement in M. truncatula following the papilionoid duplication may have been driven by factors including short generation times, high selfing rates or small effective population sizes, although these characteristics are not unique to M. truncatula.

Legumes and actinorhizal species are capable of forming a specialized organ, the root nodule, a highly differentiated structure hosting nitrogen-fixing symbionts. Phylogenetic studies suggest that nodulation may have evolved multiple times in the Fabidae, but the observation that all nodulating species are contained within this single clade indicates that a predisposition to nodulate evolved in their common ancestor12. It is unknown whether nodulation with rhizobia preceded the divergence of the three legume subfamilies or evolved on multiple occassions13. Nevertheless, rhizobial nodulation and the 58-Myr-ago WGD are features common to most papilionoid legumes and both occurred early in the emergence of the group2. Given that WGDs generate genetic redundancy that potentially facilitates the emergence of novel gene functions without compromising existing ones14, we examined the M. truncatula genome to ask whether the 58-Myr-ago WGD might have had a role in the evolution of rhizobial nodulation in M. truncatula and its relatives.

Nod factors are bacterial signalling molecules that initiate nodulation. Previous studies have shown that several of the plant components involved in the response to Nod factors also function in mycorrhizal signalling15. However, some Nod factor receptors and transcription factors have distinctly nodulation-specific functions. Among these nodulation-specific components, we found that the Nod factor receptor, NFP, and the transcription factor, ERN1, each have paralogues, LYR1 and ERN2 respectively, that trace back to the papilionoid WGD based on genome location and synonymous substitution rate values (Supplementary Fig. 10 and Supplementary Data 2). Both sets of gene pairs also show contrasting expression patterns and functional specialization. NFP and ERN1 are expressed predominantly in the nodule and are known to function in nodulation16,17, whereas LYR1 and ERN2 are highly expressed during mycorrhizal colonization (Supplementary Fig. 11). These observations indicate that two important nodulation-specific signalling components in M. truncatula might have evolved from more ancient genes originally functioning in mycorrhizal signalling and then duplicated by the 58-Myr-ago WGD. In the case of M. truncatula NFP/LYR1, this conclusion is supported by the observation that the apparent orthologue of NFP in the nodulating non-legume Parasponia andersonii functions in both nodule and mycorrhizal signalling18. Thus, the 58-Myr-ago WGD seems to have led to sub-functionalization of an ancestral gene participating in both interactions, resulting in two homeologous genes that each performs just one of the original functions.

To assess further the contribution of the WGD to M. truncatula nodulation, we analysed expression of paralogous gene pairs using RNA-seq data from six different organs (Supplementary Methods 5.1). A total of 963 WGD-derived gene pairs were found (Supplementary Data 2) with 618 pairs (1,046 genes) having RNA-seq data for one or both homeologue. We then determined the number of genes showing organ-enhanced expression (defined as genes with expression level in a single organ at least twice the level in any other) within the pseudomolecule and the WGD-derived gene sets (Supplementary Table 12). In both cases, different organs contained markedly different numbers of genes with enhanced expression (χ2 with 5 degrees of freedom, P = 10−272); however, the rank order among the organs was identical. Roots had the largest number of genes with enhanced expression followed by flower, nodule, leaf, seed/pod and bud. Among gene pairs with nodule-enhanced expression, both paralogues were nodule-enhanced in eight pairs, whereas just a single paralogue was nodule-enhanced in the other 43 pairs. This is consistent with nodulation pre-dating the WGD and further sub- and neo-functionalization emerging afterwards. We went on to examine transcription factors because they can act as regulators of plant growth and development. A total of 3,692 putative TF genes were discovered (Supplementary Data 3), representing 5.9% of all M. truncatula gene models (Supplementary Table 13). Of the 1,513 TF genes on pseudomolecules with RNA-seq data, 142 genes (9.4%) derived from the 58-Myr-ago WGD (Supplementary Fig. 12 and Supplementary Data 4), consistent with previous observations indicating greater retention of transcription factors following polyploidy19. Nodule-enhanced expression was significantly higher among transcription factors (92 out of 1,513 or 6.1%) than among all pseudomolecule genes (1,111 out of 23,478 or 4.7%) (χ2 with 1 degree of freedom, P = 0.024) (Supplementary Table 12). Nodule-enhanced expression was even higher in WGD-derived transcription factors (11 out of 142 or 7.7%), although this enrichment did not reach statistical significance (P = 0.113). As expected, ERN1 is found within this group of WGD-retained, nodule-enhanced transcription factors.

These results show that many paralogous genes retained from the 58-Myr-ago WGD, especially signalling components and regulators, have undergone sub- or neo-functionalization, including several with specialized roles in nodulation. Nevertheless, separate phylogenetic analyses (Supplementary Methods 5.5) indicate that some nodule-related genes derive from the more ancient pre-rosid WGH, with their nodule-related functions pre-dating the 58-Myr-ago WGD (Supplementary Data 5). Taken together, these results are consistent with a model where the capacity for primitive interaction with new symbionts derived from existing mycorrhizal machinery involving genes recruited from the pre-rosid WGH. This capacity would have arisen early in the Fabidae clade and led to the appearance of nodulation in multiple lineages13,20. Later, the 58-Myr-ago WGD would have resulted in additional genes, including NFP, ERN1 and the transcription factors described above, that went on to become specialized for nodule-related functions in the Papilionoideae.

Medicago contains additional amplified gene families, many nodulation-related and found in tandem clusters. M. truncatula has nine symbiotic leghaemoglobins, more than twice the number in L. japonicus or G. max (Supplementary Fig. 13). Five of these genes are located in a tight cluster on Mt5. The M. truncatula genome contains 593 nodule cysteine-rich peptides (NCRs) (Supplementary Data 6), a gene family restricted to M. truncatula and its relatives21. NCRs are noteworthy because they include members essential for terminal differentiation of rhizobia22. NCRs are tightly clustered within the M. truncatula genome (Fig. 2), with 75% found in clusters of up to 11 members. The M. truncatula genome also has 764 nucleotide-binding site and leucine-rich repeat (NBS-LRR) genes (Supplementary Data 7), more than other plant genomes that have been sequenced so far23,24,25, many with nodule-specific expression (Supplementary Fig. 14). Almost 90% of NBS-LRRs occur in clusters and genome regions showing limited macrosynteny to other species, such as Mt3 and Mt6, are locations of large NBS-LRR superclusters (Fig. 2 and Supplementary Tables 14 and 15). Finally, M. truncatula secretes flavonoid signalling molecules to induce the nod genes of Sinorhizobium meliloti26. In M. truncatula, the corresponding biosynthetic pathway has expanded markedly, with 28 M. truncatula chalcone synthase genes in clusters of up to seven members compared to just four chalcone synthases in A. thaliana27 (Supplementary Data 8). M. truncatula has ten chalcone reductases compared to none in A. thaliana28 and M. truncatula has 11 chalcone isomerase genes, including one cluster of seven members, compared to just one representative in A. thaliana29 (Supplementary Figs 15 and 16).

Analysis of the M. truncatula genome supports earlier studies indicating that the dramatic radiation of the legume family (at least the papilionoid subfamily) is partly attributed to the 58-Myr-ago WGD30. Our results indicate that the WGD early in papilionoid evolution allowed the emergence of critical components in Nod factor signalling and contributed to the complexity of rhizobial nodulation observed in this clade. As such, the WGD seems to have had a crucial role in the success of papilionoid legumes, enhancing their utility to humans.

Methods Summary

DNA sequencing

Six A17 BAC and one fosmid library were used to create Mt3.5 (Supplementary Table 1). Most were processed by Sanger paired-end sequencing of 3–6-kb shotgun libraries. Sequences were downloaded in February/March 2009 with scaffolding performed by aligning all BAC and fosmid ends against contigs and then anchored and ordered primarily by optical mapping. Separately, 25 billion base pairs (Gb) of Illumina sequence was generated using short (375 nt) inserts plus 2.1 Gb from a 5 kb mate-pair library, then assembled using CLCbio (http://www.clcbio.com) and Soap (http://soap.genomics.org.cn/).

RNA sequencing

Five tissues were used for RNA-seq analysis with 10 million Illumina 36-bp reads per library (Supplementary Table 12). Three tissues were used for small RNA analysis with 3 million reads per Illumina library (Supplementary Figs 17–18, Supplementary Table 16 and Supplementary Data 9).