BACKGROUND
organisms adapt to different habitats through natural selection, which favors the fixation of alleles that increase the fitness of the individual that bears them.
however, it is quite difficult to identify the locus/loci targeted by selection.
one reason is that the number of loci involved in a particular adaptation and their phenotypic effects vary depending on the genetic architecture underlying the adaptive trait  <cit> .
another reason is that often the sequence and/or information about the gene involved in the adaptive response is unavailable.
two main approaches are used to identify genetic signatures of adaptive evolution and link them to phenotypic traits.
the first is the candidate gene approach, whereby signatures of positive selection  are identified using a population genetics framework in those genes previously known to be involved in the phenotypic trait of interest .
the advantage of this approach is that the consequences of molecular variation on phenotypes can be inferred and their adaptive significance evaluated  <cit> .
on the other hand, this approach assumes a comprehensive knowledge of the function of such genes and may neglect other genes relevant to the phenotypic trait.
the second is the genome-wide approach, where the pattern of molecular evolution of hundreds to thousands of genes scattered across the genome are analyzed simultaneously .
as high throughput technologies are becoming cheaper and more accessible, this approach is increasingly gaining attention from the evolutionary biology research community.
an important advantage of this method is that it does not rely on prior information about gene sequence, and therefore it is well suited for studying non-model species.
the possibility of performing genome scale studies in non-model organisms is indeed a powerful approach to address the genetic basis of specific adaptations, which can only be obtained by choosing the appropriate organism and ecological context.
however, the typically poor knowledge about gene function in non-model organisms often prevents a comprehensive understanding of the adaptive significance of eventual signatures of positive selection.
a useful compromise is to use the genome-wide approach in species closely related to well-studied model organisms, so that gene function can be inferred by comparative analyses.
in this way it is possible to exploit what is known about the ecology and the life history of the species, and thus the approach is particularly suited to identifying genes involved in species-specific and habitat-specific adaptations .
environmental variation is a key factor driving adaptive evolution and determining the ecological niche of a species.
plants and other sessile organisms are particularly affected by circadian and seasonal oscillations in abiotic factors.
for example, sudden drops in temperature, high levels of solar irradiation and limited access to water are common sources of environmental stress, especially at high altitude in mountainous regions  <cit> .
these stressors affect the evolution of species living in these habitats, and their capacity to adapt to these stressors ultimately determines their distribution  <cit> .
to cope with drops in temperature, plants have developed a series of physiological adaptations that rely on the up- and down-regulation of cold responsive genes triggered by cold exposure .
similarly, cold temperatures and high irradiation are not favorable to efficient photosynthesis, and consequently, a suite of photoprotective strategies are required for survival and reproduction at high altitudes .
a fairly detailed understanding of the relevant regulatory pathways and gene function in the model species arabidopsis thaliana now exists; however, little is known about their adaptive role, particularly in relation to the diverse habitats present along an altitudinal gradient.
in order to assess the adaptive role of cold responsive genes, as well of the genes involved in photosynthesis, we compared patterns of gene evolution in congeneric species living at different altitudes.
we performed molecular evolution analyses in two closely related brassicaceae species adapted to non-overlapping altitudinal ranges: cardamine resedifolia, a perennial species usually growing under conditions of high irradiation and severe temperature oscillations between  <dig>  and  <dig>  meters above sea level; and c. impatiens, an annual nemoral species normally growing between  <dig> and  <dig>  meters above sea level  <cit> .
the distinct habitats associated with high and low altitudes make c. resedifolia and c. impatiens ideal species for studying the adaptive significance of the genes involved in altitude-related stress responses.
it should also be noted that the two species may differ in their outcrossing rates, even though conclusive evidence is still lacking: c. impatiens is mainly selfing  <cit> , although some populations have been reported as partially outcrossing  <cit> ; c. resedifolia is also a predominantly selfing species, but the precise outcrossing rate is unknown  <cit> .
the two cardamine species described above are closely related to a. thaliana  <cit> , making it possible to apply the extensive knowledge about the gene functions of this model organism to our study.
gene sequences of both c. resedifolia and c. impatiens were obtained by high-throughput sequencing technology and subsequently identified and partitioned into functional classes based on gene similarity to the a. thaliana orthologues.
then, we used complementary approaches, such as analyses of non-synonymous and synonymous substitutions and of their ratio, of levels of selective pressure, codon usage, and gene expression, to quantify the difference in adaptive evolution between functional classes and between species.
this allowed us to infer the adaptive significance of genes involved in adaptation to cold stress and photosynthesis at high altitude and put them into an ecological context.
RESULTS
we obtained sequences representing the c. resedifolia and c. impatiens transcriptomes using high-throughput sequencing.
genes and their putative function were identified by comparison of the sequences obtained here to the available annotated a. thaliana gene sequences.
to minimize the chance of mistaking a paralogue for an orthologue, we considered as triplets of putative orthologues only those consisting of reciprocal best hits, i.e.
those where the three sequences were consistently found as best hit matches of one another.
eventually, we obtained  <dig>  triplets of  nuclear genes, with a mean  length for the a. thaliana orthologues of  <dig>  ±  <dig>  base pairs , corresponding to a mean coverage of  <dig>  ±  <dig> % of the full-length gene sequence.
in c. impatiens, the mean sequence length was  <dig>  ±  <dig>  bp, while in c. resedifolia the mean sequence length was  <dig>  ±  <dig>  bp.
we partitioned these genes according to their putative function, and then focused our analyses on those functional classes that are associated with the adaptive response to high altitude .
in particular, according to the annotation in a. thaliana,  <dig> genes were involved in cold acclimation ,  <dig> were involved in uv-b and high irradiation response ,  <dig> were involved in general stress responses .
to these three classes we added a manually compiled class that included  <dig> genes functionally characterized as cold responsive genes .
rate of molecular substitution in cardamine
the analysis of the rates of nucleotide substitution revealed faster molecular evolution along the c. resedifolia lineage than along the c. impatiens lineage .
because the rates of sequence evolution were highly correlated with the length of the a. thaliana orthologous gene , we corrected these measures accordingly by using the residuals of such correlations .
while the rate of synonymous substitution, ds, was similar between lineages , the rate of non-synonymous substitution, dn, was larger in the c. resedifolia lineage than in the c. impatiens lineage .
conversely, the ratio ω = dn/ds, was significantly larger in the c. impatiens than in the c. resedifolia lineage .
c. impatiens
c. resedifolia
n
a
n
a
p
c
a number of genes.
b mean and standard error.
c wilcoxon test probability after correcting for gene length
for each functional class, both the dn and the ds rates, as well their ratio ω, were not significantly different between the two lineages after correcting for multiple testing .
selective constraints varied however across functional classes within the same lineage: dn was significantly lower in pgo than in other genes in the c. resedifolia lineage , and was significantly higher in crg than in other genes in both cardamine lineages .
the statistically significant differences were not due to random effects resulting from the small sample size of the functional classes relative to the full dataset.
in fact, a bootstrap approach revealed that at most  <dig> % of the randomized datasets produced wilcoxon tests p values lower than those observed .
since there was a significant correlation between dn and the length of the a. thaliana orthologue , we also re-analyzed the comparisons using the residuals of the correlation.
the within-lineage differences were still statistically significant, thus excluding the possibility that such differences were a result of length heterogeneity across our gene sets .
similarly, the mean values of the residuals for dn and ds still did not differ between lineages, with the exception of the pgo genes, where the dn residuals were lower in c. resedifolia than in c. impatiens .
the analysis of the levels of selective pressure, ω, provided compelling evidence that different selective forces are dominating the evolution of the four functional categories in the two cardamine lineages .
again, as a result of the correlation between substitution rates and gene length, we report the results of the residuals of such correlations.
in the c. resedifolia lineage, the mean value of ω was significantly lower in pgo than in other genes , consistent with intense purifying selection.
this is in sharp contrast with what observed in crg, where the mean value of ω was significantly higher than in other genes , indicating either relaxed or positive selection.
instead, in the c. impatiens lineage, the levels of selective pressure in genes of the four functional classes were similar to those at the genome-wide level .
in c. resedifolia, selective pressures were different between genes involved in cold response classified as cgo and those classified as crg.
the two datasets have  <dig> genes in common , whose mean ω did not differ from that of the remaining genes .
however, the genes unique to cgo and crg had a mean ω in line with what was observed in the complete datasets.
that is, genes unique to the cgo category had a mean ω significantly lower than the remaining genes , confirming strong purifying selection in cgo.
conversely, genes exclusive to the crg category had a mean ω significantly higher than the other genes , indicating either relaxed or positive selection in these genes.
the rate of evolution and the occurrence of positive selection were further investigated by analyzing the pattern of molecular substitution of each single gene.
the most sensitive test for this purpose is one that employs branch-site models, which aim to detect positive selection affecting a few sites along particular lineages.
when c. impatiens was used as the foreground lineage  the test identified one outlier with fdr <  <dig> .
this number increased substantially when the c. resedifolia lineage was used instead , with seven genes showing evidence for positive selection along this lineage .
two genes among these eight outliers were identified as being involved in general stress responses .
a likelihood ratio tests  in which the alternative model explained the pattern of codon substitution better than the null model, and with a false discovery rate lower than  <dig> .
patterns of codon substitution were tested using branch models , site models , and branch-site models .
b tests compared models m <dig> and m0'.
s tests compared model m2a to m1a  and model m <dig> to m <dig> .
bs tests were based on the branch-site model a, test  <dig>  subscripts in the b and bs tests indicate whether the c. impatiens  or the c. resedifolia  lineages were used as foreground branches .
b gene of the sgo functional class.
to increase the chance of discovering interesting candidate genes, we also employed additional maximum likelihood ratio tests based on branch or site codon substitution models  <cit> .
these tests identified several genes whose pattern was better explained by the occurrence of non-neutral evolution .
branch models identified three genes  with a pattern of nucleotide substitution that fit a model allowing two ω rates, one for the focal c. impatiens branch and a second for the rest of the phylogenetic tree .
a single gene was detected if the c. resedifolia branch was used as the focal branch instead .
using site models, the fit of the data under a neutral model was compared to that under a model of positive selection.
in the first comparison, the neutral model admitted one class of sites with  <dig> < ω <  <dig>  while the selection model allowed the presence of a second class of sites with ω >  <dig> .
the second, more powerful , approach compared more realistic models where  neutrally evolving sites were partitioned in ten classes with ω values drawn from a beta distribution .
the two likelihood ratio tests consistently identified four genes  showing evidence for the action of positive selection .
as expected, the tests based on either the branch, site or branch-site models resulted in the detection of sets of outliers with little overlap.
among the three genes detected by more than one test, two  were detected by both the b and the bci tests, and one  was detected by the s <dig>  s <dig> and bsci tests.
to exclude the possibility that substitution rate patterns were the result of orthology mis-assignment, effectively inflating interspecific divergence, we compared homologous sequences from a. thaliana, cardamine and other plant species within a phylogenetic framework.
specifically, we searched for putative paralogues of the gene triplets, for each of the  <dig> candidate genes .
in six such cases we verified the presence of putative paralogues using blast searches .
the a. thaliana and the two putative cardamine orthologues identified by the reciprocal best-hit approach always clustered together in the tree, supporting their orthology and the results of the likelihood ratio tests.
expression breadth, expression level and rate of molecular substitution
the breadth of expression of a gene, namely whether it has tissue-specific expression or is expressed broadly over more tissues, has been shown to be correlated with several patterns of molecular evolution .
to evaluate such effects in cardamine, we correlated the levels of selective pressure, ω, to the breadth of expression measured across tissues, across developmental time points, and along stress treatments.
for these analyses we used the expression patterns of a. thaliana orthologues as proxies for expression in cardamine , assuming that orthologous genes will maintain similar expression patterns  <cit> .
these analyses found strong evidence that, in both cardamine lineages, genes evolved at different rates and under different selective regimes depending on their breadth of expression .
the ratio ω decreased significantly with the increase in the temporal breadth of expression, i.e.
the number of flower and leaf developmental stages at which a gene was expressed .
the same statistically significant trend held when correlating ω with the spatial breadth of expression , i.e.
the number of tissues/organs where the gene is expressed and the associated organ specificity index τ  <cit> .
because gene length was significantly correlated with both ω  and expression breadth , we repeated our analyses using the residuals of such correlations .
the partial correlations coefficients were still statistically highly significant  and indicated that the correlations between ω and breadth of expression were independent from gene length.
thus, a broad breadth of expression imposes selective constraints limiting the possible action of positive selection.
conversely, genes expressed at one developmental time point or in specific organs evolve faster, as a result of either relaxed or positive selection.
c. impatiens
c. resedifolia
* spearman's p <  <dig> , ** p <  <dig> , *** p <  <dig> , **** p <  <dig> , ns = non significant.
^ partial correlations are significant after correcting for the length of the a. thaliana orthologue.
§ significant after holm-bonferroni multiple test correction
a breadth of expression can be either spatial  or temporal .
note that the organ-specificity index τ is inversely correlated with the number  of organs at which a gene is expressed.
when ω was compared to the duration  of the stress response, we obtained mixed results that depended on the stress type .
in general, the ratio ω was positively correlated to the persistence of gene expression during stress responses .
after a holm-bonferroni correction for multiple testing, ω was significantly correlated to the duration of the responses to salt and uv-b in the c. impatiens lineage ; and to the duration of the responses to salt, osmotic and cold stress in the c. resedifolia lineage .
in addition to the breadth of expression, the level of gene expression  was also an important determinant of ω .
for instance, there was a statistically highly significant negative correlation between the mean gene expression and ω in both cardamine lineages .
this correlation also remained highly significant when the mean expression level was controlled for the effects of gene length .
thus, we conclude that selective constraints may limit the evolution of proteins encoded by highly expressed genes in cardamine.
c. impatiens
c. resedifolia
n
b
p
d
p
d
mean expression
maximum expression
* spearman correlation p <  <dig> , ** p <  <dig> , *** p <  <dig> , **** p <  <dig> , ns = non significant.
a gene functional class.
b number of genes.
c mean/maximum expression was estimated for all genes and genes included in each functional class and tested  against mean for genes not in that functional class.
d spearman's correlation between level of gene expression and ω estimated along the two cardamine lineages.
the relationship between gene expression and rate of evolution indicates that the different selective pressures we observed in genes involved in stress response may be correlated to their expression level.
for instance, genes of the cgo and pgo functional classes had a calculated expression level higher than other genes , and also had lower ω along the c. resedifolia lineage than the genome-wide mean .
to investigate the association between gene expression and ω, we then calculated the residuals of their correlation and compared them across functional classes .
the corrected ω were no longer statistically significant for pgo and cgo genes, suggesting that the levels of selective pressure in these genes were highly associated with their high expression levels.
interestingly, removing the effect of expression level produced a statistically significant difference in ω between sgo genes and the other genes .
such differences suggests heterogeneity in the association between rates of molecular evolution and expression levels in sgo genes, likely due to the diverse nature of the stress responses to which these genes are responding.
results did not change if we analyzed the correlation between maximum, rather than mean, expression level and ω , further indicating a strong association between the expression pattern of a gene and its rate of molecular evolution.
codon usage in cardamine genes
synonymous codon usage can be under weak selection and lead to non-random use of the codons coding for the same aminoacid.
due to hill-robertson interference, such bias is expected to correlate with the occurrence of positive selection, so that it will be lower in genes that experience  episodes of positive selection  <cit> .
consistent with this prediction, in both cardamine lineages there was a statistically significant negative correlation between ω and codon bias, measured as fop .
this correlation also held when controlling for the effects of other determinants of codon usage, such as expression levels, gene length and gc content .
for instance, in both cardamine lineages, fop was significantly correlated with expression level , gene length , and gc content at the third codon position .
however, partial correlation coefficients  between codon usage and ω remained highly significant .
we then analyzed codon usage for the genes in each of the four functional classes.
genes involved in stress responses generally had larger codon bias compared to other genes .
the difference was statistically significant for cgo and sgo genes in both cardamine lineages , while it was not significantly different for crg and pgo genes .
however, such differences were no longer significant when controlling for gene expression level , suggesting that the expression pattern influenced codon usage more than the distinct selective pressures to which they are subject.
thus, codon usage bias does not support the frequent action of positive selection in genes of the functional classes under consideration.
discussion
our analyses of protein-coding sequences highlighted lineage- and gene-specific evolutionary patterns in two species of cardamine characterized by distinct habitat preferences, life-histories and, possibly, breeding systems.
c. resedifolia and c. impatiens are not sister species  <cit> ; hence, the substitution rates estimated in the present study apply to the lineages leading to the two species after the split from their common ancestor, rather than to these specific taxa as such.
however, c. resedifolia, c. impatiens, and the clades they belong to, are characterized by divergent phenotypic and life-history traits that have profoundly shaped their evolutionary histories; for instance, the clade comprising c. resedifolia  is basal to the cardamine radiation and includes only perennial species living in high altitude alpine habitats with well-developed petals indicating the existence of a mixed mating system.
in contrast, c. impatiens is, to our knowledge, the only annual selfing species of its group , which includes species that are mainly found at moderate elevations.
c. impatiens is further characterized by having very reduced or no petals, in agreement with a predominantly selfing mating system  <cit> .
thus, the observed patterns of molecular evolution reflect the differences in life history traits and habitats between the two lineages.
the effects of such differences on the rate of molecular evolution are quite complex.
for instance, the rate of non-synonymous substitution was significantly higher for the c. resedifolia lineage than for that of c. impatiens.
previous studies have reported a higher substitution rate in annual plants compared to perennial plants  <cit>  .
however, our observations are in contrast with this prediction, since the annual c. impatiens had lower substitution rate than the perennial c. resedifolia.
moreover, since the rate of synonymous substitution was similar between lineages, there is no support for the hypothesis that mutation rate is higher along the c. resedifolia lineage than along that of c. impatiens.
the interspecific differences in non-synonymous substitution rates could also be the result of differences in effective population size  <cit> .
for instance, low effective population size can affect the substitution rate by allowing slightly deleterious mutations to escape purifying selection and reach fixation .
alternatively, high effective population size can increase the efficiency of positive selection by facilitating the fixation of favorable alleles in genes undergoing adaptive evolution  <cit> .
thus, a high dispersion in the distribution of dn across the genome is expected to be associated with a large effective population size, especially when there is recombination .
the analysis of non-synonymous substitution fixation rate revealed that the variance in dn was significantly larger in c. resedifolia than in c. impatiens .
this difference in protein evolution heterogeneity is thus consistent with more efficient positive selection and higher effective population size in c. resedifolia.
both of these factors, in turn, would explain why dn is higher in this lineage than in c. impatiens.
preliminary analyses of recombination levels and polymorphism suggest indeed that selfing rates are lower in c. resedifolia, and the effective population size is slightly larger than in c. impatiens .
intriguingly, mean ω was larger in c. impatiens than in c. resedifolia.
the fact that in c. impatiens ω was considerably lower than one is compatible with either moderately positive or relaxed selection.
three lines of evidence suggest a reduced efficiency of purifying selection in c. impatiens, rather than pervasive positive selection at a genome-wide level.
firstly, the analyses on dn suggest that positive selection was more efficient along the c. resedifolia lineage than the c. impatiens lineage.
secondly, only one gene  exhibited strong evidence for site-specific positive selection along the c. impatiens branch.
by comparison, seven genes exhibited strong evidence for site-specific positive selection along the c. resedifolia branch.
finally, the analysis of codon usage failed to detect differences in the efficiency of positive selection between the two lineages.
however, the contrast between rates of non-synonymous substitution and of ω warrants some caution in drawing firm conclusions.
the most likely explanation for the opposite patterns observed in dn and dn/ds lies on the high heterogeneity in substitution rate among genes.
for instance the correlations and the coefficients of determination in the regressions between ds and ω  and between dn and ω , suggest that substitution rates are relatively weak predictors of the actual selective pressure experienced by a gene.
the statistically significant correlation between synonymous substitution rate and ω is rather unusual: previous studies established that the correlation is caused by a combination of  adjacent substitutions  <cit>  and  of a bias on substitution rate estimation due to either low divergence and/or statistical methods  <cit> .
appropriate analyses will be necessary to unravel the effects of such factors in our system.
the complex dynamics of substitution rates are also evident in the statistically significant correlation to gene length.
the influence of gene length on sequence evolution had been previously reported for populus tremula  <cit> , stressing the importance of accounting for all diverse determinants of levels of gene and protein evolution.
additional studies, including analyses on intraspecific polymorphism, will be certainly necessary to disentangle the neutral and selective forces that have shaped such patterns  <cit> .
the molecular evolution analyses in stress-related genes also revealed important lineage-specific patterns that may be associated with the distinct life-history traits and habitats of c. resedifolia and c. impatiens.
for instance, similar selective regimes affected the evolution of genes involved in stress response in the c. impatiens lineage.
conversely, there was a correlation between the type of stress response and the rate of molecular evolution along the c. resedifolia lineage; for example, genes involved in photosynthesis evolved slower than other genes, consistent with selective constraints that limited the accumulation of non-synonymous substitution.
this makes sense given their involvement in a process that is particularly relevant in the high altitude habitat of c. resedifolia.
strangely, in c. resedifolia selection acted in opposite directions in the two functional classes associated with cold responses.
that is, compared to the genome-wide levels of selective pressure, genes that were identified as involved in cold response based on functional studies  displayed significantly lower levels of selective pressure, while cold responsive  genes were under more selective constraints.
since there was little overlap between the two classes , the discrepancy in the levels of selective pressure between cgo and crg genes may have captured different selective pressures acting along the cold response pathway.
however, the difference may also stem from the approaches used to define the two cold-response functional classes.
specifically, the gene ontology annotations of the cgo genes were taken from the tair database  <cit> , which is a less accurate source of functional information than the direct assays used to define crg genes.
for instance, the tair annotations may be based solely on sequence or structural similarity, and a gene may have been annotated as cold responsive because it is indirectly up- or down-regulated in plants exposed to low temperatures .
however, without a comprehensive knowledge on the response pathway it is difficult to evaluate their role in cold response in cardamine, especially considering that their expression patterns may differ from those of a. thaliana.
the levels of selective pressure were also associated with the duration and pattern of gene expression.
for instance, in c. impatiens and c. resedifolia, ω was lower for genes that were expressed only briefly during specific stress responses than for genes that were up-regulated over a longer period of time .
most importantly, the correlation existed only for specific stress responses in each species, suggesting habitat-specific selective pressures.
in particular, in c. resedifolia, ω was significantly affected by the duration of the responses to osmotic, salt and cold stress.
the responses to these three stresses involve partially overlapping pathways  <cit> , as cold stress causes membrane leakage and, as a consequence, the activation of the physiological responses also observed in high salt and osmotic stresses  <cit> .
these results support the adaptive relevance of these genes in response to the severe temperature changes that this species experiences in alpine habitats.
on the other hand, in c. impatiens ω was significantly correlated with the extent of up-regulation in those genes involved in uv-b stress response.
c. impatiens is a nemoral species that is likely to be particularly sensitive to the effect of exposure to uv-b light.
therefore, it is reasonable to expect that this species has evolved response mechanisms based on gene up-regulation to cope with discontinuous uv stress.
it will be extremely interesting to experimentally measure the expression of the cardamine genes in the functional classes considered, as species-specific adaptive variations in expression levels cannot be excluded.
while unequal selective pressures across the stress-related functional classes were well-documented here, we detected positive selection in only a few single genes.
in plants, previous studies have found evidence of genome-wide positive selection in some species , but in most cases there was little indication of widespread adaptive evolution .
other studies have identified positive selection of some genes involved in stress response .
the most likely explanation for the low rate of adaptive evolution is that most plants have low effective population size  <cit> , which ultimately diminishes the efficiency of positive selection.
estimation of the effective population size of c. resedifolia and c. impatiens will be useful to verify whether this parameter can explain the scarce evidence for positive selection in our dataset.
four methodological issues may have further resulted in an underestimation of the genes targets of positive selection.
the first is that the power of codon substitution models is fairly conservative compared to that of models incorporating polymorphism data .
the second is that our reciprocal best-hit approach was prone to miss genes with high sequence divergence caused by adaptive evolution, as it was intended to avoid an overestimation of divergence by comparing paralogues.
in particular, this approach may have overlooked duplicated genes, which are common in a. thaliana  <cit>  and in other plants , and which can undergo sub- or neo-functionalization driven by positive selection  <cit> .
for instance, our dataset included only ~11% of all a. thaliana genes, and we analyzed between 17% and 48% of the genes annotated as involved in the stress responses considered in the present study.
an examination of our dataset revealed that many of the well-characterized transcription factors involved in cold response  were absent.
this indicates either that we are missing many of the genes upstream of the stress response pathways, or that in cardamine these transcription factors may be differently regulated than in arabidopsis thaliana or do not participate to the cold response.
however, it is unlikely that this pathway is missing, since cold response is well-conserved across distantly related taxa .
a third issue is related to the use of partial genes, which may have reduced the power of the likelihood ratio tests as a result of an insufficient number of informative sites.
finally, genes expressed at a low level, including the aforementioned transcription factors, may be missing from our dataset because of insufficient coverage or normalization.
because genes expressed at a low level are those evolving faster, this bias in gene representation could contribute to the relatively low number of rapidly evolving genes identified in this study.
deeper sequencing efforts could undoubtedly improve the situation by increasing both coverage and average transcriptome length.
one c. impatiens gene and seven c. resedifolia genes showed signatures of positive selection under the branch-sites model of codon substitution.
the c. impatiens gene orthologue of at4g <dig> has not been functionally characterized in a. thaliana, although it is known that the protein encoded by this genes displays homology to the hyaluronan/mrna binding protein family, a group of proteins binding both specific rna and a high-molecular-mass polysaccharide extremely abundant in the connective tissue and extracellular matrix of animals  <cit> .
specific studies will be necessary to understand its role in adaptive processes in the c. impatiens lineage.
as for the c. resedifolia candidate genes, no functional information is available for the orthologue of at1g <dig>  which codes for a protein of unknown function with homology to tolb, a protein involved in outer membrane stability and uptake of biomolecules in e. coli  <cit> .
instead, at3g <dig> is known to code for a protein putatively involved in metal ion binding, and is similar to proteins of the heavy metal transport/detoxification superfamily.
heavy metal hyperaccumulation has been documented in several brassicaceae species  <cit>  and c. resedifolia has been reported to accumulate large quantities of nickel from the nickel-rich debris of the glacial till, where this species usually grows  <cit> .
therefore, it is possible that the signature of positive selection identified in this orthologue is related to high contents of heavy metals in soil experienced by c. resedifolia.
a third gene, orthologue of at1g <dig>  is a valine-trna synthetase.
given the multiple metabolic pathways in which valine is involved , it would be highly speculative to associate this gene to adaptive processes in the c. resedifolia lineage.
based on functional evidence from arabidopsis and other plant species, two other genes identified as putative targets of positive selection in c. resedifolia may play major roles in the response against bacterial pathogens and insect herbivory.
the first gene, at5g <dig>  is part of the jasmonate-zim-domain protein family  <cit> , whose members are involved in response to wounding and herbivory  <cit> .
the other candidate is the orthologue of at1g <dig>  which codes for the arabidopsis epithiospecifier protein .
esp catalyzes the formation of simpler nitriles and epithionitriles from glucosinolates, thus modulating the release of isothiocyanate, a metabolite involved in herbivory and pathogen defense in brassicaceae  <cit> .
in two closely related boechera species , the level of glucosinolates is negatively correlated with elevation preferences and growth rates, but positively correlated with drought tolerance  <cit> .
however, the response to herbivory as a function of elevation is not uniform across plants, with some species experiencing more  and other less  damage with increase in elevation.
interestingly, a positive correlation with levels of insect herbivory was found for c. cordifolia exposed to full sunlight, possibly resulting from moderate water stress associated with a different insect guild compared to shadowy environments  <cit> .
these observations provide the framework to experimentally test whether the signature of positive selection identified in at5g <dig> and at1g <dig> orthologues might be related to the higher exposition to sun, lower water availability or slower growth rate characterizing c. resedifolia as compared to c. impatiens.
interestingly, the orthologues of at1g <dig> and at3g <dig>  may also be related to the light regimes characterizing c. resedifolia and c. impatiens habitats.
at1g <dig> codes for a cytosolic ascorbate peroxidase  that scavenges hydrogen peroxide in plant cells  <cit> , thus reducing the accumulation of reactive oxygen species  that cause cellular damage through protein oxidation  <cit> .
the product of at3g <dig> is a transcriptional activator of the growth regulating factor family.
genes from this family have been demonstrated to be involved in the regulation of cell expansion and division in leaves, cotyledons and petals  <cit> .
strikingly, the product of at3g <dig> is one of the proteins oxidized in apx <dig> mutant plants in response to moderate light stress  <cit> , indicating that it could also be part of the signaling cascade activated by ros stress.
the site codon substitution models also identified putative targets of positive selection that deserve further characterization.
in particular, the orthologue of at2g <dig> codes for a ribosomal protein  that is involved in response to salt and genotoxic stress in a. thaliana  <cit> .
this gene is part of the ribosomal protein s <dig> family, which includes three paralogues  with similar sequences and function.
gene families can undergo relaxed selection immediately following duplication  <cit> ; however, this does not seem to be the case for at2g <dig>  as the phylogenetic tree clearly indicates that the duplication occurred before the split between the arabidopsis and cardamine lineages .
additional studies will be necessary to corroborate these findings and link the evolutionary pattern of each gene to its phenotypic effects.
in particular, the use of intraspecific variation and functional analyses in the model species a. thaliana will be crucial to ascertaining whether positive selection or relaxed selection accelerated the evolution of these genes and their relevance in adaptive processes in cardamine.
CONCLUSIONS
overall, our results highlight the importance of employing complementary approaches to studying the genetic bases of adaptation in non-model species.
our use of comparative genomics on congeneric species identified evolutionary patterns that aid the understanding of the extrinsic and intrinsic factors driving plant adaptation.
in the case of cardamine, intrinsic factors  most likely contribute to the different levels of selective pressure in c. resedifolia compared to c. impatiens lineages.
in addition, extrinsic factors  seem to be the primary drivers of heterogeneity in the levels of selective pressure observed among genes in c. resedifolia.

 organisms adapt to different habitats through natural selection, which favors the fixation of alleles that increase the fitness of the individual that bears them. however, it is quite difficult to identify the locus/loci targeted by selection. one reason is that the number of loci involved in a particular adaptation and their phenotypic effects vary depending on the genetic architecture underlying the adaptive trait. another reason is that often the sequence and/or information about the gene involved in the adaptive response is unavailable. two main approaches are used to identify genetic signatures of adaptive evolution and link them to phenotypic traits. the first is the candidate gene approach, whereby signatures of positive selection are identified using a population genetics framework in those genes previously known to be involved in the phenotypic trait of interest. the advantage of this approach is that the consequences of molecular variation on phenotypes can be inferred and their adaptive significance evaluated. on the other hand, this approach assumes a comprehensive knowledge of the function of such genes and may neglect other genes relevant to the phenotypic trait. the second is the genome-wide approach, where the pattern of molecular evolution of hundreds to thousands of genes scattered across the genome are analyzed simultaneously. as high throughput technologies are becoming cheaper and more accessible, this approach is increasingly gaining attention from the evolutionary biology research community. an important advantage of this method is that it does not rely on prior information about gene sequence, and therefore it is well suited for studying non-model species. the possibility of performing genome scale studies in non-model organisms is indeed a powerful approach to address the genetic basis of specific adaptations, which can only be obtained by choosing the appropriate organism and ecological context. however, the typically poor knowledge about gene function in non-model organisms often prevents a comprehensive understanding of the adaptive significance of eventual signatures of positive selection. a useful compromise is to use the genome-wide approach in species closely related to well-studied model organisms, so that gene function can be inferred by comparative analyses. in this way it is possible to exploit what is known about the ecology and the life history of the species, and thus the approach is particularly suited to identifying genes involved in species-specific and habitat-specific adaptations. environmental variation is a key factor driving adaptive evolution and determining the ecological niche of a species. plants and other sessile organisms are particularly affected by circadian and seasonal oscillations in abiotic factors. for example, sudden drops in temperature, high levels of solar irradiation and limited access to water are common sources of environmental stress, especially at high altitude in mountainous regions. these stressors affect the evolution of species living in these habitats, and their capacity to adapt to these stressors ultimately determines their distribution. to cope with drops in temperature, plants have developed a series of physiological adaptations that rely on the up- and down-regulation of cold responsive genes triggered by cold exposure. similarly, cold temperatures and high irradiation are not favorable to efficient photosynthesis, and consequently, a suite of photoprotective strategies are required for survival and reproduction at high altitudes. a fairly detailed understanding of the relevant regulatory pathways and gene function in the model species arabidopsis thaliana now exists; however, little is known about their adaptive role, particularly in relation to the diverse habitats present along an altitudinal gradient. in order to assess the adaptive role of cold responsive genes, as well of the genes involved in photosynthesis, we compared patterns of gene evolution in congeneric species living at different altitudes. we performed molecular evolution analyses in two closely related brassicaceae species adapted to non-overlapping altitudinal ranges: cardamine resedifolia, a perennial species usually growing under conditions of high irradiation and severe temperature oscillations between and meters above sea level; and c. impatiens, an annual nemoral species normally growing between and meters above sea level. the distinct habitats associated with high and low altitudes make c. resedifolia and c. impatiens ideal species for studying the adaptive significance of the genes involved in altitude-related stress responses. it should also be noted that the two species may differ in their outcrossing rates, even though conclusive evidence is still lacking: c. impatiens is mainly selfing, although some populations have been reported as partially outcrossing; c. resedifolia is also a predominantly selfing species, but the precise outcrossing rate is unknown. the two cardamine species described above are closely related to a. thaliana, making it possible to apply the extensive knowledge about the gene functions of this model organism to our study. gene sequences of both c. resedifolia and c. impatiens were obtained by high-throughput sequencing technology and subsequently identified and partitioned into functional classes based on gene similarity to the a. thaliana orthologues. then, we used complementary approaches, such as analyses of non-synonymous and synonymous substitutions and of their ratio, of levels of selective pressure, codon usage, and gene expression, to quantify the difference in adaptive evolution between functional classes and between species. this allowed us to infer the adaptive significance of genes involved in adaptation to cold stress and photosynthesis at high altitude and put them into an ecological context.  we obtained sequences representing the c. resedifolia and c. impatiens transcriptomes using high-throughput sequencing. genes and their putative function were identified by comparison of the sequences obtained here to the available annotated a. thaliana gene sequences. to minimize the chance of mistaking a paralogue for an orthologue, we considered as triplets of putative orthologues only those consisting of reciprocal best hits, i.e. those where the three sequences were consistently found as best hit matches of one another. eventually, we obtained triplets of nuclear genes, with a mean length for the a. thaliana orthologues of ± base pairs, corresponding to a mean coverage of ± % of the full-length gene sequence. in c. impatiens, the mean sequence length was ± bp, while in c. resedifolia the mean sequence length was ± bp. we partitioned these genes according to their putative function, and then focused our analyses on those functional classes that are associated with the adaptive response to high altitude. in particular, according to the annotation in a. thaliana, genes were involved in cold acclimation, were involved in uv-b and high irradiation response, were involved in general stress responses. to these three classes we added a manually compiled class that included genes functionally characterized as cold responsive genes. rate of molecular substitution in cardamine the analysis of the rates of nucleotide substitution revealed faster molecular evolution along the c. resedifolia lineage than along the c. impatiens lineage. because the rates of sequence evolution were highly correlated with the length of the a. thaliana orthologous gene, we corrected these measures accordingly by using the residuals of such correlations. while the rate of synonymous substitution, ds, was similar between lineages, the rate of non-synonymous substitution, dn, was larger in the c. resedifolia lineage than in the c. impatiens lineage. conversely, the ratio ω = dn/ds, was significantly larger in the c. impatiens than in the c. resedifolia lineage. c. impatiens c. resedifolia n a n a p c a number of genes. b mean and standard error. c wilcoxon test probability after correcting for gene length for each functional class, both the dn and the ds rates, as well their ratio ω, were not significantly different between the two lineages after correcting for multiple testing. selective constraints varied however across functional classes within the same lineage: dn was significantly lower in pgo than in other genes in the c. resedifolia lineage, and was significantly higher in crg than in other genes in both cardamine lineages. the statistically significant differences were not due to random effects resulting from the small sample size of the functional classes relative to the full dataset. in fact, a bootstrap approach revealed that at most % of the randomized datasets produced wilcoxon tests p values lower than those observed. since there was a significant correlation between dn and the length of the a. thaliana orthologue, we also re-analyzed the comparisons using the residuals of the correlation. the within-lineage differences were still statistically significant, thus excluding the possibility that such differences were a result of length heterogeneity across our gene sets. similarly, the mean values of the residuals for dn and ds still did not differ between lineages, with the exception of the pgo genes, where the dn residuals were lower in c. resedifolia than in c. impatiens. the analysis of the levels of selective pressure, ω, provided compelling evidence that different selective forces are dominating the evolution of the four functional categories in the two cardamine lineages. again, as a result of the correlation between substitution rates and gene length, we report the results of the residuals of such correlations. in the c. resedifolia lineage, the mean value of ω was significantly lower in pgo than in other genes, consistent with intense purifying selection. this is in sharp contrast with what observed in crg, where the mean value of ω was significantly higher than in other genes, indicating either relaxed or positive selection. instead, in the c. impatiens lineage, the levels of selective pressure in genes of the four functional classes were similar to those at the genome-wide level. in c. resedifolia, selective pressures were different between genes involved in cold response classified as cgo and those classified as crg. the two datasets have genes in common, whose mean ω did not differ from that of the remaining genes. however, the genes unique to cgo and crg had a mean ω in line with what was observed in the complete datasets. that is, genes unique to the cgo category had a mean ω significantly lower than the remaining genes, confirming strong purifying selection in cgo. conversely, genes exclusive to the crg category had a mean ω significantly higher than the other genes, indicating either relaxed or positive selection in these genes. the rate of evolution and the occurrence of positive selection were further investigated by analyzing the pattern of molecular substitution of each single gene. the most sensitive test for this purpose is one that employs branch-site models, which aim to detect positive selection affecting a few sites along particular lineages. when c. impatiens was used as the foreground lineage the test identified one outlier with fdr <. this number increased substantially when the c. resedifolia lineage was used instead, with seven genes showing evidence for positive selection along this lineage. two genes among these eight outliers were identified as being involved in general stress responses. a likelihood ratio tests in which the alternative model explained the pattern of codon substitution better than the null model, and with a false discovery rate lower than. patterns of codon substitution were tested using branch models, site models, and branch-site models. b tests compared models m and m0'. s tests compared model m2a to m1a and model m to m. bs tests were based on the branch-site model a, test subscripts in the b and bs tests indicate whether the c. impatiens or the c. resedifolia lineages were used as foreground branches. b gene of the sgo functional class. to increase the chance of discovering interesting candidate genes, we also employed additional maximum likelihood ratio tests based on branch or site codon substitution models. these tests identified several genes whose pattern was better explained by the occurrence of non-neutral evolution. branch models identified three genes with a pattern of nucleotide substitution that fit a model allowing two ω rates, one for the focal c. impatiens branch and a second for the rest of the phylogenetic tree. a single gene was detected if the c. resedifolia branch was used as the focal branch instead. using site models, the fit of the data under a neutral model was compared to that under a model of positive selection. in the first comparison, the neutral model admitted one class of sites with < ω < while the selection model allowed the presence of a second class of sites with ω >. the second, more powerful, approach compared more realistic models where neutrally evolving sites were partitioned in ten classes with ω values drawn from a beta distribution. the two likelihood ratio tests consistently identified four genes showing evidence for the action of positive selection. as expected, the tests based on either the branch, site or branch-site models resulted in the detection of sets of outliers with little overlap. among the three genes detected by more than one test, two were detected by both the b and the bci tests, and one was detected by the s s and bsci tests. to exclude the possibility that substitution rate patterns were the result of orthology mis-assignment, effectively inflating interspecific divergence, we compared homologous sequences from a. thaliana, cardamine and other plant species within a phylogenetic framework. specifically, we searched for putative paralogues of the gene triplets, for each of the candidate genes. in six such cases we verified the presence of putative paralogues using blast searches. the a. thaliana and the two putative cardamine orthologues identified by the reciprocal best-hit approach always clustered together in the tree, supporting their orthology and the results of the likelihood ratio tests. expression breadth, expression level and rate of molecular substitution the breadth of expression of a gene, namely whether it has tissue-specific expression or is expressed broadly over more tissues, has been shown to be correlated with several patterns of molecular evolution. to evaluate such effects in cardamine, we correlated the levels of selective pressure, ω, to the breadth of expression measured across tissues, across developmental time points, and along stress treatments. for these analyses we used the expression patterns of a. thaliana orthologues as proxies for expression in cardamine, assuming that orthologous genes will maintain similar expression patterns. these analyses found strong evidence that, in both cardamine lineages, genes evolved at different rates and under different selective regimes depending on their breadth of expression. the ratio ω decreased significantly with the increase in the temporal breadth of expression, i.e. the number of flower and leaf developmental stages at which a gene was expressed. the same statistically significant trend held when correlating ω with the spatial breadth of expression, i.e. the number of tissues/organs where the gene is expressed and the associated organ specificity index τ. because gene length was significantly correlated with both ω and expression breadth, we repeated our analyses using the residuals of such correlations. the partial correlations coefficients were still statistically highly significant and indicated that the correlations between ω and breadth of expression were independent from gene length. thus, a broad breadth of expression imposes selective constraints limiting the possible action of positive selection. conversely, genes expressed at one developmental time point or in specific organs evolve faster, as a result of either relaxed or positive selection. c. impatiens c. resedifolia * spearman's p <, ** p <, *** p <, **** p <, ns = non significant. ^ partial correlations are significant after correcting for the length of the a. thaliana orthologue. § significant after holm-bonferroni multiple test correction a breadth of expression can be either spatial or temporal. note that the organ-specificity index τ is inversely correlated with the number of organs at which a gene is expressed. when ω was compared to the duration of the stress response, we obtained mixed results that depended on the stress type. in general, the ratio ω was positively correlated to the persistence of gene expression during stress responses. after a holm-bonferroni correction for multiple testing, ω was significantly correlated to the duration of the responses to salt and uv-b in the c. impatiens lineage; and to the duration of the responses to salt, osmotic and cold stress in the c. resedifolia lineage. in addition to the breadth of expression, the level of gene expression was also an important determinant of ω. for instance, there was a statistically highly significant negative correlation between the mean gene expression and ω in both cardamine lineages. this correlation also remained highly significant when the mean expression level was controlled for the effects of gene length. thus, we conclude that selective constraints may limit the evolution of proteins encoded by highly expressed genes in cardamine. c. impatiens c. resedifolia n b p d p d mean expression maximum expression * spearman correlation p <, ** p <, *** p <, **** p <, ns = non significant. a gene functional class. b number of genes. c mean/maximum expression was estimated for all genes and genes included in each functional class and tested against mean for genes not in that functional class. d spearman's correlation between level of gene expression and ω estimated along the two cardamine lineages. the relationship between gene expression and rate of evolution indicates that the different selective pressures we observed in genes involved in stress response may be correlated to their expression level. for instance, genes of the cgo and pgo functional classes had a calculated expression level higher than other genes, and also had lower ω along the c. resedifolia lineage than the genome-wide mean. to investigate the association between gene expression and ω, we then calculated the residuals of their correlation and compared them across functional classes. the corrected ω were no longer statistically significant for pgo and cgo genes, suggesting that the levels of selective pressure in these genes were highly associated with their high expression levels. interestingly, removing the effect of expression level produced a statistically significant difference in ω between sgo genes and the other genes. such differences suggests heterogeneity in the association between rates of molecular evolution and expression levels in sgo genes, likely due to the diverse nature of the stress responses to which these genes are responding. results did not change if we analyzed the correlation between maximum, rather than mean, expression level and ω, further indicating a strong association between the expression pattern of a gene and its rate of molecular evolution. codon usage in cardamine genes synonymous codon usage can be under weak selection and lead to non-random use of the codons coding for the same aminoacid. due to hill-robertson interference, such bias is expected to correlate with the occurrence of positive selection, so that it will be lower in genes that experience episodes of positive selection. consistent with this prediction, in both cardamine lineages there was a statistically significant negative correlation between ω and codon bias, measured as fop. this correlation also held when controlling for the effects of other determinants of codon usage, such as expression levels, gene length and gc content. for instance, in both cardamine lineages, fop was significantly correlated with expression level, gene length, and gc content at the third codon position. however, partial correlation coefficients between codon usage and ω remained highly significant. we then analyzed codon usage for the genes in each of the four functional classes. genes involved in stress responses generally had larger codon bias compared to other genes. the difference was statistically significant for cgo and sgo genes in both cardamine lineages, while it was not significantly different for crg and pgo genes. however, such differences were no longer significant when controlling for gene expression level, suggesting that the expression pattern influenced codon usage more than the distinct selective pressures to which they are subject. thus, codon usage bias does not support the frequent action of positive selection in genes of the functional classes under consideration. discussion our analyses of protein-coding sequences highlighted lineage- and gene-specific evolutionary patterns in two species of cardamine characterized by distinct habitat preferences, life-histories and, possibly, breeding systems. c. resedifolia and c. impatiens are not sister species; hence, the substitution rates estimated in the present study apply to the lineages leading to the two species after the split from their common ancestor, rather than to these specific taxa as such. however, c. resedifolia, c. impatiens, and the clades they belong to, are characterized by divergent phenotypic and life-history traits that have profoundly shaped their evolutionary histories; for instance, the clade comprising c. resedifolia is basal to the cardamine radiation and includes only perennial species living in high altitude alpine habitats with well-developed petals indicating the existence of a mixed mating system. in contrast, c. impatiens is, to our knowledge, the only annual selfing species of its group, which includes species that are mainly found at moderate elevations. c. impatiens is further characterized by having very reduced or no petals, in agreement with a predominantly selfing mating system. thus, the observed patterns of molecular evolution reflect the differences in life history traits and habitats between the two lineages. the effects of such differences on the rate of molecular evolution are quite complex. for instance, the rate of non-synonymous substitution was significantly higher for the c. resedifolia lineage than for that of c. impatiens. previous studies have reported a higher substitution rate in annual plants compared to perennial plants. however, our observations are in contrast with this prediction, since the annual c. impatiens had lower substitution rate than the perennial c. resedifolia. moreover, since the rate of synonymous substitution was similar between lineages, there is no support for the hypothesis that mutation rate is higher along the c. resedifolia lineage than along that of c. impatiens. the interspecific differences in non-synonymous substitution rates could also be the result of differences in effective population size. for instance, low effective population size can affect the substitution rate by allowing slightly deleterious mutations to escape purifying selection and reach fixation. alternatively, high effective population size can increase the efficiency of positive selection by facilitating the fixation of favorable alleles in genes undergoing adaptive evolution. thus, a high dispersion in the distribution of dn across the genome is expected to be associated with a large effective population size, especially when there is recombination. the analysis of non-synonymous substitution fixation rate revealed that the variance in dn was significantly larger in c. resedifolia than in c. impatiens. this difference in protein evolution heterogeneity is thus consistent with more efficient positive selection and higher effective population size in c. resedifolia. both of these factors, in turn, would explain why dn is higher in this lineage than in c. impatiens. preliminary analyses of recombination levels and polymorphism suggest indeed that selfing rates are lower in c. resedifolia, and the effective population size is slightly larger than in c. impatiens. intriguingly, mean ω was larger in c. impatiens than in c. resedifolia. the fact that in c. impatiens ω was considerably lower than one is compatible with either moderately positive or relaxed selection. three lines of evidence suggest a reduced efficiency of purifying selection in c. impatiens, rather than pervasive positive selection at a genome-wide level. firstly, the analyses on dn suggest that positive selection was more efficient along the c. resedifolia lineage than the c. impatiens lineage. secondly, only one gene exhibited strong evidence for site-specific positive selection along the c. impatiens branch. by comparison, seven genes exhibited strong evidence for site-specific positive selection along the c. resedifolia branch. finally, the analysis of codon usage failed to detect differences in the efficiency of positive selection between the two lineages. however, the contrast between rates of non-synonymous substitution and of ω warrants some caution in drawing firm conclusions. the most likely explanation for the opposite patterns observed in dn and dn/ds lies on the high heterogeneity in substitution rate among genes. for instance the correlations and the coefficients of determination in the regressions between ds and ω and between dn and ω, suggest that substitution rates are relatively weak predictors of the actual selective pressure experienced by a gene. the statistically significant correlation between synonymous substitution rate and ω is rather unusual: previous studies established that the correlation is caused by a combination of adjacent substitutions and of a bias on substitution rate estimation due to either low divergence and/or statistical methods. appropriate analyses will be necessary to unravel the effects of such factors in our system. the complex dynamics of substitution rates are also evident in the statistically significant correlation to gene length. the influence of gene length on sequence evolution had been previously reported for populus tremula, stressing the importance of accounting for all diverse determinants of levels of gene and protein evolution. additional studies, including analyses on intraspecific polymorphism, will be certainly necessary to disentangle the neutral and selective forces that have shaped such patterns. the molecular evolution analyses in stress-related genes also revealed important lineage-specific patterns that may be associated with the distinct life-history traits and habitats of c. resedifolia and c. impatiens. for instance, similar selective regimes affected the evolution of genes involved in stress response in the c. impatiens lineage. conversely, there was a correlation between the type of stress response and the rate of molecular evolution along the c. resedifolia lineage; for example, genes involved in photosynthesis evolved slower than other genes, consistent with selective constraints that limited the accumulation of non-synonymous substitution. this makes sense given their involvement in a process that is particularly relevant in the high altitude habitat of c. resedifolia. strangely, in c. resedifolia selection acted in opposite directions in the two functional classes associated with cold responses. that is, compared to the genome-wide levels of selective pressure, genes that were identified as involved in cold response based on functional studies displayed significantly lower levels of selective pressure, while cold responsive genes were under more selective constraints. since there was little overlap between the two classes, the discrepancy in the levels of selective pressure between cgo and crg genes may have captured different selective pressures acting along the cold response pathway. however, the difference may also stem from the approaches used to define the two cold-response functional classes. specifically, the gene ontology annotations of the cgo genes were taken from the tair database, which is a less accurate source of functional information than the direct assays used to define crg genes. for instance, the tair annotations may be based solely on sequence or structural similarity, and a gene may have been annotated as cold responsive because it is indirectly up- or down-regulated in plants exposed to low temperatures. however, without a comprehensive knowledge on the response pathway it is difficult to evaluate their role in cold response in cardamine, especially considering that their expression patterns may differ from those of a. thaliana. the levels of selective pressure were also associated with the duration and pattern of gene expression. for instance, in c. impatiens and c. resedifolia, ω was lower for genes that were expressed only briefly during specific stress responses than for genes that were up-regulated over a longer period of time. most importantly, the correlation existed only for specific stress responses in each species, suggesting habitat-specific selective pressures. in particular, in c. resedifolia, ω was significantly affected by the duration of the responses to osmotic, salt and cold stress. the responses to these three stresses involve partially overlapping pathways, as cold stress causes membrane leakage and, as a consequence, the activation of the physiological responses also observed in high salt and osmotic stresses. these results support the adaptive relevance of these genes in response to the severe temperature changes that this species experiences in alpine habitats. on the other hand, in c. impatiens ω was significantly correlated with the extent of up-regulation in those genes involved in uv-b stress response. c. impatiens is a nemoral species that is likely to be particularly sensitive to the effect of exposure to uv-b light. therefore, it is reasonable to expect that this species has evolved response mechanisms based on gene up-regulation to cope with discontinuous uv stress. it will be extremely interesting to experimentally measure the expression of the cardamine genes in the functional classes considered, as species-specific adaptive variations in expression levels cannot be excluded. while unequal selective pressures across the stress-related functional classes were well-documented here, we detected positive selection in only a few single genes. in plants, previous studies have found evidence of genome-wide positive selection in some species, but in most cases there was little indication of widespread adaptive evolution. other studies have identified positive selection of some genes involved in stress response. the most likely explanation for the low rate of adaptive evolution is that most plants have low effective population size, which ultimately diminishes the efficiency of positive selection. estimation of the effective population size of c. resedifolia and c. impatiens will be useful to verify whether this parameter can explain the scarce evidence for positive selection in our dataset. four methodological issues may have further resulted in an underestimation of the genes targets of positive selection. the first is that the power of codon substitution models is fairly conservative compared to that of models incorporating polymorphism data. the second is that our reciprocal best-hit approach was prone to miss genes with high sequence divergence caused by adaptive evolution, as it was intended to avoid an overestimation of divergence by comparing paralogues. in particular, this approach may have overlooked duplicated genes, which are common in a. thaliana and in other plants, and which can undergo sub- or neo-functionalization driven by positive selection. for instance, our dataset included only ~11% of all a. thaliana genes, and we analyzed between 17% and 48% of the genes annotated as involved in the stress responses considered in the present study. an examination of our dataset revealed that many of the well-characterized transcription factors involved in cold response were absent. this indicates either that we are missing many of the genes upstream of the stress response pathways, or that in cardamine these transcription factors may be differently regulated than in arabidopsis thaliana or do not participate to the cold response. however, it is unlikely that this pathway is missing, since cold response is well-conserved across distantly related taxa. a third issue is related to the use of partial genes, which may have reduced the power of the likelihood ratio tests as a result of an insufficient number of informative sites. finally, genes expressed at a low level, including the aforementioned transcription factors, may be missing from our dataset because of insufficient coverage or normalization. because genes expressed at a low level are those evolving faster, this bias in gene representation could contribute to the relatively low number of rapidly evolving genes identified in this study. deeper sequencing efforts could undoubtedly improve the situation by increasing both coverage and average transcriptome length. one c. impatiens gene and seven c. resedifolia genes showed signatures of positive selection under the branch-sites model of codon substitution. the c. impatiens gene orthologue of at4g has not been functionally characterized in a. thaliana, although it is known that the protein encoded by this genes displays homology to the hyaluronan/mrna binding protein family, a group of proteins binding both specific rna and a high-molecular-mass polysaccharide extremely abundant in the connective tissue and extracellular matrix of animals. specific studies will be necessary to understand its role in adaptive processes in the c. impatiens lineage. as for the c. resedifolia candidate genes, no functional information is available for the orthologue of at1g which codes for a protein of unknown function with homology to tolb, a protein involved in outer membrane stability and uptake of biomolecules in e. coli. instead, at3g is known to code for a protein putatively involved in metal ion binding, and is similar to proteins of the heavy metal transport/detoxification superfamily. heavy metal hyperaccumulation has been documented in several brassicaceae species and c. resedifolia has been reported to accumulate large quantities of nickel from the nickel-rich debris of the glacial till, where this species usually grows. therefore, it is possible that the signature of positive selection identified in this orthologue is related to high contents of heavy metals in soil experienced by c. resedifolia. a third gene, orthologue of at1g is a valine-trna synthetase. given the multiple metabolic pathways in which valine is involved, it would be highly speculative to associate this gene to adaptive processes in the c. resedifolia lineage. based on functional evidence from arabidopsis and other plant species, two other genes identified as putative targets of positive selection in c. resedifolia may play major roles in the response against bacterial pathogens and insect herbivory. the first gene, at5g is part of the jasmonate-zim-domain protein family, whose members are involved in response to wounding and herbivory. the other candidate is the orthologue of at1g which codes for the arabidopsis epithiospecifier protein. esp catalyzes the formation of simpler nitriles and epithionitriles from glucosinolates, thus modulating the release of isothiocyanate, a metabolite involved in herbivory and pathogen defense in brassicaceae. in two closely related boechera species, the level of glucosinolates is negatively correlated with elevation preferences and growth rates, but positively correlated with drought tolerance. however, the response to herbivory as a function of elevation is not uniform across plants, with some species experiencing more and other less damage with increase in elevation. interestingly, a positive correlation with levels of insect herbivory was found for c. cordifolia exposed to full sunlight, possibly resulting from moderate water stress associated with a different insect guild compared to shadowy environments. these observations provide the framework to experimentally test whether the signature of positive selection identified in at5g and at1g orthologues might be related to the higher exposition to sun, lower water availability or slower growth rate characterizing c. resedifolia as compared to c. impatiens. interestingly, the orthologues of at1g and at3g may also be related to the light regimes characterizing c. resedifolia and c. impatiens habitats. at1g codes for a cytosolic ascorbate peroxidase that scavenges hydrogen peroxide in plant cells, thus reducing the accumulation of reactive oxygen species that cause cellular damage through protein oxidation. the product of at3g is a transcriptional activator of the growth regulating factor family. genes from this family have been demonstrated to be involved in the regulation of cell expansion and division in leaves, cotyledons and petals. strikingly, the product of at3g is one of the proteins oxidized in apx mutant plants in response to moderate light stress, indicating that it could also be part of the signaling cascade activated by ros stress. the site codon substitution models also identified putative targets of positive selection that deserve further characterization. in particular, the orthologue of at2g codes for a ribosomal protein that is involved in response to salt and genotoxic stress in a. thaliana. this gene is part of the ribosomal protein s family, which includes three paralogues with similar sequences and function. gene families can undergo relaxed selection immediately following duplication; however, this does not seem to be the case for at2g as the phylogenetic tree clearly indicates that the duplication occurred before the split between the arabidopsis and cardamine lineages. additional studies will be necessary to corroborate these findings and link the evolutionary pattern of each gene to its phenotypic effects. in particular, the use of intraspecific variation and functional analyses in the model species a. thaliana will be crucial to ascertaining whether positive selection or relaxed selection accelerated the evolution of these genes and their relevance in adaptive processes in cardamine.  overall, our results highlight the importance of employing complementary approaches to studying the genetic bases of adaptation in non-model species. our use of comparative genomics on congeneric species identified evolutionary patterns that aid the understanding of the extrinsic and intrinsic factors driving plant adaptation. in the case of cardamine, intrinsic factors most likely contribute to the different levels of selective pressure in c. resedifolia compared to c. impatiens lineages. in addition, extrinsic factors seem to be the primary drivers of heterogeneity in the levels of selective pressure observed among genes in c. resedifolia. 
high-throughput sequencing was used to obtain gene sequences from c. resedifolia and c. impatiens. using the available a. thaliana gene sequences and annotation, we identified nearly triplets of putative orthologues, including genes involved in cold response, photosynthesis or in general stress responses. we conclude that the selective pressures associated with the habitats typical of c. resedifolia may have caused the rapid evolution of genes involved in cold response. overall, our comparative genomic analyses revealed that differences in effective population size might contribute to the differences in the rate of protein evolution and in the levels of selective pressure between the c. impatiens and c. resedifolia lineages. here we studied molecular evolution on a genome-wide scale in two species of cardamine with distinct habitat preferences: c. resedifolia, found at high altitudes, and c. impatiens, found at low altitudes. our analyses revealed a statistically significant higher rate of molecular substitution in c. resedifolia than in c. impatiens, compatible with more efficient positive selection in the former. moreover, levels of selective pressure were heterogeneous between functional classes and between species, with cold responsive genes evolving particularly fast in c. resedifolia, but not in c. impatiens. by comparing estimated rates of molecular substitution, codon usage, and gene expression in these species with those of arabidopsis, we were able to evaluate the role of positive and relaxed selection in driving the evolution of cardamine genes. the within-species analyses also revealed evolutionary patterns associated with habitat preference of two cardamine species. conversely, the genome-wide level of selective pressure is compatible with more relaxed selection in c. impatiens. elucidating the selective and neutral forces underlying molecular evolution is fundamental to understanding the genetic basis of adaptation. our analyses focussed on genes that are involved in stress responses to two factors that differentiate the high- and low-altitude habitats, namely temperature and irradiation. molecular evolution next generation sequencingadaptive traitshabitat preferencecardamine
 elucidating the selective and neutral forces underlying molecular evolution is fundamental to understanding the genetic basis of adaptation. plants have evolved a suite of adaptive responses to cope with variable environmental conditions, but relatively little is known about which genes are involved in such responses. here we studied molecular evolution on a genome-wide scale in two species of cardamine with distinct habitat preferences: c. resedifolia, found at high altitudes, and c. impatiens, found at low altitudes. our analyses focussed on genes that are involved in stress responses to two factors that differentiate the high- and low-altitude habitats, namely temperature and irradiation.  high-throughput sequencing was used to obtain gene sequences from c. resedifolia and c. impatiens. using the available a. thaliana gene sequences and annotation, we identified nearly triplets of putative orthologues, including genes involved in cold response, photosynthesis or in general stress responses. by comparing estimated rates of molecular substitution, codon usage, and gene expression in these species with those of arabidopsis, we were able to evaluate the role of positive and relaxed selection in driving the evolution of cardamine genes. our analyses revealed a statistically significant higher rate of molecular substitution in c. resedifolia than in c. impatiens, compatible with more efficient positive selection in the former. conversely, the genome-wide level of selective pressure is compatible with more relaxed selection in c. impatiens. moreover, levels of selective pressure were heterogeneous between functional classes and between species, with cold responsive genes evolving particularly fast in c. resedifolia, but not in c. impatiens.  overall, our comparative genomic analyses revealed that differences in effective population size might contribute to the differences in the rate of protein evolution and in the levels of selective pressure between the c. impatiens and c. resedifolia lineages. the within-species analyses also revealed evolutionary patterns associated with habitat preference of two cardamine species. we conclude that the selective pressures associated with the habitats typical of c. resedifolia may have caused the rapid evolution of genes involved in cold response. molecular evolution next generation sequencingadaptive traitshabitat preferencecardamine 