Adaptive evolution in the Arabidopsis MADS-box gene family inferred from its complete resolved phylogeny

Results

Annotation, Nomenclature, Gene Structure, and Duplications. The list of 104 MADS-box sequences found in the Arabidopsis genome database (www.ncbi.nlm.nih.gov/blast/Genome/ara.html), given agamous-like (AGL) number, accession number, chromosome location, intron-exon structure, and type of duplication are shown in Table 3.

Type I and II genes have contrasting gene structure and chromosome location. Type I genes have always one or two exons, whereas type II genes have more than five, and typically six to eight, suggesting that these two types of genes have different predisposition to gain or loss of introns or, alternatively, a difference in exon shuffling in the building of both types of genes, perhaps since their origin. All type II genes have a clear MIKC structure, except AGL33 with a very short transcript and AGL30-related sequences that do not have a clear IKC structure but share conserved motifs in their MADS-boxes with the rest of the type II genes. Also, a coiled-coil domain similar to the K domain is inferred for at least AGL104 (data not shown).

We used chromosome map locations to make qualitative inferences on past duplication events during MADS-box gene family evolution. The distribution of type I sequences among the five chromosomes is distinct to that of type II genes. Whereas the former are concentrated in chromosome I and V (χ² = 26.77; P < 0.001 rejects uniform distribution among chromosomes correcting with chromosome size), type II are uniformly distributed (χ² = 2.62; P > 0.1) among chromosomes. Also, most type I genes can be traced to intrachromosomal duplications, whereas approximately half of type II genes seem to have originated from interchromosomal duplications. Interestingly, Lynch and Conery (39) have found that recent duplications happened more frequently within than between two chromosomes, suggesting that most type I genes diverged more recently than type II genes (40). Moreover, a survey of the A. thaliana paralogous blocks database (http://wolfe.gen.tcd.ie/athal/dup) indicates that there are more duplicates from the type II group that seem to have persisted than those from type I. Eighteen out of 26 type II sequences that are found in a nonspurious duplicated region had a close paralog in a sister region, whereas among type I genes only 2 out of 22 did. Among type II genes found in nonspurious duplicated blocks, 14 are found in interchromosomal duplications and 12 in intrachromosomal duplications, whereas among type I genes the corresponding numbers are 15 and 7. But when we consider gene pairs of terminal clades of the trees, 13 out of 18 type I gene pairs involve intrachromosomal duplications, whereas in the type II we found this to be true of 6 out of 15 gene pairs (Fig. 1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

A. thaliana MADS-box gene family Bayesian phylogeny. Numbers above or below branches represent Bayesian posterior probabilities of finding a given clade. Branch lengths are proportional to number of nucleotide substitutions. AGL105 was excluded.

MADS-Box Gene Family Phylogeny. To corroborate the monophyletic origin of the two lineages that we had previously proposed (14), we obtained trees that included 103 of the MADS-box sequences found in the A. thaliana genome (Fig. 1). The global phylogeny recovers the two lineages (types I and II) of MADS-box sequences as two monophyletic groups with both alignments used (see supporting information) if the tree root is placed between type I and II genes. However, type I genes are more divergent among them than type II genes (Tables 4 and 5, which are published as supporting information on the PNAS web site). Nonetheless conserved motifs after the MADS suggest that type I sequences are not pseudogenes, an idea supported by the recent characterization of PHERES1 (16).

AGL30 had been incorrectly assigned to type I. However, this and related genes seem to be divergent type II. This is supported by their affiliation to type II-like moss genes that bear K domains as well as by their exon number (27, 41), and by conserved MADS-box motifs with respect to other type II genes. In our global tree (Fig. 1), these genes are resolved in a different position to that in the type II tree (Fig. 2b). Nonetheless, this remaining ambiguity does not affect PDS analyses. We ran PDS analyses for MIKC genes with an alternative type II topology similar to that in Fig. 1 and results recover the same sites with high PP as those obtained with the topology shown in Fig. 2b (data available on request).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Type I and type II A. thaliana MADS-box Bayesian phylogenies. (a) Type I tree polarized with type II sequences. (b) Type II tree polarized with type I sequences. Numbers above or below branches represent posterior probabilities. Branch lengths are proportional to the number of nucleotide substitutions. Boxed letters identify clades in which site-specific tests of positive selection yielded statistically significant LRT results for at least two model comparisons and in which at least one of the models detected sites under PDS with PP > 0.90. Branches underlined with a broken line identify cases in which the branch-site analyses yielded significant PDS results. #, Excluded from the site-specific analyses.

Given the ambiguity of the alignment of sequences after the MADS-boxes of both lineages, we resolved the internal phylogenetic relationships of each lineage separately. For the type I lineage tree, we used type II sequences as outgroups and vice versa for type II lineage. In both cases, tree topologies were very similar if regions downstream of the MADS-boxes of outgroup sequences were assumed to be homologous to those of ingroup sequences or were displaced. In contrast to previously published phylogenies, the trees shown here for both lineages resolve the branching sequence of the monophyletic groups (Fig. 2).

In the type I lineage, several strongly supported monophyletic clades are resolved and these are confirmed in trees with various outgroups (data not shown). It is noteworthy that no type I sequences bear the IKC region typical of lineage II and that several of the previously identified (14) privative amino acids within the putative MADS-domain of this lineage are found in most available sequences. Two main monophyletic groups are distinguished within the type I, suggesting at least one ancestral duplication within this lineage. DNA sequences beyond the MADS-box in the AGL23-like sublineage are conserved within each small monophyletic group resolved within this sublineage but are very divergent among groups (Fig. 2a and Table 4). Within the AGL26-like genes, two groups are resolved [AGL26-like genes themselves and PHERES1-like genes that includes the first type I gene functionally characterized (16)] and within each a very high degree of conservation is found in the putative domains beyond the MADS (Fig. 2a and Table 4). Sisters to these groups, AGL47 and AGL82, are resolved but have divergent putative domains beyond the MADS.

In the type II lineage, which includes all of the MIKC genes functionally characterized up to now, several clades are resolved and all are well supported (Fig. 2b). It is noteworthy that the AGAMOUS clade has AGL12, which had been reported as root-specific (42), as its sister gene. Another important finding is the strong association of AGL79, expressed in roots (data not shown), with AP1, CAL, and FUL, which are well characterized flower development genes (43). Therefore, this tree suggests that not all monophyletic groups resolved include genes with similar expression patterns and functions as previously thought (15, 21), but formal and robust inferences on evolution of MADS-box gene expression and function will have to await more experimental data and the inclusion of genes from additional taxa.

Positive Selection in MADS-Box Gene Evolution. We compared Models M0 and M3 to evaluate whether there had been dN/dS ratio variation among codon positions below each node of type I and II trees from Fig. 2 (Tables 1 and 2). We found rate variation at deep, intermediate, and recent duplications of both type I and II lineages (data available on request).

Secondly, we applied the LRT to compare data fit to models M1 vs. M2 and M7 vs. M8 to address whether PDS promoted divergence of MADS-box genes below nodes and whether the action of selection has been heterogeneous among protein domains codified by these genes, using the trees of Fig. 2 (Tables 1 and 2 and Fig. 3). Below many of the deep nodes of type II tree, model M3 and at least one of M2 or M8 had significant LRT results (Fig. 2b). However, none of the sites with high PP were detected by more than one model comparison with significant LRT results (data available on request). In contrast, a similar analysis for type I lineage reveals several nodes below which specific sites appear to have been under PDS (Tables 1 and 2). For instance, in nodes AH and AL, positions 72-74 appear to have been under PDS with high PP, with position 72 showing five different amino acids for the seven sequences involved (Fig. 2a and Table 4). Positions 73 and 74 are part of the otherwise highly conserved “RQVTF” motif, and in the human serum response factor, position 72 has been shown to be involved in DNA contact (44). Below nodes AH and AR, position 82 was also found to be under PDS in two of the model pairs compared, although model M3 collapsed to only two rate classes. At this position, amino acid diversity is very high. For example, 10 different residues can be found for 20 sequences analyzed below node AR. In contrast, the homologous position of type II genes (position 26 in Table 4) shows only six different amino acids for 45 sequences.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Schematic representation of the distribution of sites under PDS in type I (a) and type II (b) sequences. MADS, I, K, and COOH domains are indicated. In a, the upper row corresponds to branch-site analysis along the branch leading to PHERES1, and the rest of the rows correspond, from top to bottom, to the site analyses performed below the nodes marked with the boxed letters AG, AH, AL, AP, AQ, and AR in Fig. 2a. In b, the upper row corresponds to branch-site analysis along the FLC branch and the lower row corresponds to branch-site analysis along the SVP branch. In a, the upper scale corresponds to amino acid position along the MADS domain and the lower scale corresponds to amino acid position along our alignment. In b, the scale corresponds to amino acid position along our alignment. Sites with PP > 0.70 are included. Shaded regions in b were excluded from PDS to avoid false positives. All sites are listed in Tables 1 and 2.

Strong evidence for positive selection in type I evolution was also found in less variable positions. This is the case for position 123 of node AG and position 58 of nodes AP and AQ. In position 123, there are only two different amino acids, although their distribution suggests that this site mutated twice during the history of descendants of node AG. In position 58, there are only five variable sites out of 14 sequences compared.

In the above analyses, PDS is detected at individual sites only if the average dN rate across lineages is higher than the average dS rate. This is observed mainly when recurrent positive selection occurs. However, PDS may change a few key residues of a protein but only at particular moments during its evolutionary history. In the latter case, detecting a significantly elevated dN would be hard if an average across-lineages estimate is considered. Thus, we applied the branch-site model of Yang and Nielsen (33) to test for PDS affecting individual sites along the branches leading to the flowering-time FLC- and SVP-like genes in type II lineage and along the branch leading to the only functionally characterized type I protein (PHERES1).

Parameter estimates suggest that in the basal branch of the FLC-like lineage, four residues were fixed by PDS with PP > 0.95 (Tables 1 and 2 and Fig. 3). Two of the residues were found within the MADS and correspond to amino acids that in other MADS-domain proteins participate in interactions between subunits (44-46). For example, site 42 is homologous to a site that in the myocyte enhancer factor-2 has been shown to intervene in subunit folding (45). Another site with high PP is 154. This site is within the K-box and has been shown to intervene in AP3/PI dimerization and determines functional specificity in AP1 and AG (47).

PDS seems to have been important also along the branch leading to the SVP-like genes, and we found two sites that appear to have been fixed by PDS (PP > 0.95). Site 16 is within the MADS-domain, and its homologous position in MEF2A plays a role in DNA-protein complex stabilization (46). Position 144 is found within the K-domain. SVP and AGL24 are the only type II proteins that have a lysine at that position.

Branch-site models also detected strong (PP > 0.95; see Tables 1 and 2) support for PDS along the branch leading to PHE1 at sites 92 and 72, which have been reported (44) to be important for α-β folding, and for position 105, which is involved in both dimerization and α-β folding.

Our analyses suggest that gene conversion or concerted evolution has not been prevalent during MADS-box genes evolution and does not bias PDS inferences. The overall modified Nei-Gojobori (48) means of synonymous-nonsynonymous differences were relatively high: 0.18 and 0.09 (transition/transversion ratio of 1.72 and 1.7, respectively) for types I and II, respectively (49). Gene conversion does not seem to have played significant roles in MADS-box evolution because only three conversion events were detected among type II genes, but none of these included genes for which we infer PDS and no conversion events were detected among type I sequences.

Discussion

We presented an annotated list of 104 MADS-box sequences from the complete A. thaliana genome database. Our phylogenetic analyses provide a resolved evolutionary hypothesis for the A. thaliana MADS-box gene family. This will be a useful reference for establishing orthology relationships, postulating functional hypotheses for uncharacterized MADS-box genes, and evaluating the role of MADS-box genes in plant morphological evolution.

The monophyly of the type II lineage is strongly supported in the present analyses, and type I comprises several sublineages with divergent putative domains after the MADS. However, previous analyses (14), as well as contrasting exon-intron structure and chromosomal distribution between type I and II sequences, still support the placement of the root between the type I and II genes in the tree of the complete gene family (Fig. 1). This tree hence resolves type I and II sequences in two monophyletic lineages. Nonetheless, genes from other plant, animal, and fungal species should be included in future analyses to trace MADS-box gene duplications with respect to taxa divergence and to specifically reevaluate the number of MADS-box gene duplications that occurred before the divergence of plants and animals (14). Such analyses will provide further evidence to reevaluate the monophyly of type I and II lineages.

Gene family structure has to be understood in the context of extensive gene duplications that have occurred in the evolutionary history of A. thaliana (40). Duplications leading to the chromosome stretches identified in the Arabidopsis Genome Initiative occurred 65 million years ago or before (refs. 50 and 51, but see ref. 52). Most retained groups within these stretches belong to type II genes, and duplications among type II seem to have been more ancient than those among type I, as suggested by their differential distribution among chromosomes. Retention due to a balance between genetic drift and mutation (5, 53) would depend on population characteristics (mainly effective population size) and hence would affect sequences of type I and II equally. But contrasting roles of selection between these two lineages could underlie the contrasting retention rates observed between them. The different evolutionary histories could have been determined by the fact that genes from these two lineages were recruited for different functions.

Interestingly, although duplications of type I seem to have occurred more recently than those leading to the type II lineage, type I sequences are more divergent among them in comparison to type II. This finding would suggest that whereas type II genes have been affected by sporadic PDS at the origin of new functions, followed by strong functional constraint, type I genes have been subject to recurrent events of PDS. In turn, this could be an indication that the functional roles of type I genes are overall distinct to those of type II genes. This pattern also allows us to put forward the hypothesis that type I orthologues from other taxa are less conserved than most type II orthologues.

Indeed, PDS analyses presented here suggest that type I and II lineages have been subject to overall contrasting selection regimes. We found that recurrent positive selection could have played a role in fixing specific amino acids after several duplication events during the evolution of type I genes. In contrast, analyses of site models did not provide strong evidence for PDS selection among type II genes. We had to use site-branch models to detect a role for PDS among type II genes. Indeed, we found evidence for PDS along the branches leading to the groups of genes that control flowering time that evolved a new function with respect to most other genes characterized up to now that are involved in cell- or organ-type specification. Indeed, probably by their control of life-history traits, flowering-time genes may have directly impacted plant fitness, and this could also explain the prevalence of positive selection during protein evolution among them.

Sites with high PP of having been fixed by natural selection in both lineages were found mainly in the MADS and K domains. However, our analyses are biased toward these domains because we only focused on conserved stretches that may be unambiguously aligned and excluded most variable domains. Future studies focusing on particular closely related genes for several species will be useful to address the role of PDS within COOH and other divergent domains. Indeed, a recent study showed that regions within the C-terminal domain determine functional specificity in AP3 and PI and may be relevant for floral organ evolution (54). The localization of the sites with high PP identified here suggest a role for PDS in MADS-domain protein diversification through interactions with protein partners and changes in affinity to binding motifs (46, 47, 55). Sites and genes identified here to have been under PDS become interesting targets for functional evaluations.

Evaluations of assumptions and predictions made by models of gene duplication and persistence will require phylogenetically driven analyses of functional and population level data for related genes in different species. The MADS-box gene family might become a good “model family” for such a purpose. For example, our site model analyses did not find that PDS played a role in the divergence of redundant AP1, CAL, and FUL (56). However, population-level data do suggest a role for positive selection in the divergence of these genes (57). More powerful analyses should be performed to rule out false negatives in our analyses due to low gene number (32). Moreover, our conclusions were based on conservative analyses and unambiguously aligned sequences to avoid false positives. Additional tests (data available on request; Fig. 2b) suggest that the role of PDS in MADS-box gene evolution might be more widespread. Other approaches (58) and the inclusion of sequences for additional taxa should be considered when further investigating the role of PDS in MADS evolution.

Our results suggest a role for positive selection during MADS-box evolution in plants. Previous studies have emphasized the role of changes in cis-regulatory regions of transcriptional regulators during plant evolution (59). Fewer recent studies, however, have also demonstrated that the evolution of transcriptional regulators' cDNA sequences played important roles in plant evolution (54). The detection of positive selection in MADS protein sequences that are developmentally important also indicates that changes in cDNA, and not only in the regulatory regions of these genes, have played a role in the evolution of plant body plans.