BACKGROUND
the structure of the genetic code imposes constraints on the evolutionary variability of a nucleotide in a gene sequence, depending on its codon position.
the third codon positions are almost freely variable, because of the redundancy in the genetic code.
the excess variability at the third codon positions, after accounting for the replication-caused effects could be explained by mrna secondary structure, in mitochondrial genes  <cit> .
therefore, a third codon position nucleotide in a loop is more variable than that in a helix.
this suggests additional evolutionary constraints imposed by the mrna secondary structure, which motivated us to study its relationship with evolutionary variability.
note that the specific structure of the genetic code maximizes the capacity of coding regions to form helices in their secondary structure, while simultaneously also maximizing occurrences of off-frame stops  <cit> .
while maximizing numbers of off frame stops minimizes costs of ribosomal slippages  <cit> , the rationale for maximizing propensity for secondary structure is not known.
we suggest that it might protect from mutations, and test this hypothesis.
we chose a dataset of chloroplast mrnas for our secondary structure analyses.
this is because regulation of plastid gene expression during development and in response to light is of substantial interest in plant biology.
understanding the key factors underlying these processes can yield important insights that can be used to engineer and transform plastids  <cit> .
we use mfold, a free-energy minimization algorithm that uses sequence-dependent thermodynamic parameters, for evaluating the sub-optimal range of secondary structures for chloroplast mrnas.
we found that the secondary structure stability correlated negatively with the percentage of variable sites in the gene.
in order to gain further insight on the mrna secondary structure and sequence relationship, we decided to simulate sequences by mutating mrnas at each site according to their most likely substitution rate category.
this categorization for each site was calculated using the 'dnaml' tool from phylip v.  <dig> .
<dig>  package.
our simulation analyses revealed mrna size as an important factor, which acts as a flip-flop switch control on how these mutations affected the secondary structure.
transcript length has already been reported to be negatively correlated with gene expression level , positively with protein divergence in drosophila  <cit>  and metabolic cost minimization of protein synthesis  <cit> , and other factors relevant to molecular evolutionary ecology  <cit> .
we found that for shorter sequences, a greater percentage of the secondary structure of mutated mrnas is composed of helices, than that for the natural mrnas.
whereas for longer sequences, the secondary structure of mutated mrnas is composed of lesser percentage of sites in helices, than that for natural mrnas.
it is possible that secondary structure  usually protects functionally important sites from mutations.
this possibility is also ascertained by the fact that spontaneous mutations occur more rapidly on single-stranded dna than double-stranded dna  <cit> .
this could bias the mutational spectrum of a gene towards micro-adaptive changes, as mutations at less important sites are more likely to produce functional proteins with slightly altered optima.
the opposite is also a potential adaptation, to increase adaptability to more drastic environmental changes and in life-history.
the former fits better species considered as overall k-strategists because it increases micro-adaptations and niche-specialization, while the latter fits overall r-strategists, as its costs of decreasing survival are more likely to be bearable by organisms with high numbers of off-spring.
our results fit these general principles, and justify further adequate analyses.
this approach towards molecular processes and evolution converges with that described for r- and k-strategists at whole organism level by skulachev  <cit> .
RESULTS
secondary structural folding stability of nucleic acid sequences modulates their evolutionary variability
we tested for correlations between the length-adjusted residual stability measure for each gene  and that gene's relative fit to four pre-defined substitution rate categories .
this estimate of stability correlated positively  with the corresponding gene's percent fit to the slowest substitution rate category m <dig> and negatively with the percent fits to faster substitution rate categories m <dig>  and m <dig> .
correlations of residual rna stability with percent fit to the highest rate category m <dig> were also positive but not statistically significant  because only  <dig> out of  <dig> genes contain sites fitting m <dig> category.
overall, these results indicate that the chloroplast sense strand rna sequences with highly variable sites fold into less stable secondary structures, and corroborate similar findings in primate mitochondrial protein coding genes  <cit> .
the slight differences between the correlation coefficients for parametric and non-parametric analyses are mainly due to minimization of extreme data point effects by non-parametric analyses.
all p-values are according to 2-tailed t-tests.
the dna stability of the sense-strand sequence bore similar significant results with percent fits of the genes to two among four rate categories .
the correlation coefficients of dna stability with percent fits to m <dig>  m <dig>  m <dig> and m <dig> were  <dig> , p =  <dig> ; - <dig> , p =  <dig> ; - <dig> , p =  <dig>  and  <dig> , p =  <dig>  respectively.
non-parametric spearman's rank correlation analyses yielded significant results for percent fits to m <dig> , m <dig>  and m <dig>  but not to m <dig> .
it is not clear whether these subtle differences in the strength and significance of correlations for the rna and dna strands are due to an inherent difference in the robustness of the respective folding programs or indeed reflects different genuine evolutionary interactions between rna and dna secondary structures and sequences.
however, the similarity between the trends for rna and dna strands confirms that this correlation is due to effects specific to underlying single-strandedness.
dna-specific effects might associate with protection from mutations during single stranded periods, such as during dna replication  <cit>  or mrna transcription  <cit> .
the notably stronger rna-associated effects could relate to protection from decay  <cit> , as well as alterations of the sequence that would change the coding properties of the mrna.
variability distribution at synonymous and non-synonymous sites
in order to clearly interpret the above results in the light of synonymous/non-synonymous site variation, we analyzed the distribution of sites at all three codon positions among the four chosen substitution rate categories .
the percentage fits to these four categories were averaged across genes, individually for each codon position.
we have plotted these averages along with the standard deviation bars of these estimates.
the non-synonymous sites, i.e.
sites at first or second codon position, are more prevalent as m1-fitting sites.
second codon positions fit the m <dig> category better than first codon positions, as expected by its greater functional relevance  <cit> .
the third codon positions are relatively more prevalent as m2-, m3- and m4-fitting sites.
chloroplast mrnas are gc-rich, more at first and second codon positions, than at third codon positions, which contain more a's and t's .
the average gc content at the first and second codon position of mrnas correlates positively with their residual length-adjusted stability ) as shown in figure  <dig>  the gc content at the third codon position does not correlate with the residual stability.
overall, this is compatible with previous finding that a negative correlation exists between substitution rates at the first two codon positions versus the third codon position, for gc content  <cit> .
these results confirm that variability that negatively affects mrna stability primarily arises at synonymous sites because these sites tend to fit m <dig>  m <dig> and m <dig> more than the non-synonymous sites.
the argument becomes circular: the synonymous sites being gc-poor, a priori contributed lesser to mrna stability and thereby, have lesser evolutionary conservation constraints.
there is considerable variation among genes in the distribution of variability levels at the synonymous and non-synonymous sites, as indicated by the standard deviation bars.
our conclusion qualitatively remains the same even after separately analyzing this average distribution for individual genes .
however, the variability distributions for individual genes reveal slightly distinct patterns between the photosystem subunits and the ribosomal protein and atpase subunits.
in photosystem genes, 71% and 86% of first and second codon position sites fit m <dig>  whereas in ribosomal protein and atpase subunits, these are 55% and 63% respectively.
as a result, fits to the higher variability categories m <dig> and m <dig> were approximately  <dig> – <dig> and  <dig> – <dig> fold respectively, at the first and second codon positions, for the atpase and ribosomal protein subunits than for the photosystem genes .
local changes in secondary structure upon mutations in psilotum nudum and nephroselmis olivacea: effects of mrna size
we compared helix-forming propensities from nine mrnas of different lengths in psilotum nudum and nephroselmis olivacea among natural, mutated and randomized sequences .
these comparisons were performed after considering all the sites together as well as separately, based on the rate category of sites: m1-, m2-, m3-, and m4-fitting .
only three genes  contained sites that fit the m <dig> category, and therefore, it was not possible to show the trend in this case.
in general, for all types of sequences , the helix-forming propensities of mrnas increase with length of the mrna up to a given length and then saturate for longer mrnas, as indicated by the gradual logarithmic trend .
the overall integral stabilities  of mrnas follow a linear trend versus the size of the mrna  but the average percentage of sites in these mrnas that form helices correlate logarithmically with integral mrna stability .
this explains the saturation effect we observe in figure  <dig> .
upon mutating mrna sequences under the four site-specific substitution rate categories  and imposing transition: transversion bias of 2: <dig>  we found that the difference between the overall helix-forming propensities of mutated and natural mrnas varies according to their size.
the mutated mrna helix-forming propensity is greater than that of the natural mrnas for shorter transcript sizes, and the trend reverses for longer mrnas thus generating a crossover point in the correlation trend .
the difference between the natural and mutated helix-forming propensity curves is significant for all sites considered together  as well as for m1-fitting sites , but not for the sites fitting the remaining rate categories.
the differences between natural and random or between random and mutated helix-forming propensities are not significant for any of the cases, not even for the entire mrna.
this is particularly interesting, considering that mutating sequences results on average, 17% nucleotide changes, whereas randomizations result in 73% nucleotide changes  as compared to the original mrna sequence.
these length-dependent randomization effects on rna secondary structures might explain contradictions in the field on effects of sequence randomization on rna folding stabilities.
some report greater stability for natural than randomized mrna sequences  <cit> , while some report contrarily  <cit> .
protection from mutations at gene levels
the correlations between evolutionary variability of sites and stemminess of those sites for the  <dig> homologous genes in all seventeen chloroplast genomes are described in a table .
in a statistically significant majority of genes , conserved sites form more secondary structure than variable sites.
this tendency was significant for  <dig> specific genes , which is more than  <dig> times the  <dig>  statistically significant results expected according to the multiplicity of tests at p <  <dig> .
in four genes , the opposite was true: in most plant species, evolutionary conserved sites had lesser helix-forming propensity than other sites.
this number is also twice more than the expected number of false positive results at p <  <dig> .
hence, results suggest that in a majority of cases, sites coding for functionally important  protein regions are protected from mutations by mrna secondary structure , but that in specific cases, also the opposite occurs.
statistically, positive and negative associations between site conservation and formation of stems can be considered as existing and are not artifacts.
we did not detect significant differences between these associations for genes coded on the sense strand versus those coded on the anti-sense strand, not even for any particular species, despite the differences between the nucleotide contents for these strands.
the correlations between site variability and helix-formation was significant at p <  <dig>   in  <dig> among  <dig> ,  <dig> cases positive ) mrna-species combinations, which is twice the  <dig> % expected by pure chance for negative, but barely more  for positive correlations.
if hypotheses were available for expecting the direction of these correlations and one tailed tests were used,  <dig> and  <dig> cases would be significantly positive and negative, respectively.
this number is more than expected by pure chance for negative correlations, but not for positive ones.
however, it is clear that some mrnas have more significant correlations than expected, and others avoid them, including for positive ones.
indeed, a binomial distribution would expect in these conditions,  <dig> % of the mrnas to have no significant correlation in any of the  <dig> species, but in fact  <dig> %  of the mrnas had no significant correlations in any species.
on the other hand, there were several specific mrnas, for which more species bore significant  correlations than expected by a binomial distribution.
this result suggests that in some mrnas, evolutionary pressures tend to create such correlations, and in others, pressures prevent them.
considering that results suggest the existence of two forces, one occasionally causing positive, and another occasionally negative associations at the level of whole mrna sequences, one should not be surprised that in a larger than expected number of cases, balance between forces results in no or low, non-significant correlations between site-specific mrna secondary structure and variability levels.
these principles were also qualitatively correct while considering numbers of significant cases according to one tailed tests, separately for positive and negative associations.
hence, some specific mrnas tend towards negative , a smaller group towards positive , and a third group towards no association between secondary structure and site conservation .
investigating the rationale behind this trimodal distribution of mrnas is beyond the scope of this project, but suggests different adaptive strategies for different protein groups, and for different types of sites or amino acids.
results also suggest that the pressures causing associations between mrna secondary structure formation and site conservation are ubiquitous, and exist also in mrnas where no association is detected due to balancing forces.
protection from mutations at species levels
at the level of species, evolutionary variability of sites associated negatively with stemminess in a statistically significant majority of mrnas for  <dig> among  <dig> species .
at the level of single species, this tendency was significant at p <  <dig>   for chlamydomonas reinhardtii, chlorella vulgaris, mesostigma viride, physcomitella patens patens, triticum aestivum and psilotum nudum, which is more than the number of species expected by pure chance.
there was no species for which a significant number of mrnas had the opposite association.
calycanthus floridus was the most extreme in this direction, with 40% of its mrnas showing the positive association between evolutionary variability at a site and stem formation by that site.
there was no species for which no mrna had a significant correlation between site variability and helix-formation, which is less than what would be expected.
there were more species with no significant positive correlations between site variability and helix-formation than expected, which suggests that in some species, strong positive correlations are avoided .
there were lesser species with no significant negative correlations between site variability and stem formation than expected, which suggests that some species tend to have many strong negative correlations .
these results suggest that protein types and plant species vary with respect to correlations between site conservation and helix-formation.
in some cases, exposure to, rather than protection from mutations, is evolutionary adaptive, at least in some mrnas, such as rnapolb, where there was a positive correlation in  <dig> among  <dig> species, statistically significant at p =  <dig>   for  <dig> species.
we found that overall, protection from mutations decreases with gene length by calculating the correlation between gene length and the correlation coefficient between site variability and helix-formation for that gene.
this correlation was positive for  <dig> among  <dig> mrnas , and significant according to a two tailed test for  <dig> mrnas, among which none is negative.
hence, overall, for homologous genes, the extent of protection of conserved sites decreases with gene length.
a superficial look at these results could suggest that the lower correlations for longer mrnas are due to lower accuracy of mfold predictions for longer sequences.
according to this rationale, the longer mrna homologues would have the least negative t values for correlations between site variability and helix-formation, and these would be close to zero .
however, this rationale does not hold; in several mrna species where the decrease with gene length was significant, the longer mrna homologues do not have t statistics that are closer to zero or have low statistical significances, but they are often positive and sometimes, even statistically significant .
this result is not compatible with random errors due to folding inaccuracies, and hence the hypothesis that prediction inaccuracies increasing with sequence length can explain these results seems less likely than the one indicating that these results reflect biological phenomena.
it is possible that the decreased correlations for longer homologues reflect that their greater lengths are not due to additions of functionally necessary regions, such as catalytic sites, but perhaps, due to addition of sites with species-specific regulatory functions, or with no or little functional importance.
at the level of single species, across genes, the level of protection of conserved sites decreased with gene length in  <dig> species, significantly so, for spinacea oleracea and adiantum capillus-veneris.
this was more the case for larger "vascular" plants that can be considered broadly as k-strategists as compared to "micro" plants .
within vascular plants, these tendencies decreased with levels of domestication for  <dig> among  <dig> genes, significantly so at p <  <dig>  for  <dig> genes and increased significantly only in one mrna .
this suggests that artificial selection made genes of these relatively k-strategy oriented plants resemble more those of r-strategists.
overall, these tendencies are clearer while considering data across species for specific mrnas than while analyzing across mrnas for specific species, but the same principles are valid at both levels, although weaker for species-specific analyses.
hence, results suggest more mrna-specific adaptations than species-specific ones.
discussion
accuracy and inaccuracy of folding prediction programs
the assumption that secondary structure predictions are essentially accurate is central to this study, and hence it is worth dwelling on some evidence justifying it.
we will also discuss subsequently how the known inaccuracies of folding prediction programs are unlikely to bias our results.
real functional rnas are known to and also predicted to have fewer alternative structures than randomized sequences  <cit> .
this already establishes a level of realism in the details of the predictions.
a prediction accuracy of 90% reported for intermolecular hybridization of rnas, indicates at least a similar accuracy for our simpler case of intra-molecular rna hybridization  <cit> .
combinations of free energy minimization and comparative sequence analysis as used here, are more effective in enhancing structure prediction accuracy by finding a common low free energy structure  <cit> .
inaccuracy in secondary structure prediction algorithms is more prominent while having to predict a unique functional structure and is in fact minimized while considering a distribution of sub-optimal secondary structures.
nevertheless, at least in some cases where secondary structure was experimentally authenticated, specific local stem-loop structures are accurately predicted by folding programs for functionally important secondary structure components in chloroplast mrnas  <cit> .
in addition, structural predictions become more inaccurate with increase in the length of the sequences because of the concomitant increase in possible folds for longer rnas and in particular, when the contact distance between the base pairs exceeds  <dig> nucleotides  <cit> .
in theory, a boltzmann's ensemble of secondary structures, predicted by programs such as sfold  <cit>  and vienna  <cit>  would better reflect the population of structures for any given mrna.
however, the temperature of plants is widely determined by their environment.
the properties of the distribution of secondary structures are determined by temperature.
hence, precise calculations using boltzmann's distribution should integrate temperatures over the range and frequency encountered by each plant species in its natural environment.
therefore, while calculations using sfold's boltzmann distribution predictions are more precise and preferable for homeotherms, less accurate calculations apparently perform equally well for poikilotherms: predictions of trna structures by mfold and sfold in poikilothermic lizards yield similar results  <cit> .
there are other factors that are not accounted by classical secondary structure prediction algorithms.
one of them is the formation of protein-mrna complexes .
the other is the formation of non-canonical base pairs participating in the functional secondary structure of crucially important rnas  <cit> .
these considerations are important and not considering them probably adds noise to the analyses we present.
in other words, if it was possible to integrate in our analyses accurate information on secondary structure stabilization via complex formation with proteins, and/or on non-canonical base pairs, we would detect more statistically significant results than what we can detect now.
it is much less likely that including that information would prevent observing the patterns we describe.
however, if it was so, this would mean that these additional interactions  are designed to prevent the patterns we describe.
in that case, these would be two counter-balancing mechanisms, each important for rna function and worth describing.
it is also clear that even rnas that form complexes with proteins are only part of the time in such complexes, and hence our results might be relevant in this more restricted context for some of the rnas.
the relevance of non-canonical base pairings to our specific issue is lesser: such base pairings have been described until now mainly in rrnas and other rna molecules where rna and its secondary structure play a prime functional role.
it is unlikely that this is a major factor for typical mrnas, because their function is probably less confined to specific secondary structures, but rather more to group of alternative, closely related secondary structures with relatively low folding stabilities.
however, it makes sense to consider that in yet not described specific cases, non-canonical base pairings may play an important role in determining the secondary structure and therefore, function of the mrna.
the roles for mrna secondary structure that we detect, while considering the simpler scenario without non-canonical base pairings, only ascertain that further accurate secondary structure predictions would be very valuable to better understand rna biology, especially while integrating cellular and whole organism level physiologies with evolutionary perspectives.
secondary structural stability modulates sequence variability
we do not detect significant correlations between integral mrna stability and transcript's evolutionary variability.
this might be because evolutionary variability, quantified by a percentage is already size-adjusted, while integral stability is not.
however, integral and length adjusted residual stabilities increase with the absolute number of sites fitting each substitution rate categories  and  <dig>   for m <dig>  m <dig>  m <dig> and m <dig> respectively).
this result, for the least, stresses that the evolution of rna stability is in part explained by length, and mostly by the content and context-dependent component of the rna.
it is very unlikely that qualitatively our conclusions are due to artifacts , because correlations between evolutionary variability and both stabilities  yield similar results.
evolutionary variability  of nucleotide sequences decreases with residual, length-adjusted stability of the gene, when folded as rna and dna.
this would indicate that the helix-forming regions accumulate lesser mutations than loops  <cit> , which is in agreement with the positive correlation observed between presence of structural motifs and thermodynamic secondary structure stability  <cit> .
the m1-fitting sites are significantly over-represented at the first and second codon positions than at the third codon positions, whereas the m2-, m3-, and m4-fitting sites are relatively over-represented at the third codon positions .
the relatively high gc content at the first and second codon positions , explains well the positive correlation between the % m1-fitting sites and secondary structural stability.
similarly, the negative correlations of secondary structure stability with m2- and m3-fitting sites can be explained by the lower g+c content at the third codon positions.
the relative at richness at the third codon positions has been reported to be caused by codon usage bias in newly reported chloroplast genomes of nuphar advena and ranunculus macranthus species  <cit>  as well as in angiosperms like oryza sativa, zea mays, triticum aestivum and arabidopsis thaliana  <cit> .
codon usage bias is known to impose a variety of functional constraints by favoring increase mrna stability  <cit> , altered splicing  <cit> , hidden stop-codon formation  <cit>  and translational efficiency/gene expression often in relation to codon-anticodon matching frequencies from existing trna pool  <cit> .
it is also known to be controlled by genome-wide mutational constraints, probably due to underlying molecular processes  <cit> .
this suggests that codon usage bias may also be acting as a concerted component in shaping the positive association between rna secondary structural stability and evolutionary variability.
although our analyses yield similar negative correlations between sequence variability and folding stabilities for both rna as well as dna, we will be primarily discussing our results in the context of rna.
we do this because the steady state level of rna in a cell is likely to be much higher than that of single-stranded dna.
further, the dna folding program of mfold server  <cit>  is currently not as well developed as the rna version in terms of its incorporation of thermodynamic parameters and hence, stabilities obtained from the two versions are not comparable.
in a study testing the mutational robustness of over  <dig> naturally and artificially selected rna structures, mutation-proneness correlated negatively with thermodynamic stability in the selected rna molecules  <cit> .
our observed negative correlation of sequence variability with secondary structure stability is highly consistent with their results.
these non-intuitive results suggest the acquisition of dna mutations by the rna during transcription, thereby affecting its intrinsic stability, half-life, and ultimately, protein turnover.
these properties strongly suggest that rna stability is intimately connected with the functional regulation of the molecule in the cell.
effects of mrna length on mutational stability of secondary structure
we examined the effect of mrna length on its mutational stability in a randomly chosen species, psilotum nudum, and for assessing generality of results, a species distantly related to it, nephroselmis olivacea.
in both species, mrna length associates with mutational stability .
structural rnas such as ribosomal and nucleolar spliceosomal rnas and transfer rnas have a lower stability than a random rna sequence with the same dinucleotide frequency  <cit> .
in our analyses of the mutational stability of chloroplast mrnas and the non-randomness in intrinsic mrna stability, it is possible that the smaller mrnas behave like trnas or structural rnas and yield a similar stability response.
this speculation is because of the lower helix-forming propensity of the natural short sequences as compared to the mutated ones .
in drosophila, mrna abundance was found to be negatively correlated with transcript size  <cit> .
it is possible that these shorter chloroplast mrnas are also more abundant and their relatively low stability works in down-regulating their half lives, thereby enabling their pool sizes to stoichiometrically match with those of the lesser abundant but more stable longer transcripts in the same functional complexes.
another interesting component of this result is the logarithmic relationship of average helix-forming propensity of mrnas with overall mrna stabilities as well as mrna lengths .
this means that for mrnas with lengths beyond a certain limit, the excess length does not contribute towards increasing the helical component of the mrna, reflecting the need for a certain critical portion of the sequence in loops such that a dynamic structural state can be established vis-à-vis the function.
this result has to be considered with caution, because secondary structure predictions are less accurate for long sequences.
however, a similar saturating relationship between the integral stability of the mrna with the helix-forming propensity of that mrna suggests that this is less likely to be an artifact due to inaccuracies in folding prediction.
one could reason that the pattern we observe is species-specific and that other plant genomes might follow different trends, but then the probability of observing a similar response in two phylogenetically distant species is very low.
saturation levels of helix-forming propensity are consistently around 60–65% in both species.
based on this, we surmise that the logarithmic trend is likely to be general.
additionally, the thermodynamics of secondary structure folding are such that helix-formation decreases enthalpy of the system and is therefore, favored.
this is seen in the form of increasing helix-forming propensity with transcript length.
helix-formation, however, also decreases the entropy of the intermediate looped out regions, which becomes the limiting factor at certain mrna lengths.
further, the kinetics driven by the folding thermodynamics would less favor base-pairing between nucleotides distantly placed from one another on the linear sequence.
pairings between such distantly placed nucleotides might occur but would be unstable and therefore, will not hold.
this explains the saturation limit on helix-forming propensity in long mrnas.
indeed, studies have shown that conformational order present in evolved structures need not arise due to evolutionary optimization, but simply due to intrinsic folding rules of rna polymers  <cit> .
this is reiterated by the overall similarity between logarithmic trends in helix-forming propensities of mutated and natural mrnas versus their lengths, despite differences in the details of the trends.
protection from mutations
the simplest explanation for the positive correlation between secondary structure and evolutionary conservation of a site is that being part of helices protects those nucleotide regions from mutations.
this mechanism is most critical in the functionally important parts of the protein.
this protection may occur both in mrnas and, as indicated by results, single-stranded dna during replication.
there exists the possibility of a feedback loop: sites that are conserved due to functional constraints will tend to have more secondary structure.
this in turn will lead to fewer mutations and therefore, greater conservation.
this would result from stabilizing selection for optimal secondary structure because detrimental mutations  are selected against.
in addition, some sites should also be more constrained than others because they are complementary to functionally conserved sites.
hence, some highly conserved sites might be so, not because of their functional importance in the protein, but as part of packing-protection constraints in the mrna or dna during single stranded periods.
this probable effect at the level of dna suggests that substitutions are non-randomly distributed with respect to their functional impacts.
hence, the protecting mechanisms that we describe here bias the distribution of mutations towards favorable ones.
if correct, this mechanism results in a pattern of adaptive evolution where the secondary structure of nucleic acid sequences, prior to natural selection, filters mutations, a pattern that follows neo-lamarckian principles because the spectrum of potential mutations is a priori biased.
this likelihood of this mechanism especially increases while viewing our observation in light of the fact that spontaneous chemical changes occur at a greater rate on single- than double-stranded rna, quite analogous to mutations on dna  <cit> .
this means that secondary structure simply by virtue of causing certain regions to be duplex versus certain others, results in low and high mutation regimes.
this scenario, where secondary structure also involves evolutionary conservation, is therefore, a relatively adaptive one, as compared to one where this association was absent, where detrimental mutations are not only minimized but also certain beneficial mutations are allowed.
thereby, a mechanism evolves where structure biases the mutational distribution in a way that the positive effects of mutations are more likely.
we term this as 'lamarckian'.
it is interesting to note that this pattern could arise, as noted above, from stabilizing selection for an optimal secondary structure, which secondarily would result in a functionally biased mutation spectrum for rna or dna coding for proteins.
molecular aspects of life history strategies
we found three types of mrna species in respect to associations between site conservation and formation of secondary structure: those protecting conserved sites from mutations, a smaller group exposing them, and a group where presumably the pressures leading to the former two types are balanced, resulting in more than expected statistically non-significant cases.
this suggests that in all mrnas, and especially in the latter group, both types of strategies exist at different sites in the same mrna.
these might associate with amino acids and protein domains, as well as their function.
more detailed analyses at this level could reveal the adaptive strategies for different types of sites and proteins, whether adaptive evolution occurs there by drastic or gradual changes.
it is possible that protection of conserved sites promotes gradual adaptive evolution.
their exposure probably causes mutations that are very detrimental, majority of which are filtered out by natural selection and the remaining few which are likely to affect function, might promote saltatory evolution.
since ample sequence data are available, future explorative analyses focused on this issue will surely lead to deeper insights into mrna and protein functions, and their adaptive evolutionary strategies.
while at the level of mrna species, we do not yet understand the rationale behind the different types of correlations observed between site variability and stemminess, our preliminary analyses reveal some clues at the level of plant species.
however, note that in principle, analyses considering phylogenetic relatedness and larger species numbers are required for assessing the phenomena we preliminarily describe here below.
we find that for 66% of mrnas, 'r-strategist' species protect the conserved sites more than 'k' species.
this tendency is significant at p <  <dig>   for one mrna, psbl, but reversed in psbk.
these results indicate that minimization of mutations are higher in relatively smaller plants , following the general principles of minimization of metabolic costs  <cit> .
in the same plants, this tendency is particularly strong for large mrnas, which suggests that cost minimization is enhanced for synthesis of costly  proteins.
the opposite is true for plants which are relatively more k-strategists.
if one considers that exposing functionally conserved sites to mutations is a strategy that favors saltatory evolution, and that small proteins have more regulatory and large ones more catalytic, basic maintenance functions, then this suggests that r-strategists maximize evolutionary potential for regulatory functions, but are relatively conservative at the level of house keeping functions.
we find that the opposite is true for k-strategists.
these results are in line with the view that the large numbers of offspring in r strategists enable many unsuccessful evolutionary "experiments" to be done at the level of regulatory functions, while this approach is not tenable for large plants with low offspring numbers.
experimenting with large proteins is however more cost effective in k-strategists with large body sizes, because they have the reserves that enable them to cope with such costs, including the death of body parts, without killing the individual.
the results on associations with domestication level, which suggest that the molecular evolutionary strategy of naturally k-selected plants becomes more r-like by domestication, presumably because domestication increases reproductive output, are more speculative.
this is not only because our rankings of domestication levels between different species are difficult to justify in a formal way, but also because these analyses should be done separately on different groups of plants, such as eudicotyledons and others, or plants grown for their seeds versus those grown for other edible body parts.
however, as the majority of domesticated plants are selected for increased reproductive output, the convergence between natural r-strategists and domestication in k-strategists seems plausible.
it also strengthens the hypothesis that molecular adaptive evolutionary strategies  associates with relative investments in reproduction.
CONCLUSIONS
<dig>  high sequence variability, particularly at third codon positions, inversely correlates with length adjusted messenger rna stability, highlighting the functional and non-neutral aspects of synonymous site variation.
this is also in line with synonymous codon usage bias which favors at-richness at third codon positions.
<dig>  in silico mutations of mrnas results in destabilization of longer mrnas but not of shorter mrnas, which could explain apparently contradictory results in change in folding stabilities upon randomization of rnas.
<dig>  our analyses also reveal transcript length as an important factor, which coupled with evolutionary variability controls secondary structure stability.
this insight could substantially contribute to rna sequence and structure optimization studies that did not consider the effects of rna length.
<dig>  some highly conserved sites might be so not because of their functional importance in the protein, but as part of packing-protection constraints at the nucleic acid level.
<dig>  protection by secondary structure from mutations results in a pattern of adaptive evolution where the secondary structure of nucleic acid sequences prevents detrimental mutations, prior to or in addition to stabilizing selection on optimal secondary structures, a pattern that follows neo-lamarckian principles.
<dig>  in specific mrnas for some plants, functionally conserved sites are exposed to smutations, which might lead to saltatory adaptive evolution.
<dig>  protection from mutations at conserved sites is more widespread in r-strategists than k-strategists, perhaps because of greater mutation minimization pressures in smaller organisms.
<dig>  the above point is enhanced in r-strategists for large proteins, and for smaller proteins in k-strategists.
this might reflect adaptive strategies, the first leading to saltatory evolution in r-strategists, the other to gradual adaptive specialization in k-strategists.
<dig>  among k-strategists, domestication causes them to resemble r-strategists in terms of molecular adaptive strategies.

 the structure of the genetic code imposes constraints on the evolutionary variability of a nucleotide in a gene sequence, depending on its codon position. the third codon positions are almost freely variable, because of the redundancy in the genetic code. the excess variability at the third codon positions, after accounting for the replication-caused effects could be explained by mrna secondary structure, in mitochondrial genes. therefore, a third codon position nucleotide in a loop is more variable than that in a helix. this suggests additional evolutionary constraints imposed by the mrna secondary structure, which motivated us to study its relationship with evolutionary variability. note that the specific structure of the genetic code maximizes the capacity of coding regions to form helices in their secondary structure, while simultaneously also maximizing occurrences of off-frame stops. while maximizing numbers of off frame stops minimizes costs of ribosomal slippages, the rationale for maximizing propensity for secondary structure is not known. we suggest that it might protect from mutations, and test this hypothesis. we chose a dataset of chloroplast mrnas for our secondary structure analyses. this is because regulation of plastid gene expression during development and in response to light is of substantial interest in plant biology. understanding the key factors underlying these processes can yield important insights that can be used to engineer and transform plastids. we use mfold, a free-energy minimization algorithm that uses sequence-dependent thermodynamic parameters, for evaluating the sub-optimal range of secondary structures for chloroplast mrnas. we found that the secondary structure stability correlated negatively with the percentage of variable sites in the gene. in order to gain further insight on the mrna secondary structure and sequence relationship, we decided to simulate sequences by mutating mrnas at each site according to their most likely substitution rate category. this categorization for each site was calculated using the'dnaml' tool from phylip v.. package. our simulation analyses revealed mrna size as an important factor, which acts as a flip-flop switch control on how these mutations affected the secondary structure. transcript length has already been reported to be negatively correlated with gene expression level, positively with protein divergence in drosophila and metabolic cost minimization of protein synthesis, and other factors relevant to molecular evolutionary ecology. we found that for shorter sequences, a greater percentage of the secondary structure of mutated mrnas is composed of helices, than that for the natural mrnas. whereas for longer sequences, the secondary structure of mutated mrnas is composed of lesser percentage of sites in helices, than that for natural mrnas. it is possible that secondary structure usually protects functionally important sites from mutations. this possibility is also ascertained by the fact that spontaneous mutations occur more rapidly on single-stranded dna than double-stranded dna. this could bias the mutational spectrum of a gene towards micro-adaptive changes, as mutations at less important sites are more likely to produce functional proteins with slightly altered optima. the opposite is also a potential adaptation, to increase adaptability to more drastic environmental changes and in life-history. the former fits better species considered as overall k-strategists because it increases micro-adaptations and niche-specialization, while the latter fits overall r-strategists, as its costs of decreasing survival are more likely to be bearable by organisms with high numbers of off-spring. our results fit these general principles, and justify further adequate analyses. this approach towards molecular processes and evolution converges with that described for r- and k-strategists at whole organism level by skulachev.  secondary structural folding stability of nucleic acid sequences modulates their evolutionary variability we tested for correlations between the length-adjusted residual stability measure for each gene and that gene's relative fit to four pre-defined substitution rate categories. this estimate of stability correlated positively with the corresponding gene's percent fit to the slowest substitution rate category m and negatively with the percent fits to faster substitution rate categories m and m. correlations of residual rna stability with percent fit to the highest rate category m were also positive but not statistically significant because only out of genes contain sites fitting m category. overall, these results indicate that the chloroplast sense strand rna sequences with highly variable sites fold into less stable secondary structures, and corroborate similar findings in primate mitochondrial protein coding genes. the slight differences between the correlation coefficients for parametric and non-parametric analyses are mainly due to minimization of extreme data point effects by non-parametric analyses. all p-values are according to 2-tailed t-tests. the dna stability of the sense-strand sequence bore similar significant results with percent fits of the genes to two among four rate categories. the correlation coefficients of dna stability with percent fits to m m m and m were, p =; -, p =; -, p = and, p = respectively. non-parametric spearman's rank correlation analyses yielded significant results for percent fits to m, m and m but not to m. it is not clear whether these subtle differences in the strength and significance of correlations for the rna and dna strands are due to an inherent difference in the robustness of the respective folding programs or indeed reflects different genuine evolutionary interactions between rna and dna secondary structures and sequences. however, the similarity between the trends for rna and dna strands confirms that this correlation is due to effects specific to underlying single-strandedness. dna-specific effects might associate with protection from mutations during single stranded periods, such as during dna replication or mrna transcription. the notably stronger rna-associated effects could relate to protection from decay, as well as alterations of the sequence that would change the coding properties of the mrna. variability distribution at synonymous and non-synonymous sites in order to clearly interpret the above results in the light of synonymous/non-synonymous site variation, we analyzed the distribution of sites at all three codon positions among the four chosen substitution rate categories. the percentage fits to these four categories were averaged across genes, individually for each codon position. we have plotted these averages along with the standard deviation bars of these estimates. the non-synonymous sites, i.e. sites at first or second codon position, are more prevalent as m1-fitting sites. second codon positions fit the m category better than first codon positions, as expected by its greater functional relevance. the third codon positions are relatively more prevalent as m2-, m3- and m4-fitting sites. chloroplast mrnas are gc-rich, more at first and second codon positions, than at third codon positions, which contain more a's and t's. the average gc content at the first and second codon position of mrnas correlates positively with their residual length-adjusted stability ) as shown in figure the gc content at the third codon position does not correlate with the residual stability. overall, this is compatible with previous finding that a negative correlation exists between substitution rates at the first two codon positions versus the third codon position, for gc content. these results confirm that variability that negatively affects mrna stability primarily arises at synonymous sites because these sites tend to fit m m and m more than the non-synonymous sites. the argument becomes circular: the synonymous sites being gc-poor, a priori contributed lesser to mrna stability and thereby, have lesser evolutionary conservation constraints. there is considerable variation among genes in the distribution of variability levels at the synonymous and non-synonymous sites, as indicated by the standard deviation bars. our conclusion qualitatively remains the same even after separately analyzing this average distribution for individual genes. however, the variability distributions for individual genes reveal slightly distinct patterns between the photosystem subunits and the ribosomal protein and atpase subunits. in photosystem genes, 71% and 86% of first and second codon position sites fit m whereas in ribosomal protein and atpase subunits, these are 55% and 63% respectively. as a result, fits to the higher variability categories m and m were approximately – and – fold respectively, at the first and second codon positions, for the atpase and ribosomal protein subunits than for the photosystem genes. local changes in secondary structure upon mutations in psilotum nudum and nephroselmis olivacea: effects of mrna size we compared helix-forming propensities from nine mrnas of different lengths in psilotum nudum and nephroselmis olivacea among natural, mutated and randomized sequences. these comparisons were performed after considering all the sites together as well as separately, based on the rate category of sites: m1-, m2-, m3-, and m4-fitting. only three genes contained sites that fit the m category, and therefore, it was not possible to show the trend in this case. in general, for all types of sequences, the helix-forming propensities of mrnas increase with length of the mrna up to a given length and then saturate for longer mrnas, as indicated by the gradual logarithmic trend. the overall integral stabilities of mrnas follow a linear trend versus the size of the mrna but the average percentage of sites in these mrnas that form helices correlate logarithmically with integral mrna stability. this explains the saturation effect we observe in figure. upon mutating mrna sequences under the four site-specific substitution rate categories and imposing transition: transversion bias of 2: we found that the difference between the overall helix-forming propensities of mutated and natural mrnas varies according to their size. the mutated mrna helix-forming propensity is greater than that of the natural mrnas for shorter transcript sizes, and the trend reverses for longer mrnas thus generating a crossover point in the correlation trend. the difference between the natural and mutated helix-forming propensity curves is significant for all sites considered together as well as for m1-fitting sites, but not for the sites fitting the remaining rate categories. the differences between natural and random or between random and mutated helix-forming propensities are not significant for any of the cases, not even for the entire mrna. this is particularly interesting, considering that mutating sequences results on average, 17% nucleotide changes, whereas randomizations result in 73% nucleotide changes as compared to the original mrna sequence. these length-dependent randomization effects on rna secondary structures might explain contradictions in the field on effects of sequence randomization on rna folding stabilities. some report greater stability for natural than randomized mrna sequences, while some report contrarily. protection from mutations at gene levels the correlations between evolutionary variability of sites and stemminess of those sites for the homologous genes in all seventeen chloroplast genomes are described in a table. in a statistically significant majority of genes, conserved sites form more secondary structure than variable sites. this tendency was significant for specific genes, which is more than times the statistically significant results expected according to the multiplicity of tests at p <. in four genes, the opposite was true: in most plant species, evolutionary conserved sites had lesser helix-forming propensity than other sites. this number is also twice more than the expected number of false positive results at p <. hence, results suggest that in a majority of cases, sites coding for functionally important protein regions are protected from mutations by mrna secondary structure, but that in specific cases, also the opposite occurs. statistically, positive and negative associations between site conservation and formation of stems can be considered as existing and are not artifacts. we did not detect significant differences between these associations for genes coded on the sense strand versus those coded on the anti-sense strand, not even for any particular species, despite the differences between the nucleotide contents for these strands. the correlations between site variability and helix-formation was significant at p < in among, cases positive ) mrna-species combinations, which is twice the % expected by pure chance for negative, but barely more for positive correlations. if hypotheses were available for expecting the direction of these correlations and one tailed tests were used, and cases would be significantly positive and negative, respectively. this number is more than expected by pure chance for negative correlations, but not for positive ones. however, it is clear that some mrnas have more significant correlations than expected, and others avoid them, including for positive ones. indeed, a binomial distribution would expect in these conditions, % of the mrnas to have no significant correlation in any of the species, but in fact % of the mrnas had no significant correlations in any species. on the other hand, there were several specific mrnas, for which more species bore significant correlations than expected by a binomial distribution. this result suggests that in some mrnas, evolutionary pressures tend to create such correlations, and in others, pressures prevent them. considering that results suggest the existence of two forces, one occasionally causing positive, and another occasionally negative associations at the level of whole mrna sequences, one should not be surprised that in a larger than expected number of cases, balance between forces results in no or low, non-significant correlations between site-specific mrna secondary structure and variability levels. these principles were also qualitatively correct while considering numbers of significant cases according to one tailed tests, separately for positive and negative associations. hence, some specific mrnas tend towards negative, a smaller group towards positive, and a third group towards no association between secondary structure and site conservation. investigating the rationale behind this trimodal distribution of mrnas is beyond the scope of this project, but suggests different adaptive strategies for different protein groups, and for different types of sites or amino acids. results also suggest that the pressures causing associations between mrna secondary structure formation and site conservation are ubiquitous, and exist also in mrnas where no association is detected due to balancing forces. protection from mutations at species levels at the level of species, evolutionary variability of sites associated negatively with stemminess in a statistically significant majority of mrnas for among species. at the level of single species, this tendency was significant at p < for chlamydomonas reinhardtii, chlorella vulgaris, mesostigma viride, physcomitella patens patens, triticum aestivum and psilotum nudum, which is more than the number of species expected by pure chance. there was no species for which a significant number of mrnas had the opposite association. calycanthus floridus was the most extreme in this direction, with 40% of its mrnas showing the positive association between evolutionary variability at a site and stem formation by that site. there was no species for which no mrna had a significant correlation between site variability and helix-formation, which is less than what would be expected. there were more species with no significant positive correlations between site variability and helix-formation than expected, which suggests that in some species, strong positive correlations are avoided. there were lesser species with no significant negative correlations between site variability and stem formation than expected, which suggests that some species tend to have many strong negative correlations. these results suggest that protein types and plant species vary with respect to correlations between site conservation and helix-formation. in some cases, exposure to, rather than protection from mutations, is evolutionary adaptive, at least in some mrnas, such as rnapolb, where there was a positive correlation in among species, statistically significant at p = for species. we found that overall, protection from mutations decreases with gene length by calculating the correlation between gene length and the correlation coefficient between site variability and helix-formation for that gene. this correlation was positive for among mrnas, and significant according to a two tailed test for mrnas, among which none is negative. hence, overall, for homologous genes, the extent of protection of conserved sites decreases with gene length. a superficial look at these results could suggest that the lower correlations for longer mrnas are due to lower accuracy of mfold predictions for longer sequences. according to this rationale, the longer mrna homologues would have the least negative t values for correlations between site variability and helix-formation, and these would be close to zero. however, this rationale does not hold; in several mrna species where the decrease with gene length was significant, the longer mrna homologues do not have t statistics that are closer to zero or have low statistical significances, but they are often positive and sometimes, even statistically significant. this result is not compatible with random errors due to folding inaccuracies, and hence the hypothesis that prediction inaccuracies increasing with sequence length can explain these results seems less likely than the one indicating that these results reflect biological phenomena. it is possible that the decreased correlations for longer homologues reflect that their greater lengths are not due to additions of functionally necessary regions, such as catalytic sites, but perhaps, due to addition of sites with species-specific regulatory functions, or with no or little functional importance. at the level of single species, across genes, the level of protection of conserved sites decreased with gene length in species, significantly so, for spinacea oleracea and adiantum capillus-veneris. this was more the case for larger "vascular" plants that can be considered broadly as k-strategists as compared to "micro" plants. within vascular plants, these tendencies decreased with levels of domestication for among genes, significantly so at p < for genes and increased significantly only in one mrna. this suggests that artificial selection made genes of these relatively k-strategy oriented plants resemble more those of r-strategists. overall, these tendencies are clearer while considering data across species for specific mrnas than while analyzing across mrnas for specific species, but the same principles are valid at both levels, although weaker for species-specific analyses. hence, results suggest more mrna-specific adaptations than species-specific ones. discussion accuracy and inaccuracy of folding prediction programs the assumption that secondary structure predictions are essentially accurate is central to this study, and hence it is worth dwelling on some evidence justifying it. we will also discuss subsequently how the known inaccuracies of folding prediction programs are unlikely to bias our results. real functional rnas are known to and also predicted to have fewer alternative structures than randomized sequences. this already establishes a level of realism in the details of the predictions. a prediction accuracy of 90% reported for intermolecular hybridization of rnas, indicates at least a similar accuracy for our simpler case of intra-molecular rna hybridization. combinations of free energy minimization and comparative sequence analysis as used here, are more effective in enhancing structure prediction accuracy by finding a common low free energy structure. inaccuracy in secondary structure prediction algorithms is more prominent while having to predict a unique functional structure and is in fact minimized while considering a distribution of sub-optimal secondary structures. nevertheless, at least in some cases where secondary structure was experimentally authenticated, specific local stem-loop structures are accurately predicted by folding programs for functionally important secondary structure components in chloroplast mrnas. in addition, structural predictions become more inaccurate with increase in the length of the sequences because of the concomitant increase in possible folds for longer rnas and in particular, when the contact distance between the base pairs exceeds nucleotides. in theory, a boltzmann's ensemble of secondary structures, predicted by programs such as sfold and vienna would better reflect the population of structures for any given mrna. however, the temperature of plants is widely determined by their environment. the properties of the distribution of secondary structures are determined by temperature. hence, precise calculations using boltzmann's distribution should integrate temperatures over the range and frequency encountered by each plant species in its natural environment. therefore, while calculations using sfold's boltzmann distribution predictions are more precise and preferable for homeotherms, less accurate calculations apparently perform equally well for poikilotherms: predictions of trna structures by mfold and sfold in poikilothermic lizards yield similar results. there are other factors that are not accounted by classical secondary structure prediction algorithms. one of them is the formation of protein-mrna complexes. the other is the formation of non-canonical base pairs participating in the functional secondary structure of crucially important rnas. these considerations are important and not considering them probably adds noise to the analyses we present. in other words, if it was possible to integrate in our analyses accurate information on secondary structure stabilization via complex formation with proteins, and/or on non-canonical base pairs, we would detect more statistically significant results than what we can detect now. it is much less likely that including that information would prevent observing the patterns we describe. however, if it was so, this would mean that these additional interactions are designed to prevent the patterns we describe. in that case, these would be two counter-balancing mechanisms, each important for rna function and worth describing. it is also clear that even rnas that form complexes with proteins are only part of the time in such complexes, and hence our results might be relevant in this more restricted context for some of the rnas. the relevance of non-canonical base pairings to our specific issue is lesser: such base pairings have been described until now mainly in rrnas and other rna molecules where rna and its secondary structure play a prime functional role. it is unlikely that this is a major factor for typical mrnas, because their function is probably less confined to specific secondary structures, but rather more to group of alternative, closely related secondary structures with relatively low folding stabilities. however, it makes sense to consider that in yet not described specific cases, non-canonical base pairings may play an important role in determining the secondary structure and therefore, function of the mrna. the roles for mrna secondary structure that we detect, while considering the simpler scenario without non-canonical base pairings, only ascertain that further accurate secondary structure predictions would be very valuable to better understand rna biology, especially while integrating cellular and whole organism level physiologies with evolutionary perspectives. secondary structural stability modulates sequence variability we do not detect significant correlations between integral mrna stability and transcript's evolutionary variability. this might be because evolutionary variability, quantified by a percentage is already size-adjusted, while integral stability is not. however, integral and length adjusted residual stabilities increase with the absolute number of sites fitting each substitution rate categories and for m m m and m respectively). this result, for the least, stresses that the evolution of rna stability is in part explained by length, and mostly by the content and context-dependent component of the rna. it is very unlikely that qualitatively our conclusions are due to artifacts, because correlations between evolutionary variability and both stabilities yield similar results. evolutionary variability of nucleotide sequences decreases with residual, length-adjusted stability of the gene, when folded as rna and dna. this would indicate that the helix-forming regions accumulate lesser mutations than loops, which is in agreement with the positive correlation observed between presence of structural motifs and thermodynamic secondary structure stability. the m1-fitting sites are significantly over-represented at the first and second codon positions than at the third codon positions, whereas the m2-, m3-, and m4-fitting sites are relatively over-represented at the third codon positions. the relatively high gc content at the first and second codon positions, explains well the positive correlation between the % m1-fitting sites and secondary structural stability. similarly, the negative correlations of secondary structure stability with m2- and m3-fitting sites can be explained by the lower g+c content at the third codon positions. the relative at richness at the third codon positions has been reported to be caused by codon usage bias in newly reported chloroplast genomes of nuphar advena and ranunculus macranthus species as well as in angiosperms like oryza sativa, zea mays, triticum aestivum and arabidopsis thaliana. codon usage bias is known to impose a variety of functional constraints by favoring increase mrna stability, altered splicing, hidden stop-codon formation and translational efficiency/gene expression often in relation to codon-anticodon matching frequencies from existing trna pool. it is also known to be controlled by genome-wide mutational constraints, probably due to underlying molecular processes. this suggests that codon usage bias may also be acting as a concerted component in shaping the positive association between rna secondary structural stability and evolutionary variability. although our analyses yield similar negative correlations between sequence variability and folding stabilities for both rna as well as dna, we will be primarily discussing our results in the context of rna. we do this because the steady state level of rna in a cell is likely to be much higher than that of single-stranded dna. further, the dna folding program of mfold server is currently not as well developed as the rna version in terms of its incorporation of thermodynamic parameters and hence, stabilities obtained from the two versions are not comparable. in a study testing the mutational robustness of over naturally and artificially selected rna structures, mutation-proneness correlated negatively with thermodynamic stability in the selected rna molecules. our observed negative correlation of sequence variability with secondary structure stability is highly consistent with their results. these non-intuitive results suggest the acquisition of dna mutations by the rna during transcription, thereby affecting its intrinsic stability, half-life, and ultimately, protein turnover. these properties strongly suggest that rna stability is intimately connected with the functional regulation of the molecule in the cell. effects of mrna length on mutational stability of secondary structure we examined the effect of mrna length on its mutational stability in a randomly chosen species, psilotum nudum, and for assessing generality of results, a species distantly related to it, nephroselmis olivacea. in both species, mrna length associates with mutational stability. structural rnas such as ribosomal and nucleolar spliceosomal rnas and transfer rnas have a lower stability than a random rna sequence with the same dinucleotide frequency. in our analyses of the mutational stability of chloroplast mrnas and the non-randomness in intrinsic mrna stability, it is possible that the smaller mrnas behave like trnas or structural rnas and yield a similar stability response. this speculation is because of the lower helix-forming propensity of the natural short sequences as compared to the mutated ones. in drosophila, mrna abundance was found to be negatively correlated with transcript size. it is possible that these shorter chloroplast mrnas are also more abundant and their relatively low stability works in down-regulating their half lives, thereby enabling their pool sizes to stoichiometrically match with those of the lesser abundant but more stable longer transcripts in the same functional complexes. another interesting component of this result is the logarithmic relationship of average helix-forming propensity of mrnas with overall mrna stabilities as well as mrna lengths. this means that for mrnas with lengths beyond a certain limit, the excess length does not contribute towards increasing the helical component of the mrna, reflecting the need for a certain critical portion of the sequence in loops such that a dynamic structural state can be established vis-à-vis the function. this result has to be considered with caution, because secondary structure predictions are less accurate for long sequences. however, a similar saturating relationship between the integral stability of the mrna with the helix-forming propensity of that mrna suggests that this is less likely to be an artifact due to inaccuracies in folding prediction. one could reason that the pattern we observe is species-specific and that other plant genomes might follow different trends, but then the probability of observing a similar response in two phylogenetically distant species is very low. saturation levels of helix-forming propensity are consistently around 60–65% in both species. based on this, we surmise that the logarithmic trend is likely to be general. additionally, the thermodynamics of secondary structure folding are such that helix-formation decreases enthalpy of the system and is therefore, favored. this is seen in the form of increasing helix-forming propensity with transcript length. helix-formation, however, also decreases the entropy of the intermediate looped out regions, which becomes the limiting factor at certain mrna lengths. further, the kinetics driven by the folding thermodynamics would less favor base-pairing between nucleotides distantly placed from one another on the linear sequence. pairings between such distantly placed nucleotides might occur but would be unstable and therefore, will not hold. this explains the saturation limit on helix-forming propensity in long mrnas. indeed, studies have shown that conformational order present in evolved structures need not arise due to evolutionary optimization, but simply due to intrinsic folding rules of rna polymers. this is reiterated by the overall similarity between logarithmic trends in helix-forming propensities of mutated and natural mrnas versus their lengths, despite differences in the details of the trends. protection from mutations the simplest explanation for the positive correlation between secondary structure and evolutionary conservation of a site is that being part of helices protects those nucleotide regions from mutations. this mechanism is most critical in the functionally important parts of the protein. this protection may occur both in mrnas and, as indicated by results, single-stranded dna during replication. there exists the possibility of a feedback loop: sites that are conserved due to functional constraints will tend to have more secondary structure. this in turn will lead to fewer mutations and therefore, greater conservation. this would result from stabilizing selection for optimal secondary structure because detrimental mutations are selected against. in addition, some sites should also be more constrained than others because they are complementary to functionally conserved sites. hence, some highly conserved sites might be so, not because of their functional importance in the protein, but as part of packing-protection constraints in the mrna or dna during single stranded periods. this probable effect at the level of dna suggests that substitutions are non-randomly distributed with respect to their functional impacts. hence, the protecting mechanisms that we describe here bias the distribution of mutations towards favorable ones. if correct, this mechanism results in a pattern of adaptive evolution where the secondary structure of nucleic acid sequences, prior to natural selection, filters mutations, a pattern that follows neo-lamarckian principles because the spectrum of potential mutations is a priori biased. this likelihood of this mechanism especially increases while viewing our observation in light of the fact that spontaneous chemical changes occur at a greater rate on single- than double-stranded rna, quite analogous to mutations on dna. this means that secondary structure simply by virtue of causing certain regions to be duplex versus certain others, results in low and high mutation regimes. this scenario, where secondary structure also involves evolutionary conservation, is therefore, a relatively adaptive one, as compared to one where this association was absent, where detrimental mutations are not only minimized but also certain beneficial mutations are allowed. thereby, a mechanism evolves where structure biases the mutational distribution in a way that the positive effects of mutations are more likely. we term this as'lamarckian'. it is interesting to note that this pattern could arise, as noted above, from stabilizing selection for an optimal secondary structure, which secondarily would result in a functionally biased mutation spectrum for rna or dna coding for proteins. molecular aspects of life history strategies we found three types of mrna species in respect to associations between site conservation and formation of secondary structure: those protecting conserved sites from mutations, a smaller group exposing them, and a group where presumably the pressures leading to the former two types are balanced, resulting in more than expected statistically non-significant cases. this suggests that in all mrnas, and especially in the latter group, both types of strategies exist at different sites in the same mrna. these might associate with amino acids and protein domains, as well as their function. more detailed analyses at this level could reveal the adaptive strategies for different types of sites and proteins, whether adaptive evolution occurs there by drastic or gradual changes. it is possible that protection of conserved sites promotes gradual adaptive evolution. their exposure probably causes mutations that are very detrimental, majority of which are filtered out by natural selection and the remaining few which are likely to affect function, might promote saltatory evolution. since ample sequence data are available, future explorative analyses focused on this issue will surely lead to deeper insights into mrna and protein functions, and their adaptive evolutionary strategies. while at the level of mrna species, we do not yet understand the rationale behind the different types of correlations observed between site variability and stemminess, our preliminary analyses reveal some clues at the level of plant species. however, note that in principle, analyses considering phylogenetic relatedness and larger species numbers are required for assessing the phenomena we preliminarily describe here below. we find that for 66% of mrnas,'r-strategist' species protect the conserved sites more than'k' species. this tendency is significant at p < for one mrna, psbl, but reversed in psbk. these results indicate that minimization of mutations are higher in relatively smaller plants, following the general principles of minimization of metabolic costs. in the same plants, this tendency is particularly strong for large mrnas, which suggests that cost minimization is enhanced for synthesis of costly proteins. the opposite is true for plants which are relatively more k-strategists. if one considers that exposing functionally conserved sites to mutations is a strategy that favors saltatory evolution, and that small proteins have more regulatory and large ones more catalytic, basic maintenance functions, then this suggests that r-strategists maximize evolutionary potential for regulatory functions, but are relatively conservative at the level of house keeping functions. we find that the opposite is true for k-strategists. these results are in line with the view that the large numbers of offspring in r strategists enable many unsuccessful evolutionary "experiments" to be done at the level of regulatory functions, while this approach is not tenable for large plants with low offspring numbers. experimenting with large proteins is however more cost effective in k-strategists with large body sizes, because they have the reserves that enable them to cope with such costs, including the death of body parts, without killing the individual. the results on associations with domestication level, which suggest that the molecular evolutionary strategy of naturally k-selected plants becomes more r-like by domestication, presumably because domestication increases reproductive output, are more speculative. this is not only because our rankings of domestication levels between different species are difficult to justify in a formal way, but also because these analyses should be done separately on different groups of plants, such as eudicotyledons and others, or plants grown for their seeds versus those grown for other edible body parts. however, as the majority of domesticated plants are selected for increased reproductive output, the convergence between natural r-strategists and domestication in k-strategists seems plausible. it also strengthens the hypothesis that molecular adaptive evolutionary strategies associates with relative investments in reproduction.  high sequence variability, particularly at third codon positions, inversely correlates with length adjusted messenger rna stability, highlighting the functional and non-neutral aspects of synonymous site variation. this is also in line with synonymous codon usage bias which favors at-richness at third codon positions. in silico mutations of mrnas results in destabilization of longer mrnas but not of shorter mrnas, which could explain apparently contradictory results in change in folding stabilities upon randomization of rnas. our analyses also reveal transcript length as an important factor, which coupled with evolutionary variability controls secondary structure stability. this insight could substantially contribute to rna sequence and structure optimization studies that did not consider the effects of rna length. some highly conserved sites might be so not because of their functional importance in the protein, but as part of packing-protection constraints at the nucleic acid level. protection by secondary structure from mutations results in a pattern of adaptive evolution where the secondary structure of nucleic acid sequences prevents detrimental mutations, prior to or in addition to stabilizing selection on optimal secondary structures, a pattern that follows neo-lamarckian principles. in specific mrnas for some plants, functionally conserved sites are exposed to smutations, which might lead to saltatory adaptive evolution. protection from mutations at conserved sites is more widespread in r-strategists than k-strategists, perhaps because of greater mutation minimization pressures in smaller organisms. the above point is enhanced in r-strategists for large proteins, and for smaller proteins in k-strategists. this might reflect adaptive strategies, the first leading to saltatory evolution in r-strategists, the other to gradual adaptive specialization in k-strategists. among k-strategists, domestication causes them to resemble r-strategists in terms of molecular adaptive strategies. 
messenger rna secondary structure restricts this freedom, as revealed by previous findings in mitochondrial genes that mutations at third codon position nucleotides in helices are more selected against than those in loops. this motivated us to explore the constraints imposed by mrna secondary structure on evolutionary variability at all codon positions in general, in chloroplast systems. we found that the evolutionary variability and intrinsic secondary structure stability of these sequences share an inverse relationship. the positive correlation between secondary structure protection and functional importance of sites suggests that some sites might be conserved due to packing-protection constraints at the nucleic acid level in addition to protein level constraints. consequently, nucleic acid secondary structure a priori biases mutations. simulations of most likely single nucleotide evolution in psilotum nudum and nephroselmis olivacea mrnas, indicate that helix-forming propensities of mutated mrnas are greater than those of the natural mrnas for short sequences and vice-versa for long sequences. protection levels of functionally important sites vary across plants and proteins: r-strategists minimize mutation costs in large genes; k-strategists do the opposite. the differences between the protection levels of functionally important sites for r- and k-strategists reflect their respective molecular adaptive strategies. synonymous sites are freer to vary because of redundancy in genetic code. these converge with increasing domestication levels of k-strategists, perhaps because domestication increases reproductive output. moreover, helix-forming propensity estimated by the percentage of total mrna in helices increases gradually with mrna length, saturating beyond nucleotides. mrna length presumably predisposes shorter mrnas to evolve under different constraints than longer mrnas.
 synonymous sites are freer to vary because of redundancy in genetic code. messenger rna secondary structure restricts this freedom, as revealed by previous findings in mitochondrial genes that mutations at third codon position nucleotides in helices are more selected against than those in loops. this motivated us to explore the constraints imposed by mrna secondary structure on evolutionary variability at all codon positions in general, in chloroplast systems.  we found that the evolutionary variability and intrinsic secondary structure stability of these sequences share an inverse relationship. simulations of most likely single nucleotide evolution in psilotum nudum and nephroselmis olivacea mrnas, indicate that helix-forming propensities of mutated mrnas are greater than those of the natural mrnas for short sequences and vice-versa for long sequences. moreover, helix-forming propensity estimated by the percentage of total mrna in helices increases gradually with mrna length, saturating beyond nucleotides. protection levels of functionally important sites vary across plants and proteins: r-strategists minimize mutation costs in large genes; k-strategists do the opposite.  mrna length presumably predisposes shorter mrnas to evolve under different constraints than longer mrnas. the positive correlation between secondary structure protection and functional importance of sites suggests that some sites might be conserved due to packing-protection constraints at the nucleic acid level in addition to protein level constraints. consequently, nucleic acid secondary structure a priori biases mutations. the converse apparently occurs in a smaller number of cases, indicating a different evolutionary adaptive strategy in these plants. the differences between the protection levels of functionally important sites for r- and k-strategists reflect their respective molecular adaptive strategies. these converge with increasing domestication levels of k-strategists, perhaps because domestication increases reproductive output. 