The Evolutionary Implications of knox-I Gene Duplications in Conifers: Correlated Evidence from Phylogeny, Gene Mapping, and Analysis of Fun
http://www.100md.com
分子生物学进展 2004年第12期
* Arborea and Chaire de Recherche du Canada en Génomique Forestière et Environnementale, Centre de Recherche en Biologie Forestière, Université Laval, Sainte-Foy, Québec Canada; and Natural Resources Canada, Canadian Forest Service, Laurentian Forestry Centre, Sainte-Foy, Québec Canada
E-mail: bousquet@rsvs.ulaval.ca.
Abstract
Class I knox genes code for transcription factors that play an essential role in plant growth and development as central regulators of meristem cell identity. Based on the analysis of new cDNA sequences from various tissues and genomic DNA sequences, we identified a highly diversified group of class I knox genes in conifers. Phylogenetic analyses of complete amino acid sequences from various seed plants indicated that all conifer sequences formed a monophyletic group. Within conifers, four subgroups here named genes KN1 to KN4 were well delineated, each regrouping pine and spruce sequences. KN4 was sister group to KN3, which was sister group to KN1 and KN2. Genetic mapping on the genomes of two divergent Picea species indicated that KN1 and KN2 are located close to each other on the same linkage group, whereas KN3 and KN4 mapped on different linkage groups, correlating the more ancient divergence of these two genes. The proportion of synonymous and nonsynonymous substitutions suggested intense purifying selection for the four genes. However, rates of substitution per year indicated an evolution in two steps: faster rates were noted after gene duplications, followed subsequently by lower rates. Positive directional selection was detected for most of the internal branches harboring an accelerated rate of evolution. In addition, many sites with highly significant amino acid rate shift were identified between these branches. However, the tightly linked KN1 and KN2 did not diverge as much from each other. The implications of the correlation between phylogenetic, structural, and functional information are discussed in relation to the diversification of the knox-I gene family in conifers.
Key Words: conifer knox-I genes ? evolutionary rates ? functional divergence ? gene duplication ? genome mapping ? phylogenetic analysis
Introduction
Gene duplications are generally considered as a mechanism for increasing diversity and for functional innovation. After duplication, one gene copy may maintain the original function, whereas the other can accumulate amino acid changes and, thereby, could acquire a novel function by neofunctionalization (for a review, see Zhang [2003]). Evolutionary forces driving functional divergence between duplicated genes or paralogs can be diverse. There is increasing evidence that functional divergence is related to the translocation of gene duplicates (Lynch et al. 2001). Periods of relaxed positive evolution and evolution under strong functional constraints can also fluctuate over time (e.g., Lynch and Conery 2003). Change in the ratio of nonsynonymous to synonymous substitutions provides information for understanding the process of protein evolution (Yang and Nielsen 2000). However this approach is limited by the rapid rate of synonymous substitutions (Gaucher et al. 2002). To assess more broadly functional protein divergence among paralogs, new statistical methods that consider only amino acid replacements have emerged (Gu and Velden 2002; Knudsen et al. 2003). The central tenet of theses approaches is that the identification of rate change at nonsynonymous sites might reflect shifts in the underlying selective constraints, thus, highlighting positions that are probably involved in functional divergence (for a review, see Gaucher et al. [2002]).
Knox genes are a large group of transcription factors that belong to the homeobox gene family. KNOX proteins have a common structural organization consisting of six regions, including a conserved KNOX domain, a highly conserved ELK domain, and the homeodomain (HD) responsible for DNA binding (Ito, Sato, and Matsuoka 2002). All KNOX proteins have three conserved amino acids in the loop between helix I and helix II of HD, and, therefore, belong to the TALE (three amino acid loop extension) superfamily. Plant knox genes are defined by homology to the maize knotted1 (kn1) gene and are divided into two classes based on HD sequence similarity and expression patterns (Bharathan et al. 1999). Class II genes are found in all organs and no specific function has been ascribed to these genes. In contrast, class I genes are mainly expressed in embryos and in the shoot apical meristem (SAM), which is essential for organ formation in higher plants. Class I knox genes are known as transcriptional regulators that play an important role in plant architecture (Hake and Ori 2002). Some knox-I genes have been reported to control the development of the SAM during embryogenesis (Ito, Sato, and Matsuoka 2002). Loss-of-function mutations in the SHOOTMERISTEMLESS (STM) gene of Arabidopsis has led to the loss of embryogenic SAM, whereas other embryonic organs developed normally (reviewed in Takada and Tasaka [2002]). Thus, one of the functional attributes of knox-I genes is to maintain indeterminate cell fate and suppress determination processes within the SAM (Veit 2004). In agreement with these observations, the duplication and diversification of ancestral knox genes has resulted in gene families with related but unique functions that appeared to represent essential steps in the evolution of higher plant SAMs (Reiser, Sanchez-Baracaldo, and Hake 2000).
In conifers and more specifically in Norway spruce (Picea abies) and black spruce (Picea mariana), a total of three distinct knox-I paralogs (KN1, KN2, and KN3) have been reported from various studies (Sundas-Larsson et al. 1998; Hjortswang et al. 2002; Robert K. Rutledge, Canadian Forest Service, GenBank numbers U90091 and U90092). Phylogenetic analysis has shown that two of the paralogs (PmKN1 and PmKN2) form a monophyletic group and, thus, have evolved independently from angiosperm knox-I genes (Champagne and Ashton 2001). Recently, a distinct fourth member, KN4, has been identified from a large screening of cDNAs from diverse tissues in Picea mariana (Robert K. Rutledge, Canadian Forest Service, personal communication). Picea abies knox-I genes are expressed in embryogenic cultures, stems, roots, and cone buds but not in needles (Hjortswang et al. 2002). In Welwitschia mirabilis, knox-I genes are good indicators of meristem activity and maintenance (Pham and Sinha 2003). In addition, one of Picea abies knox-I genes, HBK2 (named PaKN3 in this study), appears to be involved in the development of somatic embryos (Hjortswang et al. 2002). Hence, the function of knox-I genes does not seem to be entirely redundant in conifers: some genes appear to be regulators during embryogenesis, and others are involved in postembryonic SAM maintenance. Therefore, it is likely that part of the diversity of conifer knox-I genes is functional and that footprints of functional divergence can be detected at the molecular level.
In this study, we used phylogenetic analysis of protein-coding sequences in combination with chromosomal localization to investigate on the expansion of knox-I gene family in conifers. Relationships with evolutionary rates were also examined to identify the tempos of evolution and the driving forces leading to the diversification of conifer knox-I genes.
Materials and Methods
Sequence Sampling
Seventy-five sequences representing class I knox genes of land plants and two sequences of class II knox genes were retrieved by screening GenBank annotations or were obtained after sequencing cDNAs or genomic fragments amplified by PCR (table 1). Only protein-coding sequences were kept for analysis after discarding the noncoding regions. For conifers, our sampling strategy was aimed at sampling and delineating the various paralogs of knox-I. Six complete cDNA sequences were retrieved or given (PmKN4), which were obtained from cDNA libraries constructed for various tissues in Picea mariana (three genes) and Picea abies (three genes) and representing a total of four distinct paralogs (KN1 to KN4). To further extend the sampling of knox-I genes in a divergent conifer taxon, the Pinus taeda EST database (http://pine.ccgb.umn.edu; Center for Computational Genomics and Bioinformatics, University of Minnesota) was screened using BlastN. Five contigs were identified, corresponding to the genes KN1 to KN3 previously found from spruce cDNAs (E-value e–48; scores ranging from 1,520 to 192). The complete cDNA sequence was further obtained for each of these pine genes by RT-PCR (see below). Using gene-specific primers designed from the alignment of spruce and pine cDNAs, genomic sequences were also obtained for additional taxa in Pinus and Picea to help delineate the four paralogs: for Pinus strobus, which represents a different subgenus than for Pinus taeda, and for Picea glauca and Picea abies, which represent with Picea mariana divergent taxa in the genus.
Table 1 List of knox I gene sequences used in this study.
DNA and RNA Extraction
DNA was isolated from individual haploid megagametophytes using a Dneasy Plant Mini Kit (Qiagen, Mississauga, Ontario). For Pinus taeda and the poplar hybrid HII-II (P. trichocarpa x P. deltoides), RNA was isolated from secondary xylem wood tissues using a CTAB procedure (Chang, Puryear, and Cairney 1993).
RT-PCR
The missing 5' and 3' ends of partial Pinus taeda cDNA clones (Center for Computational Genomics and Bioinformatics, University of Minnesota), and partial Populus trichocarpa x P. deltoides genomic contigs (International Populus Genome Consortium) were amplified by RACE using the SMART RACE cDNA Amplification Kit and the Advantage 2 PCR Kit according to the manufacturer's instructions (BD Biosciences, Palo Alto, Calif.). For RT-PCR analysis, total RNA was reverse-transcribed using the Superscript First Strand Synthesis System (Invitrogen, Carlsbad, Calif.). cDNAs were then PCR-amplified using primers specific to PtKN1 (forward, 5'-GGTTTTGAAACTAGTGTTTCATGAATATGG-3' and reverse, 5'-TGTGCTGTTTTCAGCGGCCGCTTCTTTCTA-3'), to PtKN2 (forward, 5'-GATGATGCAGTAAAGGCAGAC-3' and reverse, 5'-ACGGAAAATTACAGCTTCCCT-3'), to PtKN3 (forward, 5'- GAAATATAACTAGTTGGGATCATATG-3' and reverse, 5'-CAAATCTACAGCGGCCGCCTAAAACTG-3'), to PtdKN3 (forward, 5'-ATGGAGGACTACAATCAAATGAG-3' and reverse, 5'- TCATGGACCGAGCCGATAAT-3'), and to PtdKN2 (forward, 5'-ATGGAGGGTGGTGGTGGTGA-3' and reverse, 5'-TCAAAGCAGGGTGGGAGAGAT-3'). The complete cDNA sequence of KN4 was first determined in Picea mariana (PmKN4) by RT-PCR of midstage embryos and midstage female cones (a gift from Robert K. Rutledge, Canadian Forest Service).
Genomic PCR
Gene-specific primer pairs used to amplify spruce and pine knox-1 paralogs KN1, KN2, KN3, and KN4 are presented in tables 1 and 2, respectively, in Supplementary Material online. Because of the large size of the third intron in knox-1 genes (5.056 kb for the maize Kn1 gene [Vollbrecht et al. 1991], 5.5 kb for the rice OSH3 gene [Sato, Fukuda, and Hirano 2001], and 4.16 kb for the Arabidopsis KNAT2 gene [AC010796]), each gene was amplified in two fragments using gene-specific primers located in the more variable 5'UTR (untranslated region) and 3'UTR and designed from the alignment of complete cDNA sequences from pine and spruce (see above). For KN1, a region of 1.6 kb that encompassed the 5'UTR, the exon1, intron1, exon2, intron2, and the beginning of exon3 (5' part) and a region of 0.6 kb that encompassed the exon4, intron4, exon5, and 3'UTR (3' part) were amplified. For the other three conifer knox-I genes, the two amplicons had the same exon and intron structures with respective lengths of 2.1 kb and 0.75 kb for KN2, 1.6 kb and 0.57 kb for KN3, and 1.3 kb and 0.65 kb for KN4. PCR reactions were performed in 30 μl containing 20mM Tris-HCl (pH 8.4), 50 mM KCl, 1.5 to 2.0 mM MgCl2, 200 μM of each dNTP, 200 μM of both 5' and 3' primers, and 1.0 U of platinum Taq DNA polymerase (Invitrogen, Carlsbad, Calif.). About 5 to 20 ng of genomic DNA was used as template. A Peltier thermal cycler (DYADTM DNA Engine, MJResearch, Waltham, Mass.) was used, with the following thermal cycling profile: 4 min at 94°C, followed by 35 cycles of 30 s at 94°C, 30 s at annealing temperature optimized between 54°C and 58°C for each pair of primers, and 1 min at 72°C, followed by 10 min at 72°C. Each PCR fragment was directly sequenced in both directions with a Perkin-Elmer ABI 3730 XL DNA sequencer (Applied Biosystems, Foster City, Calif.), using BigDye Terminator cycle sequencing kits version 3.1. Contigs were constructed using Windows 32 SeqMan version 5.05 (DNASTAR Inc., Madison, Wis.) and using the BioEdit sequence alignment editor version 5.0.9 (Tom Hall, Department. of Microbiology, North Carolina State University). There was no evidence for artificial recombinants: all amplifications resulted in expected single sequencing products, alignment with known conifer cDNA sequences was congruent with expectations for each group of orthologs, and each gene was located at the same locus position on black spruce or white spruce linkage maps, whenever the 5' or the 3' part of the gene was used. Sequence data from this study (table 1) have been deposited in the GenBank database under accession numbers AY680380 to AY680405 (Pinus and Picea) and AY684937 to AY684938 (Populus trichocarpa x P. deltoides).
Table 2 dN, dS, and (dN/dS) Statistics
Phylogenetic Analyses
The 77 amino acid sequences were aligned using ClustalX (Thompson et al. 1997). ClustalX multiple-alignment parameters were gap penalty 10 and gap extension penalty 0.3, using the Gonnet250 protein weight matrix. The final alignment was then adjusted by eye using BioEdit (Hall 1999). Gene trees were inferred from amino acid sequences using the complete alignment or by considering the conserved domains only. The average length of the protein-coding sequence (without gaps) was 378 (± 51) amino acids and that of conserved domains only was 224 (± 8) amino acids. Phylogenetic relationships were inferred using maximum parsimony (MP) as implemented in PAUP* UNIX version 4.0b10 (Swofford 2002) on an Origin2000 calculator (Centre for Bioinformatics, University Laval), and the neighbor-joining (NJ) method as implemented in MEGA version 2.1 (Kumar et al. 2001). For MP analyses, heuristic searches were conducted with 100 random stepwise additions, and different analyses were conducted by considering gaps as missing data or as an additional state. For NJ, gap pairwise deletion was enforced and pairwise divergence values were corrected for multiple hits (Jukes and Cantor 1969). The robustness of tree topologies was assessed by means of 500 bootstrap replicates for NJ and 300 bootstrap replicates for the computationally more demanding MP analyses. Plant knox-I gene sequences were previously shown to be monophyletic (Champagne and Ashton 2001). Thus, in the phylogenetic analyses of the 75 land plant class I knox sequences, trees were rooted with two Arabidopsis thaliana class II knox genes here used as outgroups (AtKNAT7, GenBank number AF308451 and AtKNAT3, GenBank number X92392).
To study the evolution of the conifer group, gene trees were estimated from conifer sequences (excluding the partial Cryptomeria sequence) using MP as implemented in PAUP* UNIX version 4.0b10 (Swofford 2002). Gene trees based on DNA sequences had to be estimated for subsequent substitution rate calculations for various codon positions (see below). The aligned amino acid sequences were used as a guide for the alignment of the corresponding cDNA or genomic DNA sequences in the TRANALIGN program (EMBOSS). Then, distinct MP analyses were applied to first and second codon positions only (corresponding mostly to nonsynonymous substitutions) and to third codon positions only (corresponding mostly to synonymous substitutions). Conifer trees were rooted with two angiosperm knox-I genes (AtSTM and MtKnox) extracted from the angiosperm sister group to conifer genes, which was identified from the large-scale analysis of 75 knox-I sequences (see Results). The various trees were visualized using TREEVIEW version 1.6.6 (Page 2001; http://taxonomy.zoology.gla.ac.uk/rod/treeview.html).
Mapping knox-I Genes in Picea
The four conifer knox-I genes were mapped on the genome of Picea mariana and also on the genome of Picea glauca (white spruce) as a replicate. These two sympatric species are reproductively isolated and represent highly divergent lineages in the genus. For each species, two crosses with one parent in common were used (Picea mariana, cross numbers 9920002 and S11991V; Picea glauca, cross numbers 2856 and 2872), each consisting of 118 progeny, except 85 progeny for S11991V (table 3 in Supplementary Material online). More than 300 amplified fragment length polymorphisms (AFLPs) segregating 1:1 and 50 codominant markers (expressed sequence tag polymorphisms (ESTPs) and microsatellites), as well as the four conifer knox-I genes, were used to build each linkage map. ESTP markers of orthologous loci for Picea species have been previously developed by Pelgas, Isabel, and Bousquet (2004). The linkage analysis was performed with JoinMap version 3.0 (Van Ooijen and Voorrips 2001) using the parameter set for progeny derived by cross-pollination (CP). A minimum LOD threshold of 4.5 was used for grouping markers into linkage groups. Then, for ordering markers, a minimum LOD threshold of 3.0 and a recombination frequency of 0.35 were used. Parents and progeny were genotyped for the four conifer knox-I genes using DGGE (denaturing gradient gel electrophoresis) or CAPS (cleaved amplified polymorphic sequence) methodologies (see Pelgas, Isabel, and Bousquet [2004]), using specific primers for each gene (table 3 in Supplementary Material online). For comparison between Picea mariana and Picea glauca, at least three codominant markers per linkage group were used.
Table 3 Compositional Analysis of Conifer knox-I Genes
Compositional and Functional Diversity Analysis of knox-I Genes from Conifers
All the following analyses were based on the 19 conifer coding-region sequences (excluding the partial Cryptomeria sequence), whether at the DNA or amino acid level (table 1). Base content as well as transition/transversion ratios were estimated using MEGA version 2.1 (Kumar et al. 2001). Codon usage bias, measured using the effective number of codons (ENC, [Wright 1990]), was assessed using DnaSP version 4 (Rozas et al. 2003). Pairwise synonymous (dS) and nonsynonymous (dN) rates of substitution were estimated according to the method of Nei and Gojobori (1986), using YAM in the PAML program package (Yang 1997).
Rates of substitution per year were estimated using the penalized-likelihood method (Sanderson 2002) as implemented in the program r8s version 1.50 (Sanderson 2003). Optimal values of smoothing were determined by a cross-validation procedure implemented in the r8s program. Two distinct r8s analyses were conducted: the first one by considering the conifer MP tree constructed with first and second codon positions only (corresponding mostly to nonsynonymous substitutions) and the second one by considering the conifer MP tree constructed with third codon positions only (corresponding mostly to synonymous substitutions).
Adaptive evolution between the knox-1 conifer paralogs was investigated based on a relative rate ratio test (RRRT) developed by Creevey and McInerney (2002) as implemented in CRANN (Creevey and McInerney 2003). For each internal branch of the conifer tree, changes in the descendant clade described by that internal branch were counted. This procedure results in four types of substitutions: replacement (nonsynonymous) invariable (RI), replacement variable (RV), silent (synonymous) invariable (SI), and silent variable (SV). A significant difference between the ratio of RI to RV and SI to SV is indicative of positive selection.
To further analyze functional divergence between knox-1 conifer paralogs, sites with evolutionary rate shift involved in type-I or type-II divergence were identified using likelihood ratio tests, as implement in the program LRT (Knudsen et al. 2003). Type I sites represent amino acid configurations that are conserved in one paralog but highly variable in another; type II sites are fixed for different amino acids between paralogs (Wang and Gu 2001). LRT constructs a multiple sequence alignment with two subfamilies at a time representing two clusters of duplicated genes and infers a tree using the NJ method. The branch lengths are estimated with maximum likelihoods under the Jones, Taylor, and Thornton (JTT) model (Jones, Taylor, and Thornton 1992).
Results
Diversification of knox-1 Genes in Angiosperms and Conifers
To sample as extensively as possible the diversity of knox-I gene paralogs in dicots, monocots, and conifers, various gene searches were implemented, in addition to the new conifer sequences determined herein (table 1). In parallel to the screening of the GenBank annotations, an exhaustive search using BlastN was conducted over the complete Arabidopsis genome (TAIR database) to get all members of the knox-I family and to confirm the absence of pseudogenes. Because partial sequences of the maize knox3 and knox4 genes were retrieved from a protein database, TBlastN against the database GRAMENE (http://www.gramene.org [Ware et al. 2002]) and GenBank was also used to find their corresponding cDNA sequences. Additional BlastN searches on the TIGR Pinus Gene Index (http://www.tigr.org/tdb/tgi/pgi/; The Institute for Genomic Research) and the Picea glauca EST database (J. Mackay, J. Bousquet, et al., Arborea Project, http://www.arborea.ulaval.ca/en/results/est-sequencing.php) have not resulted in the retrieval of additional knox-I paralogs. In fact, highly significant matches (E-values e–48; scores ranging from 6,482 to 412) were obtained for contigs corresponding to the four conifer knox-I genes. The other matches were not significant (E-values > e–6). Given the large divergence between the conifer genera Pinus and Picea (140 Myr [Florin 1963; Miller 1988]), and given the vast array of tissues from which these various EST databases are derived, it is likely that the four knox-1 paralogs identified so far in pines and spruces are representative of most of the family in the conifers. BlastN searches on the 6,589 Cryptomeria japonica ESTs (Ujino-Ihara et al. 2000), available in dbEST (http://www.ncbi.nlm.nih.gov/dbEST), retrieved four ESTs corresponding to the same knox-I gene (E-values < 10–16). Given this sampling, phylogenetic trees of knox-I protein-coding regions were constructed from 31 dicot sequences, 21 monocot sequences, 21 conifer sequences, and two bryophyte sequences. Previous analyses had shown that class I and II knox genes from land plants form two monophyletic superclades (Bharathan et al. 1999; Champagne and Ashton 2001). The gene duplication event that gave rise to class I and class II knox gene superfamilies occurred before the divergence between bryophytes and spermatophytes (Champagne and Ashton 2001). Thus, two Arabidopsis knox-II genes were used to root the land plant knox-I gene tree.
The tree reported in figure 1 represents a majority-rule consensus of the 32 most parsimonious trees based on all positions of the protein sequence. The two moss sequences, PpMKN2 and PpMKN4, were sister group to the seed plants. The angiosperm sequences were divided into four large groups, each with medium to high bootstrap support from the various analytical procedures used (fig. 1). Bootstrap support for major nodes was generally consistent between MP and NJ (fig. 1). Considering sequences from conserved domains only instead of complete sequences resulted in higher or lower bootstraps, depending of the node (fig. 1). For MP analyses, considering gaps as an additional character did not result in a consistently higher or lower bootstrap support (fig. 1).
FIG. 1.— Phylogenetic analysis of knox-I genes. The tree represents the majority-rule consensus of the 32 most-parsimonious trees estimated from the amino acid sequences of 73 knox-I plant genes and by considering all positions available (consistency index [CI] = 0.611; retention index [RI] = 0.814). Bootstrap values are given for major nodes only. Bootstrap values indicated above branches were obtained from the analysis of complete sequences, whereas values below branches were from the analysis of conserved domains only. For each node, bootstrap values on the left are from maximum parsimony with gap = newstate (red); values in the middle are from parsimony with gap = missing (green); and values on the right are from neighbor-joining analysis (blue). A1 to A4 denotes four major groups of angiosperm knox-I genes. For each angiosperm group when applicable, the subdivision between dicots and monocots is indicated. D1, D2, and D3 indicate three putative duplications events in the conifer knox-I gene family. Within the conifer group, the name of each paralogous subgroup is given and P = pines and S = spruces.
The conifer genes formed a monophyletic group, with a high bootstrap support from most of the analytical procedures used (fig. 1). In addition, rooting the spermatophyte tree with class I moss sequences always resulted in a monophyletic group for conifer knox-I genes. Only MP analyses of sequences from the conserved domain, those where the more variable N-terminal and C-terminal regions were not considered, did not result in bootstrap estimates higher than 50% for this group because of the loss of many synapomorphies relative to the analysis of complete amino acid sequences. Similarly as for MP analyses of complete sequences, NJ analyses resulted in high values of bootstrap support, even when sequences from conserved domains only were considered (fig. 1). Hence, the consensus of the various analytical approaches suggests that known conifer knox-1 genes diversified and evolved independently after the split between gymnosperms and angiosperms.
Within the conifer group, the four paralogous subgroups (KN1, KN2, KN3, and KN4) contained pine and spruce sequences (fig. 1). Hence, the divergence between the four conifer subgroups occurred before the split between Pinus and Picea. The partial cDNA from Cryptomeria japonica (CjKN) could not be considered in phylogenetic analyses limited to conserved domains because it consisted of only 118 amino acid residues in the C-terminal region. Nevertheless, in the other analyses, it clustered with the conifer group but did not cluster with any of the four KN subgroups, presumably because this sequence was partial (fig. 1).
KN1 and KN2 showed a high degree of amino acid identity (average of 85%), whereas KN4 and KN3 paralogs were less similar to KN1 and KN2 and between each other (from 58% to 68% identity). Altogether, the diversification pattern suggests three successive events of duplication (D1 to D3) occurring during a period not exceeding 160 Myr from the split between angiosperms and gymnosperms (300 Myr [Savard et al. 1994]) to the divergence of conifer genera (140 Myr [Florin 1963; Miller 1988]).
Assignment of the Four Conifer knox-1 Genes to Linkage Groups
The four conifer knox-I paralogs were mapped on the genome of Picea mariana and Picea glauca. From one individual linkage map to another within each species, and between Picea mariana and Picea glauca, the order of codominant markers remained consistent within linkage groups. Detailed results on the linkage maps developed from the four progeny arrays will be reported elsewhere (B. Pelgas, S. Beauseigle, J. Bousquet, and N. Isabel, unpublished data). The less divergent KN1 and KN2 were assigned to the same linkage group and localized close to each other, whereas KN3 and KN4 mapped on different linkage groups (fig. 2), paralleling their more remote phylogenetic position.
FIG. 2.— Assignment of the four conifer knox-1 genes (KN1, KN2, KN3, and KN4) to three different linkage groups in Picea. Linkage groups represented by vertical bars were derived from segregation data of the Picea mariana cross number 9920002. Markers in boldfaced type are codominant markers used as anchor markers for intraspecific and interspecific comparisons. Marker types are indicated by the last letter in the locus names: m = microsatellites and e = ESTP. Microsatellite markers are designated by primer pairs. ESTP markers were denoted according to Pelgas, Isabel, and Bousquet (2004), except SEM017 (Lamothe and Isabel, unpublished data). AFLP markers listed in plain types are named with primer combination followed by fragment length. Distances between markers in centiMorgans (Haldane's function) are indicated to the left of each linkage group.
Accelerated Evolution of the Conifer knox-1 Genes
Within the knox-1 conifer group, a notable heterogeneity of branch lengths was detected on the various phylogenetic trees obtained. Long branches were noted during the period of gene diversification, and short branches were observed after the split between Pinus and Picea (data not shown). To further assess this pattern, rates of substitution per site per year (absolute rates) were estimated, and their variability was monitored across the conifer gene phylogeny. An age was assigned or estimated for each internal node to permit the estimation of rates of substitution per unit of time (Sanderson 2002, 2003). An age of 300 Myr was assigned to the "root" node corresponding to the split between angiosperms and gymnosperms (Savard et al. 1994) (fig. 3). The four nodes defining the four paralogous gene clusters (KN1, KN2, KN3, and KN4) were assigned an age of 140 Myr (fig. 3), which corresponds to the diversification of the Pinaceae (Florin 1963; Miller 1988). Thus, we constrained the nodes D1, D2, and D3 (which correspond to the three duplications, [fig. 1]) to a minimum age of 140 Myr and a maximum age of 300 Myr. The nodes defining splits between sequences of subgenera haploxylon and diploxylon of Pinus (fig. 3) were assigned an age of 100 Myr (Prager, Fowler, and Wilson 1976). Finally, the internal nodes corresponding to sequence divergence between Picea lineages were given an age of 15 Myr (Bouillé and Bousquet 2004). The divergence times of each duplication event (D1, D2, and D3) were estimated by r8s analyses (Sanderson 2002, 2003). The resulting two tree topologies are shown as ratograms, where branch lengths correspond directly to rates of change per site per year (Sanderson 2003) (fig. 3). Thus, on a ratogram, branches can be directly compared for their rates per site per year, a longer branch meaning a higher rate of change per site per year. If the rate of substitution per site along the branch is needed, it can simply be obtained by multiplying the rate per site per year by the number of years corresponding to that branch (fig. 3).
FIG. 3.— The knox-I conifer tree with branch lengths representing rates of substitution per site per year (ratogram [Sanderson 2003]). The topology of both trees is derived from maximum parsimony. (A) Ratogram constructed with first and second codon positions. (B) ratogram constructed with third codon positions. The node "root" was fixed at 300 Myr and asteriks indicate the basal node for each gene cluster KN1, KN2, KN3, and KN4, which was fixed at 140 Myr (see text). Internal branches, which were found under positive directional selection according to the relative rate ratio test of Creevey and McInerney (2002), are indicated in bold (numbered 3, 7, 15, 16, and 17) (see text). For branch #17, a rate of 1.80 x 10–9 substitution/site/year was estimated over a period of 58 Myr; for branch #16, a rate of 1.63 x 10–9 substitution/site/year was estimated over a period of 29 Myr; for branch #15, a rate of 1.36 x 10–9 substitution/site/year was estimated over a period of 53 Myr; for branch #7, a rate of 0.98 x 10–9 substitution/site/year was estimated over a period of 73 Myr; for branch #3, a rate of 0.81 x 10–9 substitution/site/year was estimated over a period of 102 Myr. Branch #11 (a rate of 0.75 x 10–9 substitution/site/year estimated over a period of 20 Myr) and branch #14 (a rate of 0.69 x 10–9 substitution/site/year estimated over a period of 20 Myr) were not found under positive directional selection (see table 4). For branches after the split between Pinus and Picea (see asterisks), the average rate of substitution/site/year over the four knox-I gene clusters was 0.094 x 10–9 for first and second codon positions.
Table 4 Significant Results from Relative Rate Ratio Tests of Functional Divergence Among the Conifer knox-I Genes
Both ratograms (fig. 3A and B) indicate an evolution in two steps for conifer knox-1 genes. From the root node of the four paralogs (300 Myr) to the nodes corresponding to the split between Pinus and Picea (140 Myr), a period that includes the three duplication events, accelerated rates per site per year were observed, as indicated by the long branches on the ratograms (fig. 3). For that period, the estimated rates ranged from 0.69 x 10–9 to 1.80 x 10–9 substitutions per site per year, with a mean of 1.15 x 10–9 for first and second codon positions (fig. 3A), and from 1.34 x 10–9 to 3.20 x 10–9 substitutions per site per year, with a mean of 2.27 x 10–9 for third codon positions (fig. 3B).
In contrast, after the split between Pinus and Picea (140 Myr, see asterisks in fig. 3), the rates per site per year slowed down dramatically for all four paralogs, as indicated by the very short branches on the ratograms (fig. 3), with mean values of 0.094 x 10–9 (12-fold difference) for first and second codon positions and 0.385 x 10–9 (sixfold difference) for third codon positions. A correlated response was observed between the two categories of sites because a rate decrease of approximately the same amplitude was observed for third codon positions defining mostly synonymous changes versus first and second codon positions defining mostly nonsynonymous changes (fig. 3).
Rates of Synonymous and Nonsynonymous Substitution in Conifer knox-I Genes
The above estimates of substitution rate per site per year indicated an early period of faster evolution, followed by a period of slower evolution for the four conifer knox-I genes. To investigate whether this temporal heterogeneity in evolutionary rates was caused by selection, synonymous (dS) and nonsynonymous (dN) rates of substitution per site and the ratio dN/dS () were evaluated for all pairwise comparisons of conifer sequences (table 2). Between paralogs, the estimated dN values and the estimated dS values had a mean of 0.271 and 2.003, respectively. Between orthologs, the estimated dN values had a mean of 0.020, and the estimated dS values averaged 0.135. For all pairwise comparisons, the estimated dS values were seven times higher, on average, than dN values. Both dN and dS values were, on average, 14 times higher in pairwise comparisons of paralogues than between orthologous gene pairs within each conifer knox-I gene cluster, indicating a larger degree of sequence divergence between the four conifer paralogs. The coefficient of variation of rates among all conifer knox-I genes was 0.58 for dN and 0.63 for dS, indicating that synonymous and nonsynonymous rates had varied about to the same extant among conifer genes. The dN/dS ratios () measured between conifer paralogs had a mean of 0.143 and were similar to those between orthologous gene pairs within each conifer knox-I gene cluster, which averaged 0.130. A near homogenous among conifer knox-I genes also indicates that both synonymous and nonsynonymous sites evolved in a correlated fashion, as seen from r8s analyses. Overall, the ratios indicated strong functional constraints (purifying selection, < 1).
To further characterize the patterns of sequence divergence between duplicated genes, we examined factors other than selection that could contribute to variation in rates of synonymous and nonsynonymous substitution among genes. As GC content and codon usage bias have been shown to affect rates of synonymous substitution (e.g., Zhang, Vision, and Gaut 2002), we estimated the effective number of codons (ENCs), the GC content, and the transition/tranversion ratio for each gene. GC content was similar among paralogs. However, the ENCs were lower for KN4, the %GC at third codon positions was lower for KN3, and both genes had a less skewed transition/tranversion ratio, as compared with that of the closely linked KN1 and KN2 (table 3).
Adaptive Evolution
We used the relative rate ratio test (RRRT) developed by Creevey and Mclnerney (2002) to detect whether knox-I conifer genes underwent a period of adaptive evolution. We found far more replacement invariant (RI) substitutions, those nonsynonymous substitutions where the new character state is preserved in all subsequent lineages, than would be expected from neutrality alone for a number of internal branches of the conifer knox-I gene tree (table 4). These branches are all located before the divergence of Pinus from Picea, 140 MYA (fig. 3A). These branches include that from the conifer knox-1 "root" node to node D1 (#17), from the node D1 to D2 (#16), from the node D2 to D3 (#15), from the node D1 to the KN4 clade (#3), and from the node D2 to the KN3 clade (#7). Although the estimation of failed to detect positive selection, such an excess of RI over RV pattern should be taken as evidence for positive directional selection and neofunctionalization. Theses branches were also characterized by an accelerated rate of change per site per year in the r8s analyses (fig. 3A). No significant differences between the RI/RV and SI/SV ratios were noted for the branches leading to the KN1 clade (#11) or the KN2 clade (#14) (table 4).
The likelihood ratio tests (LRTs) developed by Knudsen et al. (2003) were used to detect significant rate shift at specific amino acid sites (table 5). The first duplication, that between KN4 and the other three gene paralogs, was analyzed by comparing amino acid changes in the multiple sequence alignment constructed with KN4 genes (here considered as subfamily 1), and KN1, KN2, and KN3 genes (here considered as subfamily 2). The LRTs identified 15 sites with U 2, those considered to have highly significant rate shift (Knudsen et al. 2003), and some of them were located in conserved functional domains (table 5). The second duplication was analyzed by comparing the amino acid changes between KN3 genes (here considered as subfamily 1) and KN1 and KN2 genes (here considered as subfamily 2). The LRTs identified 5 sites with U 2, with some of them in conserved functional domains (table 5). The third duplication was analyzed by comparing the amino acid sequences of KN1 orthologs to those of KN2 orthologs. The LRTs identified only one site with U 2, and it was located outside the conserved functional domains (table 5). These results correlate those from the RRRTs.
Table 5 Number of Amino Acid Sites with Significant Rate Shift Among Conifer knox-I Genes Using the Likelihood Ratio
Discussion
knox-1 Genes in Spermatophytes
The knox-I gene family has diversified in at least four large groups early in the evolution of spermatophytes. All these groups (A1, A2, A3, and A4) contained angiosperm sequences. The conifer group sister to A4 contained all the known conifer sequences, which have diversified into four subgroups apparently not shared by angiosperms. The knox-I gene family was also found to be composed of four paralogs in Arabidopsis (Semiarti et al. 2001). However, they did not form a monophyletic group and were remote from each other (see fig. 1). Polyphyly and large dispersion was also observed for other angiosperm taxa well represented in the sampling such as Zea, Oryza, Lycopersicon, and Helianthus (see fig. 1). Polyphyly was also observed for the conifer genera Pinus and Picea for which a diversity of tissues was sampled for cDNAs. However, the conifer group as a whole remained cohesive. This pattern suggests that knox-I genes have pursued a distinct path of evolution in conifers, presumably by gaining new paralogs after the split of their common ancestor from that of angiosperms. In comparison with other families of plant transcription factors, such a distinct evolution from angiosperms appears to be unusual. Families of plant transcription factors such as MYBR2R3 (Xue et al. 2003), HD-GL2 (Ingouff et al. 2003), HD-ZIP class III (Guillet-Claude et al., unpublished data) and diverse classes of MADS-Box genes (Nam et al. 2003) failed to reveal a monophyletic assemblage for conifer genes. However, monophyletic assemblages were observed in Larix for the genes coding for 4-coumarate:coenzyme A ligase (4CL) and the genes coding for cinnamyl alcohol dehydrogenase (CAD) (Wei and Wang 2004).
For three of the four angiosperm groups delineated by phylogenetic analyses, there were no conifer sister groups. Given that such groups were expected, this observation suggests incomplete sampling of conifer genes or gene loss. However, given the numerous clones analyzed from cDNA libraries screened from diverse tissues in Picea mariana and given the various pine and spruce EST databases screened herein (see Results), it is possible that the absence of conifer genes in some of the angiosperm groups observed is not a sampling artifact and that gene loss might have happened in conifers after new paralogs were gained. The birth and death of genes is a common theme in gene family evolution (Nei, Rogozin, and Piontkivska 2000), and the majority of genes appear to be lost after duplication (e.g. Kellis, Biren, and Lander 2004). The absence of monocot sequences in some of the knox-I angiosperm groups (fig. 1) may represent another putative case of gene loss (see figure 1), as reported by Tioni, Gonzales, and Chan (2003). It is also possible that other knox-I genes exist in conifers and in monocots, which have not been sampled because of transient or low levels of expression (Robert G. Rutledge, personal communication). It could be argued that gene loss in certain taxa would be more likely if the various groups of spermatophyte knox-I genes had high functional redundancy. Recently, sequence comparisons, expression studies, and mutant analysis indicated that the knox-I genes AtKNAT1, OsH15, ZmRS1, and ZmKnox4, all present in the group A3 (fig. 1), may encode equivalent knox functions in grasses and dicots (Veit 2004). However, evidence for functional redundancy between phylogenetically distant knox-I genes in angiosperms remains contradictory. For instance, it has been recently shown that AtSTM, an Arabidopsis knox-I gene from the group A4 (fig. 1), regulates AtKNAT1 (group A3) and AtKNAT2 (group A2), thus, suggesting divergent roles in regulating meristem function among these three angiosperm knox-I groups (Byrne, Simorowski, and Martienssen 2002). At the same time, mutant analyses revealed partial redundancy among Arabidopsis genes. AtKNAT1 was shown to assume a redundant role with AtSTM in the vegetative shoot apical meristem but could not compensate for AtSTM in floral meristem (Byrne, Simorowski, and Martienssen 2002). More detailed functional studies will be necessary before any firm conclusion can be reached on the redundancy of function among knox-I paralogs.
Mechanisms for the Expansion of the knox-1 Gene Family in Conifers
The phylogenetic structure of knox-I genes from conifers suggests that three duplication events (D1, D2, and D3 [fig. 1]) occurred between the split of the lineage leading to conifers from that leading to angioperms, around 300 MYA (Savard et al. 1994), and the more recent divergence between Pinus and Picea, here indicated by the diversification of the Pinaceae, around 140 MYA (Florin 1963; Miller 1988). KN1 and KN2, which are the result of a more recent duplication, were assigned to the same linkage group and localized close to each other. On the other hand, KN3 and KN4 mapped on different linkage groups (fig. 2), paralleling their more remote locations on the phylogenetic tree and their increased divergence at the sequence level. The remote chromosomal location of KN3 and KN4 would involve local duplication followed by a secondary mechanism bringing about rearrangements and translocation of chromosomal segments to a distant location. Such a model of structural evolution has been proposed for the homeobox-leucine zipper genes in Arabidopsis (Schena and Davis 1994). Comparative genomics studies also indicated that gene duplication is often accompanied by genetic map changes in all genomes surveyed (e.g., Riechmann et al. 2000; Lynch et al. 2001). Our results are coherent with the analysis of the Arabidopsis homeobox gene family, which showed that gene duplicates located on different chromosomes are more frequent ( 71%) than those on the same ones (Riechmann et al. 2000). Perhaps this trend is reflective of the relative abundance of ancient versus more recent gene duplication events.
Functional Divergence Among Conifer knox-I Genes
It has been shown early on that after gene duplication, one gene duplicate maintain the original function, whereas the other duplicate is free to accumulate amino acid changes as a result of functional redundancy and positive directional selection (Ohno 1970; Li 1983). Recent studies have indicated a more diverse array of potential outcomes in terms of fates of duplicated genes: pseudogeneization, functional redundancy, subfunctionalization in which each daughter gene adopts part of the functions of their parental gene, and neofunctionalization in which there is gain of novel function (for reviews, see Lawton-Rauh [2003] and Zhang [2003]). The differences between the last three categories are not so obvious in the absence of structural or biochemical data at the protein level. In our study, evidence for purifying selection was detected at the molecular level, with ratios of nonsynonymous (dN) to synonymous (dS) substitutions much smaller than 1 between sequences within each conifer knox-I gene as well as between conifer knox-I genes (dS values six times higher than dN values on average) (table 2). At their face value, these estimates indicate little evidence for directional selection. A similar trend towards low values (dS values five times higher than dN values) was observed between gene duplicates in Arabidopsis (Zhang, Vision, and Gaut 2002).
However, the ratio is based on rate averages since the split between two gene sequences. It does not take into account possible temporal heterogeneity in rates of evolution since the split between two gene duplicates, which could be indicative of periods of relaxed selection, followed by periods of more stringent selection. Such a rate shift is suggested by recent theoretical studies showing that duplicate genes experience two phases of evolution: a rather short period of relaxed selection and, thereafter, a period of purifying selection (Lynch and Conery 2003). When estimating rates of substitution per year for the diverse branches of the conifer tree (fig. 3), signatures of positive selection and functional divergence were evident for internal branches leading to the four conifer knox-I gene clades. R8s analyses showed that the calibrated rates of evolution per year were on average nine times larger before the split between Pinus and Picea than after, 140 MYA. Using the RRRT of Creevey and McInerney (2002), significant positive directional selection was also detected for several internal branches after the duplication of conifer knox-I genes. Such nonneutral divergence, early in the history of the conifer knox-I genes, is likely the result of a reduction in selective constraints on at least one duplicate, together with a commensurate increase in the amino acid replacement rate (Zhang 2003). At the opposite, short branches at the tip of the tree after the split between Pinus and Picea are indicative of a slow down in rates of evolution for the four conifer knox-I genes (fig. 3), which would imply more intense purifying selection in the recent past.
Relation Between Functional Divergence and Change in Map Position
The different chromosomal locations observed for the conifer knox-I paralogs were related to increased functional divergence. In our study, three unlinked locations were disclosed, corresponding respectively to (1) the genes KN1 and KN2, (2) the gene KN3, and (3) the gene KN4. These distinct locations parallel the significant results from RRRTs for positive selection and the large number of sites with significant rate shift from LRTs among the three groups. On the basis of model-based predictions, Lynch et al. (2001) suggested that the probability of change in map position per newly arisen gene duplicate increases with the strength of selection on neofunctional alleles. Similarly, the analysis of the genomic organization of multigene families in Arabidopsis suggested that the movement of gene duplicates to unlinked chromosomal regions may play an important role in functional divergence because of less exposure to sequence homogenization (Leister 2004).
The conifer paralogs KN1 and KN2 were positioned close to each other on the same linkage group, and there was little indication of significant positive directional selection between them, which is suggestive of more or less complete redundancy of function. At the same time, they presented a higher degree of sequence identity (average of 85% at the amino acid level), they had similar effective numbers of codon, %GC at third codon positions and transition/transversion ratios (table 3). It might be argued that this greater homogeneity is the result of a more recent gene duplication. However, the close physical relationship between KN1 and KN2 might have contributed to this greater homogeneity by exposing them to sequence homogenization effects. Indeed, the rates of substitution per site per year were lower for the branches leading to the KN1 or the KN2 clade (#11 and #14), as compared with those leading to the KN3 or KN4 clade (#7 and #3) and, more notably, to the branches representing the common ancestor of KN1 and KN2 (#15) and that of KN1, KN2, and KN3 (#16) (fig. 3).
Correlated Rates of Evolution Between Synonymous and Nonsynonymous Sites
A near homogenous was observed among conifer knox-I genes, indicating correlated rates of synonymous and nonsynonymous substitutions. Correlated temporal heterogeneity in rates of evolution was also observed between rates of change per year for first and second codon positions (mostly nonsynonymous) and rates of change per year for third codon positions (mostly synonymous): longer internal branches were noted before the divergence between Pinus, and Picea and shorter branches we noted after the divergence between Pinus and Picea (fig. 3). These trends are suggestive of factors affecting both synonymous and nonsynonymous sites in a similar way. A similar correlated trend between dS and dN was reported from the study of 10 nuclear conifer genes (Kusumi et al. 2002) and from the study of 242 duplicated genes from chromosomes 2 and 4 of Arabidopsis thaliana (Zhang, Vision, and Gaut 2002). This pattern could be indicative of background selection, which implies the removal of deleterious mutations as well as linked neutral or nearly neutral substitutions (Charlesworth, Morgan, and Charlesworth 1993). It would result in a decrease in the proportion of neutral mutations and a consequent correlated decrease in the rates of synonymous and nonsynonymous substitutions (Ohta 1995). Such a correlated selective constraint on both synonymous and nonsynonymous sites has also been observed in mammalian genes (Mouchiroud, Gautier, and Bernardi 1995).
This study demonstrates that knox-I genes did not evolve in a linear fashion in conifers. Initial periods of relaxed evolution have been followed by a period of conservatism, at least for the past 140 Myr. Taxon sampling in additional orders of the Gymnosperms would help delineate the beginning of this conservative period. A number of hypotheses at the functional level also need to be addressed by gene expression studies and investigations on biochemical properties and protein structure. For instance, it is anticipated that KN1 and KN2 would show more redundancy or overlapping in their physiological roles than with KN3 and KN4, which roles should be more distinct. These studies should assist in validating the putative functional importance of the changes observed at the molecular level and the extent of subfunctionalization or neofunctionalization among these genes.
Supplementary Material
The following information is available online: a list of primers used for the genomic amplification of spruce knox-I genes (table 1 of the Supplementary Material online) and pine knox-I genes (table 2 of the Supplementary Material online) and additional information about the mapping protocol of the four conifer knox-I genes in Picea mariana and Picea glauca (table 3 of the Supplementary Material online).
Acknowledgements
This work was carried out as part of the Arborea project lead by J. Mackay (Université Laval) and was supported by grants from Genome Québec and Genome Canada, the Natural Sciences and Engineering Research Council of Canada (NSERC, genomic program), the Canadian Biotechnology Strategy and the Canadian Forest Service. We thank R. K. Rutledge (Canadian Forest Service, Laurentian Forestry Centre, Québec) for the gift of PmKN4 from Picea mariana and for his advice and critical reading of a previous draft of this manuscript. We also thank F. Gagnon, S. Blais, and M. Ouellet for their assistance in the laboratory and their contribution to the RACE and RT-PCR experiments and two anonymous reviewers and B. Gaut for their insightful comments, which helped improve the manuscript.
References
Bharathan, G., B. J. Janssen, E. A. Kellogg, and N. Sinha. 1999. Phylogenetic relationships and evolution of the KNOTTED class of plant homeodomain proteins. Mol. Biol. Evol. 16:553–563.
Bouillé, M., and J. Bousquet. 2004. Trans-species shared polymorphisms at orthologous nuclear gene loci among distant species in the conifer Picea: implications for the maintenance of genetic diversity in trees. Am. J. Bot. (in press).
Byrne, M. E., J. Simorowski, and R. A. Martienssen. 2002. ASYMMETRIC LEAVES1 reveals knox gene redundancy in Arabidopsis. Development 129:1957–1965.
Champagne, C. E. M., and N. W. Ashton. 2001. Ancestry of KNOX genes revealed by bryophyte (Physcomytrella patens) homologs. New Phytol. 150:23–36.
Chang, S., J. Puryear, and J. Cairney. 1993. A simple and efficient method for isolating RNA from pine trees. Plant Mol. Biol. Rep. 11:113–116.
Charlesworth, B., M. T. Morgan, and D. Charlesworth. 1993. The effect of deleterious mutations on neutral molecular variation. Genetics 134:1289–1303.
Creevey, C. J., and J. O. McInerney. 2002. An algorithm for detecting directional and non-directional positive selection, neutrality and negative selection in protein coding DNA sequences. Gene 300:43–51.
Creevey, C. J., and J. O. McInerney. 2003. CRANN: detecting adaptive evolution in protein-coding DNA sequences. Bioinformatics 19:1726.
Florin, R. 1963. The distribution of conifer and taxad genera in time and space. Acta Horti. Bergiani 20:121–312.
Gaucher, E. A., X. Gu, M. M. Miyamoto, and S. A. Benner. 2002. Predicting functional divergence in protein evolution by site-specific rate shifts. Trends Biochem. Sci. 27:315–321.
Gu, X., and K. Vander Velden. 2002. DIVERGE: phylogeny-based analysis for functional-structural divergence of a protein family. Bioinformatics 18:500–501.
Hake, S., and N. Ori. 2002. Plant morphogenesis and KNOX genes. Nat. Genet. 31:121–122.
Hall, T. A. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acids Symp. Ser. 41:95–98.
Hjortswang, H. I., A. Sundas-Larsson, G. Bharathan, P. V. Bozhkov, S. von Arnold, and T. Vahala. 2002. KNOTTED1-like homeobox genes of gymnosperm, Norway spruce, expressed during somatic embryogenesis. Plant Physiol. Biochem. 40:837–843.
Ingouff, M., I. Farbos, M. Wiweger, and S. von Arnold. 2003. The molecular characterization of PaHB2, a homeobox gene of the HD-GL2 family expressed during embryo development in Norway spruce. J. Exp. Bot. 54:1343–1350.
Ito, M., Y. Sato, and M. Matsuoka. 2002. Involvement of homeobox genes in early body plan of monocot. Int. Rev. Cytol. 218:1–35.
Jones, D. T., W. R. Taylor, and J. M. Thornton. 1992. The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci. 8:275–282.
Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21–132 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.
Kellis, M., B. W. Birren, and E. S. Lander. 2004. Proof and evolutionary analysis of ancient genome duplication in the yeast Saccharomyces cerevisiae. Nature 428:617–624.
Knudsen, B., M. M. Miyamoto, P. J. Laipis, and D. N. Silverman. 2003. Using evolutionary rates to investigate protein functional divergence and conservation: a case study of the carbonic anhydrases. Genetics 164:1261–1269.
Kumar, S., K. Tamura, I. B. Jakobsen, and M. Nei. 2001. MEGA2: molecular evolutionary genetics analysis software. Bioinformatics 17:1244–1245.
Kusumi, J., Y. Tsumura, H. Yoshimaru, and H. Tachida. 2002. Molecular evolution of nuclear genes in Cupressacea, a group of conifer trees. Mol. Biol. Evol. 19:736–747.
Lawton-Rauh, A. 2003. Evolutionary dynamics of duplicated genes in plants. Mol. Phylogenet. Evol. 29:396–409.
Leister, D. 2004. Tandem and segmental gene duplication and recombination in the evolution of plant disease resistance genes. Trends Genet. 20:116–122.
Li, W. H. 1983. Evolution of duplicated genes. Pp. 14–37 in M. Nei and R. K. Koehin, eds. Evolution of genes and proteins. Sinauer, Sunderland, Mass.
Lynch, M., and J. S. Conery. 2003. The evolutionary demography of duplicate genes. J. Struct. Funct. Genomics 3:35–44.
Lynch, M., M. O'Hely, B. Walsh, and A. Force. 2001. The probability of preservation of a newly arisen gene duplicate. Genetics 159:1789–1804.
Miller, C. N. 1988. The origin of modern conifer families. Pp. 448–486 in C. B. Beck, ed. Origin and evolution of Gymnosperms. Columbia University Press, New York.
Mouchiroud, D., C. Gautier, and G. Bernardi. 1995. Frequencies of synonymous substitutions in mammals are gene-specific and correlated with frequencies of nonsynonymous substitutions. J. Mol. Evol. 40:107–113.
Nam, J., C. W. dePamphilis, H. Ma, and M. Nei. 2003. Antiquity and evolution of the MADS-box gene family controlling flower development in plants. Mol. Biol. Evol. 20:1435–1447.
Nei, M., and T. Gojobori. 1986. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol. Biol. Evol. 3:418–426.
Nei, M., I. B. Rogozin, and H. Piontkivska. 2000. Purifying selection and birth-and-death evolution in the ubiquitin gene family. Proc. Natl. Acad. Sci. USA 97:10866–10871.
Ohno, S. 1970. Evolution by gene duplication. Allen and Unwin, London.
Ohta, T. 1995. Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral theory. J. Mol. Evol. 40:56–63.
Pelgas, B., N. Isabel, and J. Bousquet. 2004. Efficient screening for expressed sequence tag polymorphisms (ESTP) by DNA pool sequencing and denaturing gradient gel electrophoresis in spruces. Mol. Breed. 13:263–279.
Pham, T., and N. Sinha. 2003. Role of KNOX genes in shoot development of Welwitschia mirabilis. Int. J. Plant Sci. 164:333–343.
Prager, E. M., D. P. Fowler, and A. C. Wilson. 1976. Rates of evolution in conifers (pinaceae). Evolution 30:637–649.
Reiser, L., P. Sanchez-Baracaldo, and S. Hake. 2000. Knots in the family tree: evolutionary relationships and functions of knox homeobox genes. Plant Mol. Biol. 42:151–166.
Riechmann, J. L., J. Heard, G. Martin et al. (17 co-authors). 2000. Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science 290:2105–2110.
Rozas, J., J. C. Sanchez-DelBarrio, X. Messeguer, and R. Rozas. 2003. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19:2496–2497.
Sanderson, M. J. 2002. Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol. Biol. Evol. 19:101–109.
———. 2003. r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. Bioinformatics 19:301–302.
Sato, Y., Y. Fukuda, and H. Y. Hirano. 2001. Mutations that cause amino acid substitutions at the invariant positions in homeodomain of OSH3.KNOX protein suggest artificial selection during rice domestication. Genes Genet. Syst. 76:381–392.
Savard, L., P. Li, S. H. Strauss, M. W. Chase, M. Michaud, and J. Bousquet. 1994. Chloroplast and nuclear gene sequences indicate late Pennsylvanian time for the last common ancestor of extant seed plants. Proc. Natl. Acad. Sci. USA 91:5163–5167.
Schena, M., and R. W. Davis. 1994. Structure of homeobox-leucine zipper genes suggests a model for the evolution of gene families. Proc. Natl. Acad. Sci. USA 91:8393–8397.
Semiarti, E., Y. Ueno, H. Tsukaya, H. Iwakawa, C. Machida, and Y. Machida. 2001. The ASYMMETRIC LEAVES2 gene of Arabidopsis thaliana regulates formation of a symmetric lamina, establishment of venation and repression of meristem-related homeobox genes in leaves. Development 128:1771–1783.
Sundas-Larsson, A., M. Svenson, H. Liao, and P. Engstrom. 1998. A homeobox gene with potential developmental control function in the meristem of the conifer Picea abies. Proc. Natl. Acad. Sci. USA 95:15118–15122.
Swofford, D. L. 2002. PAUP*: phylogeneric analysis using parsimony (* and other methods). Sinauer Associates, Sunderland, Mass.
Takada, S., and M. Tasaka. 2002. Embryonic shoot apical meristem formation in higher plants. J. Plant Res. 115:411–417.
Thompson, J. D., T. J. Gibson, F. Plewniak, F. Jeanmougin, and D. G. Higgins. 1997. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25:4876–4882.
Tioni, M. F., D. H. Gonzalez, and R. L. Chan. 2003. Knotted1-like genes are strongly expressed in differentiated cell types in sunflower. J. Exp. Bot. 54:681–690.
Ujino-Ihara, T., K. Yoshimura, Y. Ugawa, H. Yoshimaru, K. Nagasaka, and Y. Tsumura. 2000. Expression analysis of ESTs derived from the inner bark of Cryptomeria japonica. Plant Mol. Biol. 43:451–457.
Van Ooijen, J. W., and R. E. Voorrips. 2001. JoinMap: Software for the calculation of genetic linkage maps. Version 3.0. Plant Research International, Wageningen, The Netherlands.
Veit, B. 2004. Determination of cell fate in apical meristems. Curr. Opin. Plant Biol. 7:57–64.
Vollbrecht, E., B. Veit, N. Sinha, and S. Hake. 1991. The developmental gene Knotted-1 is a member of a maize homeobox gene family. Nature 350:241–243.
Wang, Y., and X. Gu. 2001. Functional divergence in the caspase gene family and altered functional constraints: statistical analysis and prediction. Genetics 158:1311–1320.
Ware, D., P. Jaiswal, J. Ni et al. (12 co-authors). 2002. Gramene: a resource for comparative grass genomics. Nucleic Acids Res. 30:103–105.
Wei, X. X., and X. Q. Wang. 2004. Evolution of 4-coumarate:coenzyme A ligase (4CL) gene and divergence of Larix (Pinaceae). Mol. Phylogenet. Evol. 31:542–553.
Wright, F. 1990. The ‘effective number of codons’ used in a gene. Gene 87:23–29.
Xue, B., P. J. Charest, Y. Devantier, and R. G. Rutledge. 2003. Characterization of a MYBR2R3 gene from black spruce (Picea mariana) that shares functional conservation with maize C1. Mol. Genet. Genomics 270:78–86.
Yang, Z. 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. 13:555–556.
Yang, Z., and R. Nielsen. 2000. Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17:32–43.
Zhang, J. 2003. Evolution by gene duplication: an update. Trends Ecol. Evol. 18:292–298.
Zhang, L., T. J. Vision, and B. S. Gaut. 2002. Patterns of nucleotide substitution among simultaneously duplicated gene pairs in Arabidopsis thaliana. Mol. Biol. Evol. 19:1464–1473.(Carine Guillet-Claude*, N)
E-mail: bousquet@rsvs.ulaval.ca.
Abstract
Class I knox genes code for transcription factors that play an essential role in plant growth and development as central regulators of meristem cell identity. Based on the analysis of new cDNA sequences from various tissues and genomic DNA sequences, we identified a highly diversified group of class I knox genes in conifers. Phylogenetic analyses of complete amino acid sequences from various seed plants indicated that all conifer sequences formed a monophyletic group. Within conifers, four subgroups here named genes KN1 to KN4 were well delineated, each regrouping pine and spruce sequences. KN4 was sister group to KN3, which was sister group to KN1 and KN2. Genetic mapping on the genomes of two divergent Picea species indicated that KN1 and KN2 are located close to each other on the same linkage group, whereas KN3 and KN4 mapped on different linkage groups, correlating the more ancient divergence of these two genes. The proportion of synonymous and nonsynonymous substitutions suggested intense purifying selection for the four genes. However, rates of substitution per year indicated an evolution in two steps: faster rates were noted after gene duplications, followed subsequently by lower rates. Positive directional selection was detected for most of the internal branches harboring an accelerated rate of evolution. In addition, many sites with highly significant amino acid rate shift were identified between these branches. However, the tightly linked KN1 and KN2 did not diverge as much from each other. The implications of the correlation between phylogenetic, structural, and functional information are discussed in relation to the diversification of the knox-I gene family in conifers.
Key Words: conifer knox-I genes ? evolutionary rates ? functional divergence ? gene duplication ? genome mapping ? phylogenetic analysis
Introduction
Gene duplications are generally considered as a mechanism for increasing diversity and for functional innovation. After duplication, one gene copy may maintain the original function, whereas the other can accumulate amino acid changes and, thereby, could acquire a novel function by neofunctionalization (for a review, see Zhang [2003]). Evolutionary forces driving functional divergence between duplicated genes or paralogs can be diverse. There is increasing evidence that functional divergence is related to the translocation of gene duplicates (Lynch et al. 2001). Periods of relaxed positive evolution and evolution under strong functional constraints can also fluctuate over time (e.g., Lynch and Conery 2003). Change in the ratio of nonsynonymous to synonymous substitutions provides information for understanding the process of protein evolution (Yang and Nielsen 2000). However this approach is limited by the rapid rate of synonymous substitutions (Gaucher et al. 2002). To assess more broadly functional protein divergence among paralogs, new statistical methods that consider only amino acid replacements have emerged (Gu and Velden 2002; Knudsen et al. 2003). The central tenet of theses approaches is that the identification of rate change at nonsynonymous sites might reflect shifts in the underlying selective constraints, thus, highlighting positions that are probably involved in functional divergence (for a review, see Gaucher et al. [2002]).
Knox genes are a large group of transcription factors that belong to the homeobox gene family. KNOX proteins have a common structural organization consisting of six regions, including a conserved KNOX domain, a highly conserved ELK domain, and the homeodomain (HD) responsible for DNA binding (Ito, Sato, and Matsuoka 2002). All KNOX proteins have three conserved amino acids in the loop between helix I and helix II of HD, and, therefore, belong to the TALE (three amino acid loop extension) superfamily. Plant knox genes are defined by homology to the maize knotted1 (kn1) gene and are divided into two classes based on HD sequence similarity and expression patterns (Bharathan et al. 1999). Class II genes are found in all organs and no specific function has been ascribed to these genes. In contrast, class I genes are mainly expressed in embryos and in the shoot apical meristem (SAM), which is essential for organ formation in higher plants. Class I knox genes are known as transcriptional regulators that play an important role in plant architecture (Hake and Ori 2002). Some knox-I genes have been reported to control the development of the SAM during embryogenesis (Ito, Sato, and Matsuoka 2002). Loss-of-function mutations in the SHOOTMERISTEMLESS (STM) gene of Arabidopsis has led to the loss of embryogenic SAM, whereas other embryonic organs developed normally (reviewed in Takada and Tasaka [2002]). Thus, one of the functional attributes of knox-I genes is to maintain indeterminate cell fate and suppress determination processes within the SAM (Veit 2004). In agreement with these observations, the duplication and diversification of ancestral knox genes has resulted in gene families with related but unique functions that appeared to represent essential steps in the evolution of higher plant SAMs (Reiser, Sanchez-Baracaldo, and Hake 2000).
In conifers and more specifically in Norway spruce (Picea abies) and black spruce (Picea mariana), a total of three distinct knox-I paralogs (KN1, KN2, and KN3) have been reported from various studies (Sundas-Larsson et al. 1998; Hjortswang et al. 2002; Robert K. Rutledge, Canadian Forest Service, GenBank numbers U90091 and U90092). Phylogenetic analysis has shown that two of the paralogs (PmKN1 and PmKN2) form a monophyletic group and, thus, have evolved independently from angiosperm knox-I genes (Champagne and Ashton 2001). Recently, a distinct fourth member, KN4, has been identified from a large screening of cDNAs from diverse tissues in Picea mariana (Robert K. Rutledge, Canadian Forest Service, personal communication). Picea abies knox-I genes are expressed in embryogenic cultures, stems, roots, and cone buds but not in needles (Hjortswang et al. 2002). In Welwitschia mirabilis, knox-I genes are good indicators of meristem activity and maintenance (Pham and Sinha 2003). In addition, one of Picea abies knox-I genes, HBK2 (named PaKN3 in this study), appears to be involved in the development of somatic embryos (Hjortswang et al. 2002). Hence, the function of knox-I genes does not seem to be entirely redundant in conifers: some genes appear to be regulators during embryogenesis, and others are involved in postembryonic SAM maintenance. Therefore, it is likely that part of the diversity of conifer knox-I genes is functional and that footprints of functional divergence can be detected at the molecular level.
In this study, we used phylogenetic analysis of protein-coding sequences in combination with chromosomal localization to investigate on the expansion of knox-I gene family in conifers. Relationships with evolutionary rates were also examined to identify the tempos of evolution and the driving forces leading to the diversification of conifer knox-I genes.
Materials and Methods
Sequence Sampling
Seventy-five sequences representing class I knox genes of land plants and two sequences of class II knox genes were retrieved by screening GenBank annotations or were obtained after sequencing cDNAs or genomic fragments amplified by PCR (table 1). Only protein-coding sequences were kept for analysis after discarding the noncoding regions. For conifers, our sampling strategy was aimed at sampling and delineating the various paralogs of knox-I. Six complete cDNA sequences were retrieved or given (PmKN4), which were obtained from cDNA libraries constructed for various tissues in Picea mariana (three genes) and Picea abies (three genes) and representing a total of four distinct paralogs (KN1 to KN4). To further extend the sampling of knox-I genes in a divergent conifer taxon, the Pinus taeda EST database (http://pine.ccgb.umn.edu; Center for Computational Genomics and Bioinformatics, University of Minnesota) was screened using BlastN. Five contigs were identified, corresponding to the genes KN1 to KN3 previously found from spruce cDNAs (E-value e–48; scores ranging from 1,520 to 192). The complete cDNA sequence was further obtained for each of these pine genes by RT-PCR (see below). Using gene-specific primers designed from the alignment of spruce and pine cDNAs, genomic sequences were also obtained for additional taxa in Pinus and Picea to help delineate the four paralogs: for Pinus strobus, which represents a different subgenus than for Pinus taeda, and for Picea glauca and Picea abies, which represent with Picea mariana divergent taxa in the genus.
Table 1 List of knox I gene sequences used in this study.
DNA and RNA Extraction
DNA was isolated from individual haploid megagametophytes using a Dneasy Plant Mini Kit (Qiagen, Mississauga, Ontario). For Pinus taeda and the poplar hybrid HII-II (P. trichocarpa x P. deltoides), RNA was isolated from secondary xylem wood tissues using a CTAB procedure (Chang, Puryear, and Cairney 1993).
RT-PCR
The missing 5' and 3' ends of partial Pinus taeda cDNA clones (Center for Computational Genomics and Bioinformatics, University of Minnesota), and partial Populus trichocarpa x P. deltoides genomic contigs (International Populus Genome Consortium) were amplified by RACE using the SMART RACE cDNA Amplification Kit and the Advantage 2 PCR Kit according to the manufacturer's instructions (BD Biosciences, Palo Alto, Calif.). For RT-PCR analysis, total RNA was reverse-transcribed using the Superscript First Strand Synthesis System (Invitrogen, Carlsbad, Calif.). cDNAs were then PCR-amplified using primers specific to PtKN1 (forward, 5'-GGTTTTGAAACTAGTGTTTCATGAATATGG-3' and reverse, 5'-TGTGCTGTTTTCAGCGGCCGCTTCTTTCTA-3'), to PtKN2 (forward, 5'-GATGATGCAGTAAAGGCAGAC-3' and reverse, 5'-ACGGAAAATTACAGCTTCCCT-3'), to PtKN3 (forward, 5'- GAAATATAACTAGTTGGGATCATATG-3' and reverse, 5'-CAAATCTACAGCGGCCGCCTAAAACTG-3'), to PtdKN3 (forward, 5'-ATGGAGGACTACAATCAAATGAG-3' and reverse, 5'- TCATGGACCGAGCCGATAAT-3'), and to PtdKN2 (forward, 5'-ATGGAGGGTGGTGGTGGTGA-3' and reverse, 5'-TCAAAGCAGGGTGGGAGAGAT-3'). The complete cDNA sequence of KN4 was first determined in Picea mariana (PmKN4) by RT-PCR of midstage embryos and midstage female cones (a gift from Robert K. Rutledge, Canadian Forest Service).
Genomic PCR
Gene-specific primer pairs used to amplify spruce and pine knox-1 paralogs KN1, KN2, KN3, and KN4 are presented in tables 1 and 2, respectively, in Supplementary Material online. Because of the large size of the third intron in knox-1 genes (5.056 kb for the maize Kn1 gene [Vollbrecht et al. 1991], 5.5 kb for the rice OSH3 gene [Sato, Fukuda, and Hirano 2001], and 4.16 kb for the Arabidopsis KNAT2 gene [AC010796]), each gene was amplified in two fragments using gene-specific primers located in the more variable 5'UTR (untranslated region) and 3'UTR and designed from the alignment of complete cDNA sequences from pine and spruce (see above). For KN1, a region of 1.6 kb that encompassed the 5'UTR, the exon1, intron1, exon2, intron2, and the beginning of exon3 (5' part) and a region of 0.6 kb that encompassed the exon4, intron4, exon5, and 3'UTR (3' part) were amplified. For the other three conifer knox-I genes, the two amplicons had the same exon and intron structures with respective lengths of 2.1 kb and 0.75 kb for KN2, 1.6 kb and 0.57 kb for KN3, and 1.3 kb and 0.65 kb for KN4. PCR reactions were performed in 30 μl containing 20mM Tris-HCl (pH 8.4), 50 mM KCl, 1.5 to 2.0 mM MgCl2, 200 μM of each dNTP, 200 μM of both 5' and 3' primers, and 1.0 U of platinum Taq DNA polymerase (Invitrogen, Carlsbad, Calif.). About 5 to 20 ng of genomic DNA was used as template. A Peltier thermal cycler (DYADTM DNA Engine, MJResearch, Waltham, Mass.) was used, with the following thermal cycling profile: 4 min at 94°C, followed by 35 cycles of 30 s at 94°C, 30 s at annealing temperature optimized between 54°C and 58°C for each pair of primers, and 1 min at 72°C, followed by 10 min at 72°C. Each PCR fragment was directly sequenced in both directions with a Perkin-Elmer ABI 3730 XL DNA sequencer (Applied Biosystems, Foster City, Calif.), using BigDye Terminator cycle sequencing kits version 3.1. Contigs were constructed using Windows 32 SeqMan version 5.05 (DNASTAR Inc., Madison, Wis.) and using the BioEdit sequence alignment editor version 5.0.9 (Tom Hall, Department. of Microbiology, North Carolina State University). There was no evidence for artificial recombinants: all amplifications resulted in expected single sequencing products, alignment with known conifer cDNA sequences was congruent with expectations for each group of orthologs, and each gene was located at the same locus position on black spruce or white spruce linkage maps, whenever the 5' or the 3' part of the gene was used. Sequence data from this study (table 1) have been deposited in the GenBank database under accession numbers AY680380 to AY680405 (Pinus and Picea) and AY684937 to AY684938 (Populus trichocarpa x P. deltoides).
Table 2 dN, dS, and (dN/dS) Statistics
Phylogenetic Analyses
The 77 amino acid sequences were aligned using ClustalX (Thompson et al. 1997). ClustalX multiple-alignment parameters were gap penalty 10 and gap extension penalty 0.3, using the Gonnet250 protein weight matrix. The final alignment was then adjusted by eye using BioEdit (Hall 1999). Gene trees were inferred from amino acid sequences using the complete alignment or by considering the conserved domains only. The average length of the protein-coding sequence (without gaps) was 378 (± 51) amino acids and that of conserved domains only was 224 (± 8) amino acids. Phylogenetic relationships were inferred using maximum parsimony (MP) as implemented in PAUP* UNIX version 4.0b10 (Swofford 2002) on an Origin2000 calculator (Centre for Bioinformatics, University Laval), and the neighbor-joining (NJ) method as implemented in MEGA version 2.1 (Kumar et al. 2001). For MP analyses, heuristic searches were conducted with 100 random stepwise additions, and different analyses were conducted by considering gaps as missing data or as an additional state. For NJ, gap pairwise deletion was enforced and pairwise divergence values were corrected for multiple hits (Jukes and Cantor 1969). The robustness of tree topologies was assessed by means of 500 bootstrap replicates for NJ and 300 bootstrap replicates for the computationally more demanding MP analyses. Plant knox-I gene sequences were previously shown to be monophyletic (Champagne and Ashton 2001). Thus, in the phylogenetic analyses of the 75 land plant class I knox sequences, trees were rooted with two Arabidopsis thaliana class II knox genes here used as outgroups (AtKNAT7, GenBank number AF308451 and AtKNAT3, GenBank number X92392).
To study the evolution of the conifer group, gene trees were estimated from conifer sequences (excluding the partial Cryptomeria sequence) using MP as implemented in PAUP* UNIX version 4.0b10 (Swofford 2002). Gene trees based on DNA sequences had to be estimated for subsequent substitution rate calculations for various codon positions (see below). The aligned amino acid sequences were used as a guide for the alignment of the corresponding cDNA or genomic DNA sequences in the TRANALIGN program (EMBOSS). Then, distinct MP analyses were applied to first and second codon positions only (corresponding mostly to nonsynonymous substitutions) and to third codon positions only (corresponding mostly to synonymous substitutions). Conifer trees were rooted with two angiosperm knox-I genes (AtSTM and MtKnox) extracted from the angiosperm sister group to conifer genes, which was identified from the large-scale analysis of 75 knox-I sequences (see Results). The various trees were visualized using TREEVIEW version 1.6.6 (Page 2001; http://taxonomy.zoology.gla.ac.uk/rod/treeview.html).
Mapping knox-I Genes in Picea
The four conifer knox-I genes were mapped on the genome of Picea mariana and also on the genome of Picea glauca (white spruce) as a replicate. These two sympatric species are reproductively isolated and represent highly divergent lineages in the genus. For each species, two crosses with one parent in common were used (Picea mariana, cross numbers 9920002 and S11991V; Picea glauca, cross numbers 2856 and 2872), each consisting of 118 progeny, except 85 progeny for S11991V (table 3 in Supplementary Material online). More than 300 amplified fragment length polymorphisms (AFLPs) segregating 1:1 and 50 codominant markers (expressed sequence tag polymorphisms (ESTPs) and microsatellites), as well as the four conifer knox-I genes, were used to build each linkage map. ESTP markers of orthologous loci for Picea species have been previously developed by Pelgas, Isabel, and Bousquet (2004). The linkage analysis was performed with JoinMap version 3.0 (Van Ooijen and Voorrips 2001) using the parameter set for progeny derived by cross-pollination (CP). A minimum LOD threshold of 4.5 was used for grouping markers into linkage groups. Then, for ordering markers, a minimum LOD threshold of 3.0 and a recombination frequency of 0.35 were used. Parents and progeny were genotyped for the four conifer knox-I genes using DGGE (denaturing gradient gel electrophoresis) or CAPS (cleaved amplified polymorphic sequence) methodologies (see Pelgas, Isabel, and Bousquet [2004]), using specific primers for each gene (table 3 in Supplementary Material online). For comparison between Picea mariana and Picea glauca, at least three codominant markers per linkage group were used.
Table 3 Compositional Analysis of Conifer knox-I Genes
Compositional and Functional Diversity Analysis of knox-I Genes from Conifers
All the following analyses were based on the 19 conifer coding-region sequences (excluding the partial Cryptomeria sequence), whether at the DNA or amino acid level (table 1). Base content as well as transition/transversion ratios were estimated using MEGA version 2.1 (Kumar et al. 2001). Codon usage bias, measured using the effective number of codons (ENC, [Wright 1990]), was assessed using DnaSP version 4 (Rozas et al. 2003). Pairwise synonymous (dS) and nonsynonymous (dN) rates of substitution were estimated according to the method of Nei and Gojobori (1986), using YAM in the PAML program package (Yang 1997).
Rates of substitution per year were estimated using the penalized-likelihood method (Sanderson 2002) as implemented in the program r8s version 1.50 (Sanderson 2003). Optimal values of smoothing were determined by a cross-validation procedure implemented in the r8s program. Two distinct r8s analyses were conducted: the first one by considering the conifer MP tree constructed with first and second codon positions only (corresponding mostly to nonsynonymous substitutions) and the second one by considering the conifer MP tree constructed with third codon positions only (corresponding mostly to synonymous substitutions).
Adaptive evolution between the knox-1 conifer paralogs was investigated based on a relative rate ratio test (RRRT) developed by Creevey and McInerney (2002) as implemented in CRANN (Creevey and McInerney 2003). For each internal branch of the conifer tree, changes in the descendant clade described by that internal branch were counted. This procedure results in four types of substitutions: replacement (nonsynonymous) invariable (RI), replacement variable (RV), silent (synonymous) invariable (SI), and silent variable (SV). A significant difference between the ratio of RI to RV and SI to SV is indicative of positive selection.
To further analyze functional divergence between knox-1 conifer paralogs, sites with evolutionary rate shift involved in type-I or type-II divergence were identified using likelihood ratio tests, as implement in the program LRT (Knudsen et al. 2003). Type I sites represent amino acid configurations that are conserved in one paralog but highly variable in another; type II sites are fixed for different amino acids between paralogs (Wang and Gu 2001). LRT constructs a multiple sequence alignment with two subfamilies at a time representing two clusters of duplicated genes and infers a tree using the NJ method. The branch lengths are estimated with maximum likelihoods under the Jones, Taylor, and Thornton (JTT) model (Jones, Taylor, and Thornton 1992).
Results
Diversification of knox-1 Genes in Angiosperms and Conifers
To sample as extensively as possible the diversity of knox-I gene paralogs in dicots, monocots, and conifers, various gene searches were implemented, in addition to the new conifer sequences determined herein (table 1). In parallel to the screening of the GenBank annotations, an exhaustive search using BlastN was conducted over the complete Arabidopsis genome (TAIR database) to get all members of the knox-I family and to confirm the absence of pseudogenes. Because partial sequences of the maize knox3 and knox4 genes were retrieved from a protein database, TBlastN against the database GRAMENE (http://www.gramene.org [Ware et al. 2002]) and GenBank was also used to find their corresponding cDNA sequences. Additional BlastN searches on the TIGR Pinus Gene Index (http://www.tigr.org/tdb/tgi/pgi/; The Institute for Genomic Research) and the Picea glauca EST database (J. Mackay, J. Bousquet, et al., Arborea Project, http://www.arborea.ulaval.ca/en/results/est-sequencing.php) have not resulted in the retrieval of additional knox-I paralogs. In fact, highly significant matches (E-values e–48; scores ranging from 6,482 to 412) were obtained for contigs corresponding to the four conifer knox-I genes. The other matches were not significant (E-values > e–6). Given the large divergence between the conifer genera Pinus and Picea (140 Myr [Florin 1963; Miller 1988]), and given the vast array of tissues from which these various EST databases are derived, it is likely that the four knox-1 paralogs identified so far in pines and spruces are representative of most of the family in the conifers. BlastN searches on the 6,589 Cryptomeria japonica ESTs (Ujino-Ihara et al. 2000), available in dbEST (http://www.ncbi.nlm.nih.gov/dbEST), retrieved four ESTs corresponding to the same knox-I gene (E-values < 10–16). Given this sampling, phylogenetic trees of knox-I protein-coding regions were constructed from 31 dicot sequences, 21 monocot sequences, 21 conifer sequences, and two bryophyte sequences. Previous analyses had shown that class I and II knox genes from land plants form two monophyletic superclades (Bharathan et al. 1999; Champagne and Ashton 2001). The gene duplication event that gave rise to class I and class II knox gene superfamilies occurred before the divergence between bryophytes and spermatophytes (Champagne and Ashton 2001). Thus, two Arabidopsis knox-II genes were used to root the land plant knox-I gene tree.
The tree reported in figure 1 represents a majority-rule consensus of the 32 most parsimonious trees based on all positions of the protein sequence. The two moss sequences, PpMKN2 and PpMKN4, were sister group to the seed plants. The angiosperm sequences were divided into four large groups, each with medium to high bootstrap support from the various analytical procedures used (fig. 1). Bootstrap support for major nodes was generally consistent between MP and NJ (fig. 1). Considering sequences from conserved domains only instead of complete sequences resulted in higher or lower bootstraps, depending of the node (fig. 1). For MP analyses, considering gaps as an additional character did not result in a consistently higher or lower bootstrap support (fig. 1).
FIG. 1.— Phylogenetic analysis of knox-I genes. The tree represents the majority-rule consensus of the 32 most-parsimonious trees estimated from the amino acid sequences of 73 knox-I plant genes and by considering all positions available (consistency index [CI] = 0.611; retention index [RI] = 0.814). Bootstrap values are given for major nodes only. Bootstrap values indicated above branches were obtained from the analysis of complete sequences, whereas values below branches were from the analysis of conserved domains only. For each node, bootstrap values on the left are from maximum parsimony with gap = newstate (red); values in the middle are from parsimony with gap = missing (green); and values on the right are from neighbor-joining analysis (blue). A1 to A4 denotes four major groups of angiosperm knox-I genes. For each angiosperm group when applicable, the subdivision between dicots and monocots is indicated. D1, D2, and D3 indicate three putative duplications events in the conifer knox-I gene family. Within the conifer group, the name of each paralogous subgroup is given and P = pines and S = spruces.
The conifer genes formed a monophyletic group, with a high bootstrap support from most of the analytical procedures used (fig. 1). In addition, rooting the spermatophyte tree with class I moss sequences always resulted in a monophyletic group for conifer knox-I genes. Only MP analyses of sequences from the conserved domain, those where the more variable N-terminal and C-terminal regions were not considered, did not result in bootstrap estimates higher than 50% for this group because of the loss of many synapomorphies relative to the analysis of complete amino acid sequences. Similarly as for MP analyses of complete sequences, NJ analyses resulted in high values of bootstrap support, even when sequences from conserved domains only were considered (fig. 1). Hence, the consensus of the various analytical approaches suggests that known conifer knox-1 genes diversified and evolved independently after the split between gymnosperms and angiosperms.
Within the conifer group, the four paralogous subgroups (KN1, KN2, KN3, and KN4) contained pine and spruce sequences (fig. 1). Hence, the divergence between the four conifer subgroups occurred before the split between Pinus and Picea. The partial cDNA from Cryptomeria japonica (CjKN) could not be considered in phylogenetic analyses limited to conserved domains because it consisted of only 118 amino acid residues in the C-terminal region. Nevertheless, in the other analyses, it clustered with the conifer group but did not cluster with any of the four KN subgroups, presumably because this sequence was partial (fig. 1).
KN1 and KN2 showed a high degree of amino acid identity (average of 85%), whereas KN4 and KN3 paralogs were less similar to KN1 and KN2 and between each other (from 58% to 68% identity). Altogether, the diversification pattern suggests three successive events of duplication (D1 to D3) occurring during a period not exceeding 160 Myr from the split between angiosperms and gymnosperms (300 Myr [Savard et al. 1994]) to the divergence of conifer genera (140 Myr [Florin 1963; Miller 1988]).
Assignment of the Four Conifer knox-1 Genes to Linkage Groups
The four conifer knox-I paralogs were mapped on the genome of Picea mariana and Picea glauca. From one individual linkage map to another within each species, and between Picea mariana and Picea glauca, the order of codominant markers remained consistent within linkage groups. Detailed results on the linkage maps developed from the four progeny arrays will be reported elsewhere (B. Pelgas, S. Beauseigle, J. Bousquet, and N. Isabel, unpublished data). The less divergent KN1 and KN2 were assigned to the same linkage group and localized close to each other, whereas KN3 and KN4 mapped on different linkage groups (fig. 2), paralleling their more remote phylogenetic position.
FIG. 2.— Assignment of the four conifer knox-1 genes (KN1, KN2, KN3, and KN4) to three different linkage groups in Picea. Linkage groups represented by vertical bars were derived from segregation data of the Picea mariana cross number 9920002. Markers in boldfaced type are codominant markers used as anchor markers for intraspecific and interspecific comparisons. Marker types are indicated by the last letter in the locus names: m = microsatellites and e = ESTP. Microsatellite markers are designated by primer pairs. ESTP markers were denoted according to Pelgas, Isabel, and Bousquet (2004), except SEM017 (Lamothe and Isabel, unpublished data). AFLP markers listed in plain types are named with primer combination followed by fragment length. Distances between markers in centiMorgans (Haldane's function) are indicated to the left of each linkage group.
Accelerated Evolution of the Conifer knox-1 Genes
Within the knox-1 conifer group, a notable heterogeneity of branch lengths was detected on the various phylogenetic trees obtained. Long branches were noted during the period of gene diversification, and short branches were observed after the split between Pinus and Picea (data not shown). To further assess this pattern, rates of substitution per site per year (absolute rates) were estimated, and their variability was monitored across the conifer gene phylogeny. An age was assigned or estimated for each internal node to permit the estimation of rates of substitution per unit of time (Sanderson 2002, 2003). An age of 300 Myr was assigned to the "root" node corresponding to the split between angiosperms and gymnosperms (Savard et al. 1994) (fig. 3). The four nodes defining the four paralogous gene clusters (KN1, KN2, KN3, and KN4) were assigned an age of 140 Myr (fig. 3), which corresponds to the diversification of the Pinaceae (Florin 1963; Miller 1988). Thus, we constrained the nodes D1, D2, and D3 (which correspond to the three duplications, [fig. 1]) to a minimum age of 140 Myr and a maximum age of 300 Myr. The nodes defining splits between sequences of subgenera haploxylon and diploxylon of Pinus (fig. 3) were assigned an age of 100 Myr (Prager, Fowler, and Wilson 1976). Finally, the internal nodes corresponding to sequence divergence between Picea lineages were given an age of 15 Myr (Bouillé and Bousquet 2004). The divergence times of each duplication event (D1, D2, and D3) were estimated by r8s analyses (Sanderson 2002, 2003). The resulting two tree topologies are shown as ratograms, where branch lengths correspond directly to rates of change per site per year (Sanderson 2003) (fig. 3). Thus, on a ratogram, branches can be directly compared for their rates per site per year, a longer branch meaning a higher rate of change per site per year. If the rate of substitution per site along the branch is needed, it can simply be obtained by multiplying the rate per site per year by the number of years corresponding to that branch (fig. 3).
FIG. 3.— The knox-I conifer tree with branch lengths representing rates of substitution per site per year (ratogram [Sanderson 2003]). The topology of both trees is derived from maximum parsimony. (A) Ratogram constructed with first and second codon positions. (B) ratogram constructed with third codon positions. The node "root" was fixed at 300 Myr and asteriks indicate the basal node for each gene cluster KN1, KN2, KN3, and KN4, which was fixed at 140 Myr (see text). Internal branches, which were found under positive directional selection according to the relative rate ratio test of Creevey and McInerney (2002), are indicated in bold (numbered 3, 7, 15, 16, and 17) (see text). For branch #17, a rate of 1.80 x 10–9 substitution/site/year was estimated over a period of 58 Myr; for branch #16, a rate of 1.63 x 10–9 substitution/site/year was estimated over a period of 29 Myr; for branch #15, a rate of 1.36 x 10–9 substitution/site/year was estimated over a period of 53 Myr; for branch #7, a rate of 0.98 x 10–9 substitution/site/year was estimated over a period of 73 Myr; for branch #3, a rate of 0.81 x 10–9 substitution/site/year was estimated over a period of 102 Myr. Branch #11 (a rate of 0.75 x 10–9 substitution/site/year estimated over a period of 20 Myr) and branch #14 (a rate of 0.69 x 10–9 substitution/site/year estimated over a period of 20 Myr) were not found under positive directional selection (see table 4). For branches after the split between Pinus and Picea (see asterisks), the average rate of substitution/site/year over the four knox-I gene clusters was 0.094 x 10–9 for first and second codon positions.
Table 4 Significant Results from Relative Rate Ratio Tests of Functional Divergence Among the Conifer knox-I Genes
Both ratograms (fig. 3A and B) indicate an evolution in two steps for conifer knox-1 genes. From the root node of the four paralogs (300 Myr) to the nodes corresponding to the split between Pinus and Picea (140 Myr), a period that includes the three duplication events, accelerated rates per site per year were observed, as indicated by the long branches on the ratograms (fig. 3). For that period, the estimated rates ranged from 0.69 x 10–9 to 1.80 x 10–9 substitutions per site per year, with a mean of 1.15 x 10–9 for first and second codon positions (fig. 3A), and from 1.34 x 10–9 to 3.20 x 10–9 substitutions per site per year, with a mean of 2.27 x 10–9 for third codon positions (fig. 3B).
In contrast, after the split between Pinus and Picea (140 Myr, see asterisks in fig. 3), the rates per site per year slowed down dramatically for all four paralogs, as indicated by the very short branches on the ratograms (fig. 3), with mean values of 0.094 x 10–9 (12-fold difference) for first and second codon positions and 0.385 x 10–9 (sixfold difference) for third codon positions. A correlated response was observed between the two categories of sites because a rate decrease of approximately the same amplitude was observed for third codon positions defining mostly synonymous changes versus first and second codon positions defining mostly nonsynonymous changes (fig. 3).
Rates of Synonymous and Nonsynonymous Substitution in Conifer knox-I Genes
The above estimates of substitution rate per site per year indicated an early period of faster evolution, followed by a period of slower evolution for the four conifer knox-I genes. To investigate whether this temporal heterogeneity in evolutionary rates was caused by selection, synonymous (dS) and nonsynonymous (dN) rates of substitution per site and the ratio dN/dS () were evaluated for all pairwise comparisons of conifer sequences (table 2). Between paralogs, the estimated dN values and the estimated dS values had a mean of 0.271 and 2.003, respectively. Between orthologs, the estimated dN values had a mean of 0.020, and the estimated dS values averaged 0.135. For all pairwise comparisons, the estimated dS values were seven times higher, on average, than dN values. Both dN and dS values were, on average, 14 times higher in pairwise comparisons of paralogues than between orthologous gene pairs within each conifer knox-I gene cluster, indicating a larger degree of sequence divergence between the four conifer paralogs. The coefficient of variation of rates among all conifer knox-I genes was 0.58 for dN and 0.63 for dS, indicating that synonymous and nonsynonymous rates had varied about to the same extant among conifer genes. The dN/dS ratios () measured between conifer paralogs had a mean of 0.143 and were similar to those between orthologous gene pairs within each conifer knox-I gene cluster, which averaged 0.130. A near homogenous among conifer knox-I genes also indicates that both synonymous and nonsynonymous sites evolved in a correlated fashion, as seen from r8s analyses. Overall, the ratios indicated strong functional constraints (purifying selection, < 1).
To further characterize the patterns of sequence divergence between duplicated genes, we examined factors other than selection that could contribute to variation in rates of synonymous and nonsynonymous substitution among genes. As GC content and codon usage bias have been shown to affect rates of synonymous substitution (e.g., Zhang, Vision, and Gaut 2002), we estimated the effective number of codons (ENCs), the GC content, and the transition/tranversion ratio for each gene. GC content was similar among paralogs. However, the ENCs were lower for KN4, the %GC at third codon positions was lower for KN3, and both genes had a less skewed transition/tranversion ratio, as compared with that of the closely linked KN1 and KN2 (table 3).
Adaptive Evolution
We used the relative rate ratio test (RRRT) developed by Creevey and Mclnerney (2002) to detect whether knox-I conifer genes underwent a period of adaptive evolution. We found far more replacement invariant (RI) substitutions, those nonsynonymous substitutions where the new character state is preserved in all subsequent lineages, than would be expected from neutrality alone for a number of internal branches of the conifer knox-I gene tree (table 4). These branches are all located before the divergence of Pinus from Picea, 140 MYA (fig. 3A). These branches include that from the conifer knox-1 "root" node to node D1 (#17), from the node D1 to D2 (#16), from the node D2 to D3 (#15), from the node D1 to the KN4 clade (#3), and from the node D2 to the KN3 clade (#7). Although the estimation of failed to detect positive selection, such an excess of RI over RV pattern should be taken as evidence for positive directional selection and neofunctionalization. Theses branches were also characterized by an accelerated rate of change per site per year in the r8s analyses (fig. 3A). No significant differences between the RI/RV and SI/SV ratios were noted for the branches leading to the KN1 clade (#11) or the KN2 clade (#14) (table 4).
The likelihood ratio tests (LRTs) developed by Knudsen et al. (2003) were used to detect significant rate shift at specific amino acid sites (table 5). The first duplication, that between KN4 and the other three gene paralogs, was analyzed by comparing amino acid changes in the multiple sequence alignment constructed with KN4 genes (here considered as subfamily 1), and KN1, KN2, and KN3 genes (here considered as subfamily 2). The LRTs identified 15 sites with U 2, those considered to have highly significant rate shift (Knudsen et al. 2003), and some of them were located in conserved functional domains (table 5). The second duplication was analyzed by comparing the amino acid changes between KN3 genes (here considered as subfamily 1) and KN1 and KN2 genes (here considered as subfamily 2). The LRTs identified 5 sites with U 2, with some of them in conserved functional domains (table 5). The third duplication was analyzed by comparing the amino acid sequences of KN1 orthologs to those of KN2 orthologs. The LRTs identified only one site with U 2, and it was located outside the conserved functional domains (table 5). These results correlate those from the RRRTs.
Table 5 Number of Amino Acid Sites with Significant Rate Shift Among Conifer knox-I Genes Using the Likelihood Ratio
Discussion
knox-1 Genes in Spermatophytes
The knox-I gene family has diversified in at least four large groups early in the evolution of spermatophytes. All these groups (A1, A2, A3, and A4) contained angiosperm sequences. The conifer group sister to A4 contained all the known conifer sequences, which have diversified into four subgroups apparently not shared by angiosperms. The knox-I gene family was also found to be composed of four paralogs in Arabidopsis (Semiarti et al. 2001). However, they did not form a monophyletic group and were remote from each other (see fig. 1). Polyphyly and large dispersion was also observed for other angiosperm taxa well represented in the sampling such as Zea, Oryza, Lycopersicon, and Helianthus (see fig. 1). Polyphyly was also observed for the conifer genera Pinus and Picea for which a diversity of tissues was sampled for cDNAs. However, the conifer group as a whole remained cohesive. This pattern suggests that knox-I genes have pursued a distinct path of evolution in conifers, presumably by gaining new paralogs after the split of their common ancestor from that of angiosperms. In comparison with other families of plant transcription factors, such a distinct evolution from angiosperms appears to be unusual. Families of plant transcription factors such as MYBR2R3 (Xue et al. 2003), HD-GL2 (Ingouff et al. 2003), HD-ZIP class III (Guillet-Claude et al., unpublished data) and diverse classes of MADS-Box genes (Nam et al. 2003) failed to reveal a monophyletic assemblage for conifer genes. However, monophyletic assemblages were observed in Larix for the genes coding for 4-coumarate:coenzyme A ligase (4CL) and the genes coding for cinnamyl alcohol dehydrogenase (CAD) (Wei and Wang 2004).
For three of the four angiosperm groups delineated by phylogenetic analyses, there were no conifer sister groups. Given that such groups were expected, this observation suggests incomplete sampling of conifer genes or gene loss. However, given the numerous clones analyzed from cDNA libraries screened from diverse tissues in Picea mariana and given the various pine and spruce EST databases screened herein (see Results), it is possible that the absence of conifer genes in some of the angiosperm groups observed is not a sampling artifact and that gene loss might have happened in conifers after new paralogs were gained. The birth and death of genes is a common theme in gene family evolution (Nei, Rogozin, and Piontkivska 2000), and the majority of genes appear to be lost after duplication (e.g. Kellis, Biren, and Lander 2004). The absence of monocot sequences in some of the knox-I angiosperm groups (fig. 1) may represent another putative case of gene loss (see figure 1), as reported by Tioni, Gonzales, and Chan (2003). It is also possible that other knox-I genes exist in conifers and in monocots, which have not been sampled because of transient or low levels of expression (Robert G. Rutledge, personal communication). It could be argued that gene loss in certain taxa would be more likely if the various groups of spermatophyte knox-I genes had high functional redundancy. Recently, sequence comparisons, expression studies, and mutant analysis indicated that the knox-I genes AtKNAT1, OsH15, ZmRS1, and ZmKnox4, all present in the group A3 (fig. 1), may encode equivalent knox functions in grasses and dicots (Veit 2004). However, evidence for functional redundancy between phylogenetically distant knox-I genes in angiosperms remains contradictory. For instance, it has been recently shown that AtSTM, an Arabidopsis knox-I gene from the group A4 (fig. 1), regulates AtKNAT1 (group A3) and AtKNAT2 (group A2), thus, suggesting divergent roles in regulating meristem function among these three angiosperm knox-I groups (Byrne, Simorowski, and Martienssen 2002). At the same time, mutant analyses revealed partial redundancy among Arabidopsis genes. AtKNAT1 was shown to assume a redundant role with AtSTM in the vegetative shoot apical meristem but could not compensate for AtSTM in floral meristem (Byrne, Simorowski, and Martienssen 2002). More detailed functional studies will be necessary before any firm conclusion can be reached on the redundancy of function among knox-I paralogs.
Mechanisms for the Expansion of the knox-1 Gene Family in Conifers
The phylogenetic structure of knox-I genes from conifers suggests that three duplication events (D1, D2, and D3 [fig. 1]) occurred between the split of the lineage leading to conifers from that leading to angioperms, around 300 MYA (Savard et al. 1994), and the more recent divergence between Pinus and Picea, here indicated by the diversification of the Pinaceae, around 140 MYA (Florin 1963; Miller 1988). KN1 and KN2, which are the result of a more recent duplication, were assigned to the same linkage group and localized close to each other. On the other hand, KN3 and KN4 mapped on different linkage groups (fig. 2), paralleling their more remote locations on the phylogenetic tree and their increased divergence at the sequence level. The remote chromosomal location of KN3 and KN4 would involve local duplication followed by a secondary mechanism bringing about rearrangements and translocation of chromosomal segments to a distant location. Such a model of structural evolution has been proposed for the homeobox-leucine zipper genes in Arabidopsis (Schena and Davis 1994). Comparative genomics studies also indicated that gene duplication is often accompanied by genetic map changes in all genomes surveyed (e.g., Riechmann et al. 2000; Lynch et al. 2001). Our results are coherent with the analysis of the Arabidopsis homeobox gene family, which showed that gene duplicates located on different chromosomes are more frequent ( 71%) than those on the same ones (Riechmann et al. 2000). Perhaps this trend is reflective of the relative abundance of ancient versus more recent gene duplication events.
Functional Divergence Among Conifer knox-I Genes
It has been shown early on that after gene duplication, one gene duplicate maintain the original function, whereas the other duplicate is free to accumulate amino acid changes as a result of functional redundancy and positive directional selection (Ohno 1970; Li 1983). Recent studies have indicated a more diverse array of potential outcomes in terms of fates of duplicated genes: pseudogeneization, functional redundancy, subfunctionalization in which each daughter gene adopts part of the functions of their parental gene, and neofunctionalization in which there is gain of novel function (for reviews, see Lawton-Rauh [2003] and Zhang [2003]). The differences between the last three categories are not so obvious in the absence of structural or biochemical data at the protein level. In our study, evidence for purifying selection was detected at the molecular level, with ratios of nonsynonymous (dN) to synonymous (dS) substitutions much smaller than 1 between sequences within each conifer knox-I gene as well as between conifer knox-I genes (dS values six times higher than dN values on average) (table 2). At their face value, these estimates indicate little evidence for directional selection. A similar trend towards low values (dS values five times higher than dN values) was observed between gene duplicates in Arabidopsis (Zhang, Vision, and Gaut 2002).
However, the ratio is based on rate averages since the split between two gene sequences. It does not take into account possible temporal heterogeneity in rates of evolution since the split between two gene duplicates, which could be indicative of periods of relaxed selection, followed by periods of more stringent selection. Such a rate shift is suggested by recent theoretical studies showing that duplicate genes experience two phases of evolution: a rather short period of relaxed selection and, thereafter, a period of purifying selection (Lynch and Conery 2003). When estimating rates of substitution per year for the diverse branches of the conifer tree (fig. 3), signatures of positive selection and functional divergence were evident for internal branches leading to the four conifer knox-I gene clades. R8s analyses showed that the calibrated rates of evolution per year were on average nine times larger before the split between Pinus and Picea than after, 140 MYA. Using the RRRT of Creevey and McInerney (2002), significant positive directional selection was also detected for several internal branches after the duplication of conifer knox-I genes. Such nonneutral divergence, early in the history of the conifer knox-I genes, is likely the result of a reduction in selective constraints on at least one duplicate, together with a commensurate increase in the amino acid replacement rate (Zhang 2003). At the opposite, short branches at the tip of the tree after the split between Pinus and Picea are indicative of a slow down in rates of evolution for the four conifer knox-I genes (fig. 3), which would imply more intense purifying selection in the recent past.
Relation Between Functional Divergence and Change in Map Position
The different chromosomal locations observed for the conifer knox-I paralogs were related to increased functional divergence. In our study, three unlinked locations were disclosed, corresponding respectively to (1) the genes KN1 and KN2, (2) the gene KN3, and (3) the gene KN4. These distinct locations parallel the significant results from RRRTs for positive selection and the large number of sites with significant rate shift from LRTs among the three groups. On the basis of model-based predictions, Lynch et al. (2001) suggested that the probability of change in map position per newly arisen gene duplicate increases with the strength of selection on neofunctional alleles. Similarly, the analysis of the genomic organization of multigene families in Arabidopsis suggested that the movement of gene duplicates to unlinked chromosomal regions may play an important role in functional divergence because of less exposure to sequence homogenization (Leister 2004).
The conifer paralogs KN1 and KN2 were positioned close to each other on the same linkage group, and there was little indication of significant positive directional selection between them, which is suggestive of more or less complete redundancy of function. At the same time, they presented a higher degree of sequence identity (average of 85% at the amino acid level), they had similar effective numbers of codon, %GC at third codon positions and transition/transversion ratios (table 3). It might be argued that this greater homogeneity is the result of a more recent gene duplication. However, the close physical relationship between KN1 and KN2 might have contributed to this greater homogeneity by exposing them to sequence homogenization effects. Indeed, the rates of substitution per site per year were lower for the branches leading to the KN1 or the KN2 clade (#11 and #14), as compared with those leading to the KN3 or KN4 clade (#7 and #3) and, more notably, to the branches representing the common ancestor of KN1 and KN2 (#15) and that of KN1, KN2, and KN3 (#16) (fig. 3).
Correlated Rates of Evolution Between Synonymous and Nonsynonymous Sites
A near homogenous was observed among conifer knox-I genes, indicating correlated rates of synonymous and nonsynonymous substitutions. Correlated temporal heterogeneity in rates of evolution was also observed between rates of change per year for first and second codon positions (mostly nonsynonymous) and rates of change per year for third codon positions (mostly synonymous): longer internal branches were noted before the divergence between Pinus, and Picea and shorter branches we noted after the divergence between Pinus and Picea (fig. 3). These trends are suggestive of factors affecting both synonymous and nonsynonymous sites in a similar way. A similar correlated trend between dS and dN was reported from the study of 10 nuclear conifer genes (Kusumi et al. 2002) and from the study of 242 duplicated genes from chromosomes 2 and 4 of Arabidopsis thaliana (Zhang, Vision, and Gaut 2002). This pattern could be indicative of background selection, which implies the removal of deleterious mutations as well as linked neutral or nearly neutral substitutions (Charlesworth, Morgan, and Charlesworth 1993). It would result in a decrease in the proportion of neutral mutations and a consequent correlated decrease in the rates of synonymous and nonsynonymous substitutions (Ohta 1995). Such a correlated selective constraint on both synonymous and nonsynonymous sites has also been observed in mammalian genes (Mouchiroud, Gautier, and Bernardi 1995).
This study demonstrates that knox-I genes did not evolve in a linear fashion in conifers. Initial periods of relaxed evolution have been followed by a period of conservatism, at least for the past 140 Myr. Taxon sampling in additional orders of the Gymnosperms would help delineate the beginning of this conservative period. A number of hypotheses at the functional level also need to be addressed by gene expression studies and investigations on biochemical properties and protein structure. For instance, it is anticipated that KN1 and KN2 would show more redundancy or overlapping in their physiological roles than with KN3 and KN4, which roles should be more distinct. These studies should assist in validating the putative functional importance of the changes observed at the molecular level and the extent of subfunctionalization or neofunctionalization among these genes.
Supplementary Material
The following information is available online: a list of primers used for the genomic amplification of spruce knox-I genes (table 1 of the Supplementary Material online) and pine knox-I genes (table 2 of the Supplementary Material online) and additional information about the mapping protocol of the four conifer knox-I genes in Picea mariana and Picea glauca (table 3 of the Supplementary Material online).
Acknowledgements
This work was carried out as part of the Arborea project lead by J. Mackay (Université Laval) and was supported by grants from Genome Québec and Genome Canada, the Natural Sciences and Engineering Research Council of Canada (NSERC, genomic program), the Canadian Biotechnology Strategy and the Canadian Forest Service. We thank R. K. Rutledge (Canadian Forest Service, Laurentian Forestry Centre, Québec) for the gift of PmKN4 from Picea mariana and for his advice and critical reading of a previous draft of this manuscript. We also thank F. Gagnon, S. Blais, and M. Ouellet for their assistance in the laboratory and their contribution to the RACE and RT-PCR experiments and two anonymous reviewers and B. Gaut for their insightful comments, which helped improve the manuscript.
References
Bharathan, G., B. J. Janssen, E. A. Kellogg, and N. Sinha. 1999. Phylogenetic relationships and evolution of the KNOTTED class of plant homeodomain proteins. Mol. Biol. Evol. 16:553–563.
Bouillé, M., and J. Bousquet. 2004. Trans-species shared polymorphisms at orthologous nuclear gene loci among distant species in the conifer Picea: implications for the maintenance of genetic diversity in trees. Am. J. Bot. (in press).
Byrne, M. E., J. Simorowski, and R. A. Martienssen. 2002. ASYMMETRIC LEAVES1 reveals knox gene redundancy in Arabidopsis. Development 129:1957–1965.
Champagne, C. E. M., and N. W. Ashton. 2001. Ancestry of KNOX genes revealed by bryophyte (Physcomytrella patens) homologs. New Phytol. 150:23–36.
Chang, S., J. Puryear, and J. Cairney. 1993. A simple and efficient method for isolating RNA from pine trees. Plant Mol. Biol. Rep. 11:113–116.
Charlesworth, B., M. T. Morgan, and D. Charlesworth. 1993. The effect of deleterious mutations on neutral molecular variation. Genetics 134:1289–1303.
Creevey, C. J., and J. O. McInerney. 2002. An algorithm for detecting directional and non-directional positive selection, neutrality and negative selection in protein coding DNA sequences. Gene 300:43–51.
Creevey, C. J., and J. O. McInerney. 2003. CRANN: detecting adaptive evolution in protein-coding DNA sequences. Bioinformatics 19:1726.
Florin, R. 1963. The distribution of conifer and taxad genera in time and space. Acta Horti. Bergiani 20:121–312.
Gaucher, E. A., X. Gu, M. M. Miyamoto, and S. A. Benner. 2002. Predicting functional divergence in protein evolution by site-specific rate shifts. Trends Biochem. Sci. 27:315–321.
Gu, X., and K. Vander Velden. 2002. DIVERGE: phylogeny-based analysis for functional-structural divergence of a protein family. Bioinformatics 18:500–501.
Hake, S., and N. Ori. 2002. Plant morphogenesis and KNOX genes. Nat. Genet. 31:121–122.
Hall, T. A. 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucl. Acids Symp. Ser. 41:95–98.
Hjortswang, H. I., A. Sundas-Larsson, G. Bharathan, P. V. Bozhkov, S. von Arnold, and T. Vahala. 2002. KNOTTED1-like homeobox genes of gymnosperm, Norway spruce, expressed during somatic embryogenesis. Plant Physiol. Biochem. 40:837–843.
Ingouff, M., I. Farbos, M. Wiweger, and S. von Arnold. 2003. The molecular characterization of PaHB2, a homeobox gene of the HD-GL2 family expressed during embryo development in Norway spruce. J. Exp. Bot. 54:1343–1350.
Ito, M., Y. Sato, and M. Matsuoka. 2002. Involvement of homeobox genes in early body plan of monocot. Int. Rev. Cytol. 218:1–35.
Jones, D. T., W. R. Taylor, and J. M. Thornton. 1992. The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci. 8:275–282.
Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21–132 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.
Kellis, M., B. W. Birren, and E. S. Lander. 2004. Proof and evolutionary analysis of ancient genome duplication in the yeast Saccharomyces cerevisiae. Nature 428:617–624.
Knudsen, B., M. M. Miyamoto, P. J. Laipis, and D. N. Silverman. 2003. Using evolutionary rates to investigate protein functional divergence and conservation: a case study of the carbonic anhydrases. Genetics 164:1261–1269.
Kumar, S., K. Tamura, I. B. Jakobsen, and M. Nei. 2001. MEGA2: molecular evolutionary genetics analysis software. Bioinformatics 17:1244–1245.
Kusumi, J., Y. Tsumura, H. Yoshimaru, and H. Tachida. 2002. Molecular evolution of nuclear genes in Cupressacea, a group of conifer trees. Mol. Biol. Evol. 19:736–747.
Lawton-Rauh, A. 2003. Evolutionary dynamics of duplicated genes in plants. Mol. Phylogenet. Evol. 29:396–409.
Leister, D. 2004. Tandem and segmental gene duplication and recombination in the evolution of plant disease resistance genes. Trends Genet. 20:116–122.
Li, W. H. 1983. Evolution of duplicated genes. Pp. 14–37 in M. Nei and R. K. Koehin, eds. Evolution of genes and proteins. Sinauer, Sunderland, Mass.
Lynch, M., and J. S. Conery. 2003. The evolutionary demography of duplicate genes. J. Struct. Funct. Genomics 3:35–44.
Lynch, M., M. O'Hely, B. Walsh, and A. Force. 2001. The probability of preservation of a newly arisen gene duplicate. Genetics 159:1789–1804.
Miller, C. N. 1988. The origin of modern conifer families. Pp. 448–486 in C. B. Beck, ed. Origin and evolution of Gymnosperms. Columbia University Press, New York.
Mouchiroud, D., C. Gautier, and G. Bernardi. 1995. Frequencies of synonymous substitutions in mammals are gene-specific and correlated with frequencies of nonsynonymous substitutions. J. Mol. Evol. 40:107–113.
Nam, J., C. W. dePamphilis, H. Ma, and M. Nei. 2003. Antiquity and evolution of the MADS-box gene family controlling flower development in plants. Mol. Biol. Evol. 20:1435–1447.
Nei, M., and T. Gojobori. 1986. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol. Biol. Evol. 3:418–426.
Nei, M., I. B. Rogozin, and H. Piontkivska. 2000. Purifying selection and birth-and-death evolution in the ubiquitin gene family. Proc. Natl. Acad. Sci. USA 97:10866–10871.
Ohno, S. 1970. Evolution by gene duplication. Allen and Unwin, London.
Ohta, T. 1995. Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral theory. J. Mol. Evol. 40:56–63.
Pelgas, B., N. Isabel, and J. Bousquet. 2004. Efficient screening for expressed sequence tag polymorphisms (ESTP) by DNA pool sequencing and denaturing gradient gel electrophoresis in spruces. Mol. Breed. 13:263–279.
Pham, T., and N. Sinha. 2003. Role of KNOX genes in shoot development of Welwitschia mirabilis. Int. J. Plant Sci. 164:333–343.
Prager, E. M., D. P. Fowler, and A. C. Wilson. 1976. Rates of evolution in conifers (pinaceae). Evolution 30:637–649.
Reiser, L., P. Sanchez-Baracaldo, and S. Hake. 2000. Knots in the family tree: evolutionary relationships and functions of knox homeobox genes. Plant Mol. Biol. 42:151–166.
Riechmann, J. L., J. Heard, G. Martin et al. (17 co-authors). 2000. Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science 290:2105–2110.
Rozas, J., J. C. Sanchez-DelBarrio, X. Messeguer, and R. Rozas. 2003. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19:2496–2497.
Sanderson, M. J. 2002. Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol. Biol. Evol. 19:101–109.
———. 2003. r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. Bioinformatics 19:301–302.
Sato, Y., Y. Fukuda, and H. Y. Hirano. 2001. Mutations that cause amino acid substitutions at the invariant positions in homeodomain of OSH3.KNOX protein suggest artificial selection during rice domestication. Genes Genet. Syst. 76:381–392.
Savard, L., P. Li, S. H. Strauss, M. W. Chase, M. Michaud, and J. Bousquet. 1994. Chloroplast and nuclear gene sequences indicate late Pennsylvanian time for the last common ancestor of extant seed plants. Proc. Natl. Acad. Sci. USA 91:5163–5167.
Schena, M., and R. W. Davis. 1994. Structure of homeobox-leucine zipper genes suggests a model for the evolution of gene families. Proc. Natl. Acad. Sci. USA 91:8393–8397.
Semiarti, E., Y. Ueno, H. Tsukaya, H. Iwakawa, C. Machida, and Y. Machida. 2001. The ASYMMETRIC LEAVES2 gene of Arabidopsis thaliana regulates formation of a symmetric lamina, establishment of venation and repression of meristem-related homeobox genes in leaves. Development 128:1771–1783.
Sundas-Larsson, A., M. Svenson, H. Liao, and P. Engstrom. 1998. A homeobox gene with potential developmental control function in the meristem of the conifer Picea abies. Proc. Natl. Acad. Sci. USA 95:15118–15122.
Swofford, D. L. 2002. PAUP*: phylogeneric analysis using parsimony (* and other methods). Sinauer Associates, Sunderland, Mass.
Takada, S., and M. Tasaka. 2002. Embryonic shoot apical meristem formation in higher plants. J. Plant Res. 115:411–417.
Thompson, J. D., T. J. Gibson, F. Plewniak, F. Jeanmougin, and D. G. Higgins. 1997. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25:4876–4882.
Tioni, M. F., D. H. Gonzalez, and R. L. Chan. 2003. Knotted1-like genes are strongly expressed in differentiated cell types in sunflower. J. Exp. Bot. 54:681–690.
Ujino-Ihara, T., K. Yoshimura, Y. Ugawa, H. Yoshimaru, K. Nagasaka, and Y. Tsumura. 2000. Expression analysis of ESTs derived from the inner bark of Cryptomeria japonica. Plant Mol. Biol. 43:451–457.
Van Ooijen, J. W., and R. E. Voorrips. 2001. JoinMap: Software for the calculation of genetic linkage maps. Version 3.0. Plant Research International, Wageningen, The Netherlands.
Veit, B. 2004. Determination of cell fate in apical meristems. Curr. Opin. Plant Biol. 7:57–64.
Vollbrecht, E., B. Veit, N. Sinha, and S. Hake. 1991. The developmental gene Knotted-1 is a member of a maize homeobox gene family. Nature 350:241–243.
Wang, Y., and X. Gu. 2001. Functional divergence in the caspase gene family and altered functional constraints: statistical analysis and prediction. Genetics 158:1311–1320.
Ware, D., P. Jaiswal, J. Ni et al. (12 co-authors). 2002. Gramene: a resource for comparative grass genomics. Nucleic Acids Res. 30:103–105.
Wei, X. X., and X. Q. Wang. 2004. Evolution of 4-coumarate:coenzyme A ligase (4CL) gene and divergence of Larix (Pinaceae). Mol. Phylogenet. Evol. 31:542–553.
Wright, F. 1990. The ‘effective number of codons’ used in a gene. Gene 87:23–29.
Xue, B., P. J. Charest, Y. Devantier, and R. G. Rutledge. 2003. Characterization of a MYBR2R3 gene from black spruce (Picea mariana) that shares functional conservation with maize C1. Mol. Genet. Genomics 270:78–86.
Yang, Z. 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. 13:555–556.
Yang, Z., and R. Nielsen. 2000. Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17:32–43.
Zhang, J. 2003. Evolution by gene duplication: an update. Trends Ecol. Evol. 18:292–298.
Zhang, L., T. J. Vision, and B. S. Gaut. 2002. Patterns of nucleotide substitution among simultaneously duplicated gene pairs in Arabidopsis thaliana. Mol. Biol. Evol. 19:1464–1473.(Carine Guillet-Claude*, N)