Congruent evolution of genetic and environmental robustness in microRNA

0.46 MB·English

by  Gergely J. Szöllősi

Congruent evolution of genetic and environmental robustness in microRNA
Gergely J. Sz¨oll˝osi∗
E¨otv¨os University, Budapest, Hungary P´azm´any P´eter S´et´any 1/A.†
Imre Der´enyi‡
E¨otv¨os University, Budapest, Hungary P´azm´any P´eter S´et´any 1/A§
(Dated: January 19, 2009)
Genetic robustness, thepreservation of an optimal phenotypein theface of mutations, is critical
totheunderstandingofevolutionasphenotypicallyexpressedgeneticvariationisthefuelofnatural
selection. The origin of genetic robustness, whether it evolves directly by natural selection or it is
a correlated byproduct of other phenotypic traits, is, however, unresolved. Examining microRNA
9 (miRNA)genesofseveraleukaryoticspecies,BorensteinandRuppin(Borensteinetal.2006,PNAS
0 103:6593),showedthatthestructureofmiRNAprecursorstem-loopsexhibitssignificantlyincreased
0 mutationalrobustnessincomparisonwithasampleofrandomRNAsequenceswiththesamestem-
2 loop structure. The observed robustness was found to be uncorrelated with traditional measures
n of environmental robustness – implying that miRNA sequences show evidence of the direct evo-
a lution of genetic robustness. These findings are surprising as theoretical results indicate that the
J direct evolution of robustness requires high mutation rates and/or large effective population sizes
9 only found among RNA viruses, not multicellular eukaryotes. We demonstrate that the sampling
1 methodusedbyBorensteinandRuppinintroducedsignificantbiasthatleadtoanoverestimationof
robustness. Introducinganovelmeasureofenvironmentalrobustnessbasedontheequilibriumther-

] modynamic ensemble of secondary structures of the miRNA precursor sequences we demonstrate
E thatthebiophysicsofRNAfolding,inducesahighlevelofcorrelationbetweengenetic(mutational)
P and environmental (thermodynamic) robustness, as expected from thetheory of plastogenetic con-
. gruence introduced by Ancel and Fontana (Ancel et al. 2000, J. Exp. Zool. 288: 242). In light of
o
theoretical considerations we believe that this correlation strongly suggests that genetic robustness
i
b observed in miRNAsequences is thebyproduct of selection for environmental robustness.
-
q
[
Introduction [de Visser et al. 2003, Mayo and Bu¨rger 1997]: (i) the

2 most straightforward explanation, favoured by Wright,
v was that robustness evolves directly, through natural se-
8 Themagnitudeofgeneticeffectsonphenotypedepends lection [Fisher 1928]; (ii) an alternative congruent hy-
5 strongly on genetic background, the effects of the same potheses, put forward in the context of dominance by
6 mutation can be larger in one genetic background and Haldane, proposes that the evolution of genetic robust-
2
smaller in another. The idea that wild-type genotypes ness is a correlated byproduct of selection for environ-
.
0 are mutationally robust, i.e. show invariance in the face mental robustness, i.e. invariance in the face of nonher-
1 of mutations (more generally heritable perturbations), itable perturbations, e.g. temperature, salinity or inter-
8
goes back to Waddington [Waddington 1957], who origi- nal factors such as fluctuations in the concentration of
0
: nally introduced the concept as canalization. While ge- gene products during development [Haldane 1930]; (iii)
v netic robustness has been found across different levels of while a third view holds that genetic robustness is in-
Xi organization from individual genes, through simple ge- trinsic, arising simply because the buffering of a charac-
netic circuits to entire organisms (approximately 80% ter with respect to mutations is the necessary or likely
r
a of yeast single knockouts have no obvious effect in rich consequence of character adaptation, in the context of
medium [Hillenmeyer et al. 2008]), the origin of the ob- dominance Wright [Wright 1934] and later Kacser and
served robustness has remained a source of contention. Burns [Kacser and Burns 1981] argued that it arises as
Thethreemainhypothesesregardingthepotentialorigin an inevitable, passive consequence of enzyme biochem-
of genetic robustness predate the concept itself, and fall istry and selection for increased metabolic flux.
alongthelinesofthefamousdebatebetweenmembersof
Recently robustness has been a subject of re-
themodernsynthesis(inparticularWright,Haldaneand
newed interest. Several theoretical and simula-
Fisher)surroundingthe originof dominance (dominance
tion studies have addressed robustness in a wide
canbeunderstoodasasimplecaseofrobustness,adom-
range of contexts ranging from gene redundancy
inant phenotype being more robust against mutations)
[Krakauer and Plotkin 2002] to model regulatory net-
works [Siegal and Bergman 2002, Azevado et al. 2007,
Ciliberti et al. 2007a, Ciliberti et al. 2007b,
Crombach and Hogeweg 2008]. By in the first
∗Electronicaddress: [email protected]
case building on evidence of excess mutational
†URL:http://angel.elte.hu/∼ssolo
‡Electronicaddress: [email protected] robustness present in RNA secondary structure
§URL:http://angel.elte.hu/∼derenyi [Wagner and Stadler 1999] and in the second case
2
on the expectation that high mutation rates present ture they compared the single mutant neighborhood of
among RNA viruses should favour mutational robust- a given miRNA precursor sequence to the single mutant
ness [Wilke and Adami 2003] two pioneering studies, by neighborhood of the sample sequences. Calculating the
Montville et al. [Monetville et al. 2005] and Borenstein averagedistanceoftheMFEstructureofeachsinglemu-
and Ruppin [Borenstein and Ruppin 2006] have man- tant sequence to the MFE structure of the original se-
aged to step beyond computer simulations and through quencefor bothstem-loopandsamplesequences(details
using, respectively, in vitro evolution experiments on secondary structure calculations are presented be-
[Monetville et al. 2005] and microRNA sequences from low)theydemonstratedthatmiRNAprecursorsequences
diverse taxa found evidence to support the hypotheses have single mutant neighborhoods with sequences that
that genetic robustness can evolve directly. Further fold into more similar MFE structures compared to se-
work on in vitro evolution experiments has provided quences in the single mutant neighborhoods of sample
additional evidence showing that if a population is sequenceswithidenticalMFE structure. While asimilar
highly polymorphic robustness can evolve directly comparisonofthefoldingminimumfree-energyshoweda
[Sa´njuan et al. 2007, Bloom et al. 2007]. comparable,butlowerbias,thefindingthatthetwowere
only weakly correlated allowed the authors to conclude
The theoretical underpinnings of these studies is
that the observed bias is a result of direct selection for
provided by the results of van Nimwegen et al.
mutational robustness. Their results were reexamined
[van Nimwegen et al. 1999], who through solving the
by Shu et al. [Shu et al. 2007] who argued that muta-
quasispecies equations describing the evolution of a
tional robustness among miRNA precursors may be the
population on a network of phenotypically neutral se-
correlated byproduct of selection for environmental ro-
quences, were able to demonstrate, that provided a
bustness,butfoundonlyamoderatelyhighercorrelation
sufficiently polymorphic population, mutational robust-
using a different measure of mutational robustness.
ness can evolve directly. The necessary mutation rates
and/or population sizes were found to be very large in
In light of the consistently low value of uN among
e
simulation studies using RNA secondary structure as
multicellular eukaryotes the results of Borenstein and
a genotype-phenotype map [van Nimwegen et al. 1999,
Ruppin arehighly surprising. There is no knownmecha-
Forester et al. 2006, Sz¨oll˝osi and Der´enyi 2008], direct
nism which can explain the direct evolution of robust-
evolution of increased neutrality requiring the product
ness that they observe. According to the classic re-
ofthe effective populationsize N andthe mutationrate
e sults of Kimura and Maruyama the average fitness of
per nucleotide u to be well in excess of one. Such high
anasexuallyreproducingpopulation(inthelimitofvery
mutation rates can only readily be found among RNA
large populations) depends only on the mutation rate
viruses,are extraordinaryevenamong unicellular organ-
and is independent of the details of the fitness land-
isms (Prochlorococcus 2N u ≈ 2., E. coli 2N u ≈ 0.2,
e e scape [Kimura 1966]. This result, however, only holds
S. cerevisiae 4N u ≈ 0.09) and completely unheard
e under the assumption that the fittest genotype does not
of among multicellular eukaryotes possessing RNA si-
have any neutral sites. While, the extension of these
lencing mechanisms and microRNA genes (A. thaliana
resultstomoregeneralfitnesslandscapesbyvanNimwe-
4N u ≈ 0.012, D. melanogaster 4N u ≈ 0.015, C. ele-
e e gen et al. demonstrates that the presence of neutral
gans 4N u ≈ 0.013, C. intestinalis 4N u ≈ 0.012, M.
e e genotypes can lead to selective pressure to evolve mu-
musculus 4N u ≈ 0.001, H. sapiens 4N u ≈ 0.001)
e e tational robustness simulation studies using genotype-
[Lynch and Conery 2003].
phenotype maps induced by RNA secondary structure
In their study Borenstein and Ruppin examined mi- have demonstrated that uN > 1 is a necessary condi-
e
croRNA (miRNA) precursor sequences from several eu- tion [Forester et al. 2006] even in the presence of recom-
karyotic species. miRNA are small endogenous non- bination [Sz¨oll˝osi and Der´enyi 2008]. The case for the
coding RNAs that regulate the expression of protein direct evolution of genetic robustness rests on the paral-
coding genes through the RNA interference (RNAi) lelfindingsthatastrongerbiasformutationalrobustness
pathway [Lagos-Quintana et al. 2001, Lau et al. 2001, is present in miRNA precursor sequences than for envi-
Lee and Ambros 2001, Bartel 2004]. Functionally rele- ronmental robustness and that the two are only weakly
vant short (≈ 22nt) mature miRNA sequences are ex- correlated. Introducing a new measure of environmen-
cised from longer precursor sequences that fold into a tal robustness in this paper we endeavor to demonstrate
stem-loop hairpin structure. The hairpin like secondary that, indeed as previously also suggested by Shu et al.
structure of precursor stem-loops plays a crucial role in [Shu et al. 2007], the exact opposite is true: the bias
the maturationprocess[Bartel 2004] and is under evolu- for environmental robustness is stronger and it is highly
tionary constraint to conserve its structure. Borenstein correlated with mutational robustness. The correlated
and Ruppin used the novel and ingenuous method of evolutionofenvironmentalandmutationalrobustnessin
generating for each miRNA sequence a random sample RNAsequencesunderselectiontoretainsecondarystruc-
of sequences with identical minimum free-energy (MFE) ture is expected as a corollary of a general casual link
structure to uncover traces of adaptation. To compare between environmental robustness and genetic robust-
themutationalrobustnessofmiRNAprecursorsequences ness in RNA sequences proposed by Ancel and Fontana
to random sample sequences with identical MFE struc- [Ancel and Fontana 2000]. Theyarguethatplastogenetic
3
search using optimization only search using optimization only enna RNA package [Hofbacker et al. 1994] to produce a
a b 0.1 sequencewithMFEstructureidenticaltothatofthe na-
soepatrimchiz uastiionng 0.08 structuraln deiusttraanilctey bbaasseedd mmeeaassuurree tivesequencethatisstored(ii)andsubsequentlyrandom-
0.06
sequences with izing this sequence by attempting 4L random nucleotide
desired MFE structure nces 0.04 substitutionsinamiRNAprecursorsequenceoflengthL,
ue 0.02
rasnedqoumen scteart op seq 0 0 0.2 0.4 0.6 0.8 1 astcrcuecpttuinregraemsuabinstsiutuntcihoannigfetdh.eFroersuelatcinhgmsieRqNueAncper’escMurFsoEr
soepatrimcsheiz auarsticionhngr ua sn id+n o g m o i pzraatnitmdioodinmzesa siwtreieaoqdlukn Me +nFc Ee ss twruitchtu re fraction of stem-lo 000 ...0000.4681 strucsteuraarlnc dehiusr tutaraasnniidclnteyog bmb aoaisspzeeatdditm mimoieenzaaa sstuuiorreen + scafoeacrqlcuoMamellFnp4cEae3n6soyt1nriungacgevtneuoerrusaergawespsaa>ospcg8eiaer0nt0ceeodsranawmttaeiptdinhl.esS3st6ueh4qpe1uprelueonnmbcieuqessnutwtenaieslteshmsquiavdetaeenlnruciteeaiss-l
0.02 we considered.
random start
sequence 0
0 0.2 0.4 0.6 0.8 1
rank score
Measuring thermodynamic robustness
FIG. 1: a Generating a random sample of sequences
with a desired MFE structure by stochastic minimization
In order to calculate the thermodynamic robustness
of the free-energy of the desired fold (the method em-
measure η , defined below, we sampled the equilibrium
ployed the RNAinverse program used by Borenstein and t
Ruppin [Borenstein and Ruppin2006] as well as Shu et al. thermodynamic ensemble of stem-loop and sample se-
[Shu et al. 2007] ) results in a biased sample in which se- quences using the stochastic backtracking routine from
quences with lower than average neutrality (higher than av- the Vienna RNA package producing 106 suboptimal
erage number single mutant neighbors) are overrepresented. structures per sequence, using the default temperature
This can be avoided if after finding a sequence with the de- of 310 K. The average distance from the MFE struc-
sired MFE structure a random walk is performed among se- ture in the thermodynamic ensemble can be calculated
quences with the desired MFE structure. This random walk exactlywiththe helpofbase-pairingprobabilities,which
on the neutral network associated with the MFE structure
are available as a byproduct of partition function fold-
mimics the sequnce drift of a sequence evolving under the
ing inthe Vienna package,andwereused to validate the
constraint to fold in to the desired MFE structure. b Rank
sampling.
score distributions for two measures of mutational robust-
ness (ηs and ηn see text). Comparing the distributions de-
rived from sampling using only stochastic optimization (top,
r¯sbiased = 0.25,Rsbiased = 0.83,r¯nbiased = 0.37,Rnbiased = 0.66) Statistics
to that derived from sampling with subsequent randomiza-
tion (bottom, r¯s = 0.29,Rs = 0.78,r¯n = 0.44,Rn = 0.59) GivenarankscorerandsamplesizeN agoodestimate
shows that increased neutrality is predominately an artifact
for the probability of observing an equal or lower rank
of biased sampling, while the lower than average distance of
score by chance is given by (rN)/(N +1) ≈ r. Follow-
MFE structuresin themutational neighborhood to thewild-
ing[Borenstein and Ruppin 2006]rankscoresofr <0.05
typeMFEstructurebecomessomewhat lesspronounced,but
is still significant. are considered significantly robust. To determine if the
robustness of miRNA precursor sequences according to
some measure η has the same distribution as the robust-
congruence, i.e. the correlation between the set of struc- ness of sample sequences η′ for a group of sequences,
turesthermallyaccessibletoasequence,itsplastic reper- following [Borenstein and Ruppin 2006] we test against
toire, and the MFE structures of its genetic neighbour- the null hypothesis that they are drawn from identical
hood willleadtothe emergenceofmutationalrobustness distributions using the nonparametric Wilcoxon signed
inthepresenceofselectionforsomepredefinedstructure. rank test. In contrast to [Borenstein and Ruppin 2006],
however, we do not consider as paired values η and the
average of η′ over all N sample sequences, hη′i, as we
foundthistoresultinspuriouslylowp-values,butinstead
calculate the p-values for a given group of sequences by
Materials and Methods
averaging over 1000 different sets of {η,η′} pairs where
in each set the η′ values belong to a random sample se-
microRNA sequences and sampling
quence. As a complementary approach we also tested
the hypothesis that the distribution of rank scores of a
miRNA precursor sequences were downloaded from groupofsequencesforagivenrobustnessmeasureisuni-
miRBase version 9.0 [Griffiths-Jones et al. 2006]. All form–aswewouldexpectifmiRNAprecursorsequences
4361 miRNA genes were used, yielding 3641 unique were randomly sampled from the set of sequences with
miRNAprecursorsequences. ForeachmiRNAprecursor identicalMFEstructure–usingastandardKolmogorov-
sequence we produced a sample of random sequences by Smirnov goodness of fit test. We found the two signifi-
(i)usingthestochasticoptimizationroutinefromtheVi- canceanalysestobeingoodagreementindicatinghighly
4
significant bias for higher values of ηt(dth.) and ηs, but Thermodynamic ensemble Mutational neighborhood
mneslabdstoaeuameuoooqicpsWmsnksohustpstnineelolnsenysgseefagmepscqmsq3nseeusfeu6socoeioneR4iersfnetown1saNscnrsecileiudtaAlegosehintnnntvuwciphslfqetroyorielhuamtoerelhsenemcfltpnyuoamtasaavnrhtadotrsiisiemireRcaosioogdtcprsRneNinunatcsmeratAimseMafihrstqceiaceiuenuFpconanmlosgernEttndeenmastsccfcltseospoubeoiaqtnrgarsirnuasdnnuabidoesnayicyfinrrictmrfnccyouarcsaegueeronrssnaqetmdthoasc.utroiuepetouegmdirancTohatrtcnphenrhiuesrnoaaaaerasblgptleηmyaruienotisdonrpssh.feebtdlesccarneTu.ooarefshonnnooosters-----ff GMUGAGbCAfrequency in equilibrium ensembleFmCUUGUAEGaCdUUU AUCsoGUAAtGU-dUGrmAU uCA=0AUU GicAA.0rCA0 0UtG-.0UAu001CA .UAUUUr1111A GUAGeAUCCAUAGGGUAUAUCCAUA UAC0GUAUAGCGCGACAAUUAGCGUUCUfrequency in equilibrium ensembleAUGUAUGAUUAGCU A Ad0U5 U A11G.0CA0 U=eeA0.UG00 C--A.0U00000A3U.AU165111UAGUmUGAUArAGC CGaCU0iAA AUCURACU1GAnUAGACUUmAGAA0GCNdGCAUUAdoGACdrGUUCUoaAm UGi5UApsUn-GAUd mUAr GtdsCU1AeaA=UUaoA5iGCcrnAU1mAUm -GuC1AU4c1AU0ArUp UA eGUssUpGAlUA AeoGaCrfGC UCAA2A rUemrGUUUCAsAGA 0o CA1cACsUeAp5mueqlrequ sGMd uUAesGoCGC A=Ue2neU2FArGC0GnU25 UqcCEUsAUGc5UeAuUeG AUeUssAGeqCUAsAUA AnUtGAU2uGGCUGCArAAUUUGUA5CAAACAcUU AGeGAuACAU3UAeUUnU0AcCAUCsGCcUAtue 30rdGAAGUUUACeUUGGCUUUUG=ACGA UUAC3GUAUACUAG3A5GCUA UA63GCUA5UAUAAAUGUUAGAGCUUGGUAAAUAACGCUAACUCCAUA 4 400 frequency in mutational neighborhood dc 0 .00 .0G00UAG .CG1111CAUUAGCGUUfrequency in mutational neighborhoodCUAU GUA0UGAUUAGC UA0AdU UA11.0GC0 AUee0.=AU00 GC--.0AU00000A3U.AU651111UAGUU GA5UAAG CCG0CUAAAUCUACUGAUAGACUUAGAAGCGCAUUAGdCGUUCU AiU51GUAsUG0AUdUAtGCUaAA=UUAGCnAU1A UG1CcAU40AUAU eUAG1UUGAU 5AAGfCGCUrCAAAUGUUUCA oAGA1CAAC5UAm 2MGdU0AGC GC2A=UFUA0GCGU2UECUAUG5UAU GAUsUAG CU2A AUtAA2UGAUG5GrCUG5CAAAUUUGUACAAAACUUuAGGAACAUUAUUcUACAUCGCtUAu 3 0r3dGAAGUUe0UACUUGGCUUUUG=ACGAUUACGUAUACUAG3AGCUA UA36GC5UAUAUAA AUGUUAGAGCUUG3GUAAAUAACGCUAACUCC5AUA 40 40
ceptually similar to the approach used to support the distance from MFE structure distance from MFE structure
argument that the genetic code has evolved to minimize
mutational load [Di Giulio 1987, Haigh and Hurst 1991, FIG.2: Inordertoexamin therobustnessof miRNAprecur-
sorsequencestothermalfluctuationswesampledtheequilib-
Szathma´ry and Smith 1997]. In the case of the genetic
rium thermodynamic ensemble of structures. Sampling 106
codetheauthorstookthecommongeneticcode,and,for
structuresforeachmiRNAprecursorsequenceandeachmem-
each codon, calculated the change in polarity of the en-
beroftherandomsequencesample, webinnedstructuresac-
coded amino acid caused by replacing each of the three
cordingtotheirdistancefromtheMFEstructure. aFor,e.g.
nucleotides, one after the other. In order to determine
theMonodelphis domestica miRNA precursor sequence mdo-
whether the genetic code is adapted to minimize mu- mir-1 examining the distribution of structures as a function
tational load they proceeded by comparing the mean ofthebase-pairdistanceshowsthattheaveragedrandomse-
squaredchangecausedbythereplacementofasinglenu- quence sample distribution (white bars) has a much larger
cleotide in the common genetic code to 10000 randomly fraction of structures that are drastically different from the
generatedcodeswiththesameredundancies. Theyfound MFEstructure,comparedtothedistributionofstructuresfor
that only two of the random codes were more conser- the original miRNA precursor sequence (black bars). b Ex-
amining the averaged distribution of stem-loop (black bars)
vative than the common code with respect to polarity
and random corresponding random sequence sample distri-
distances between neighbouring amino acids.
butions (white bars) shows that there is a general tendency
WeundertookasimilarprograminthecaseofmiRNA
among miRNA precursor sequences for increased thermody-
precursor sequences. Each miRNA gene encodes a short
namic robustness, i.e. of avoiding structures that are highly
≈ 22 nucleotide sequence that is partially complemen- dissimilar to the MFE structure. A strikingly similar effect
tary to the mRNA of proteins regulated by the partic- canbeobservedifweexaminethedistributionofstructuresin
ular miRNA gene. For the proper short sequence to be themutationalneighborhood. Analogoustoa,incwebinned,
excised by the protein Dicer, and hence for the miRNA accordingtotheirdistancefromtheMFEstructureofthewild
gene to be functional, a larger part of the miRNA se- type,theMFEstructuresofall(3L)singlepointmutantsfor
quence, calledthe miRNA precursorsequence, mustfold theMonodelphis domestica miRNA precursor sequence mdo-
mir-1 (blackbars)aswellastheMFEstructuresforeachse-
into the proper secondary structure. In order to deter-
quenceinthesinglemutantneighborhoodofsamplesequences
minewhetheramiRNAprecursorsequenceisadaptedto
(whitebars). Thedistributionofstructuresinboththether-
minimize the effects ofmutationaland/orenvironmental
modynamic ensemble (a) and the mutational neighborhood
perturbations, i.e. to maximize mutational and/or envi-
(c) of the mdo-mir-1 miRNA precursor sequence have a sig-
ronmental robustness, we compared the mutational and nificantly smaller fraction of structures that are highly dis-
environmental robustness of each miRNA precursor se- similar then sample sequences with identical MFE structure.
quence (robustness measures we used are defined below) Comparing thetheaveraged distribution of stem-loop (black
tothemutationalandenvironmentalrobustnessofaran- bars)andrandomcorrespondingrandomsequencesampledis-
dom sample of sequences with identical structural phe- tributions (white bars) in the mutational neighborhood (d)
notype (i.e. identical MFE structure). to similar averaged distributions in the thermodynamic en-
sembles of the same sequences (b) shows that the tendency
To generate a random sample of sequences with given
among miRNA precursor sequences for increased robustness
MFE structure we first used, starting from a random
ispresentbothinthemutationalneighborhoodandthether-
sequence, stochastic minimization of the free-energy of
modynamic ensemble, i.e. miRNA precursor sequences show
the target structure to find a sequence with the de-
excessrobustnessintheface ofboththermalandmutational
sired MFE structure. This method by itself, however, perturbation.
yields a biased sample of sequences (see Fig. 1a and
b) and must be supplemented by an additional ran-
domization step (see Materials and Methods). To mea-
5
TABLE I: Phylogenetic breakdown of different measures of robustness
group / species r¯n Rn Sn r¯s Rs Ss r¯t(25) Rt(25) St(25) c(rs,rt(25)) # of seqs.
all 0.44 0.59 0.08 0.29 0.78 0.17 0.31 0.74 0.28 0.73 3641
vertebrate 0.46 0.55 0.06 0.31 0.78 0.13 0.29 0.75 0.30 0.76 2215
invertebrate 0.37 0.68 0.10 0.21 0.88 0.27 0.22 0.84 0.36 0.73 488
landplant 0.41 0.63 0.11 0.31 0.75 0.21 0.40 0.63 0.19 0.68 848
virus 0.38 0.66 0.09 0.23 0.84 0.18 0.21 0.85 0.32 0.65 82
Homo sapiens 0.48 0.53 0.04 0.32 0.75 0.11 0.28 0.76 0.33 0.74 471
Mus musculus 0.46 0.56 0.06 0.33 0.76 0.12 0.31 0.74 0.27 0.79 373
Drosophila melanogaster 0.40 0.64 0.08 0.22 0.88 0.24 0.23 0.82 0.35 0.74 78
Caenorhabditis elegans 0.30 0.78 0.18 0.20 0.89 0.37 0.23 0.82 0.34 0.75 114
Arabidopsis thaliana 0.39 0.62 0.11 0.29 0.78 0.19 0.43 0.60 0.15 0.75 131
Epstein-Barr virus 0.31 0.78 0.00 0.16 0.96 0.22 0.16 0.87 0.48 0.81 23
Average rank-scores that indicate significantly increased according to both measures discussed in theMaterials and Methods
section (p-value<10−3) are given in bold.
sure the mutational robustness of a given sequence we sequences with identical MFE structure we found that
used the measures introduced by Borenstein and Rup- miRNAprecursorsequences–incontrasttotheresultsof
pin[Borenstein and Ruppin 2006]: (i)the structuraldis- Borenstein and Ruppin – do not have significantly more
tancebasedmutationalrobustnessmeasureη ofanRNA neutralsinglemutantneighboursthansamplesequences,
s
sequenceoflengthLisdefinedbyη =1/(3L) 3L (L− but do show a statistically significant increase in robust-
s Pi=0
di)/L, where di is the base-pair distance between the ness measured according to ηs (see Fig. 1b and Table I).
secondarystructure of mutant i and the native sequence In other words, native miRNA precursor sequences have
(given by the number of base pairs present in one struc- onaveragethesamenumberofsinglemutantneighbours
ture,butnottheother),andthesumgoesoverall3Lsin- with MFE structures identical to their own, as random
gle mutant neighbours and (ii) the more stringent mea- sample sequences with the same structure. The MFE
sure η is simply defined as the fractionof neutralsingle structure of those single mutant neighbours that are not
n
mutant neighbours, i.e. those that have identical MFE identicaltotheirownare,ontheotherhand,significantly
structure to the original sequence. In order to quantify more similar than in the case of sample sequences.
the level of excess mutational robustness among miRNA
The presence of excess mutational robustness is, by
precursor sequences we counted, for each miRNA pre-
itself, insufficient to determine whether mutational ro-
cursor sequence, the number of sample sequences that
bustness has evolved as a result of direct selection or in
have higher mutational robustness according to a given
congruence with selection for environmental robustness.
measure (see Materials and Methods) and used this to
As established previously [Borenstein and Ruppin 2006]
calculate the rank scores r and r , defined as the frac-
s n there is evidence for excess thermodynamic robustness,
tion ofsample sequences with identicalor higher robust-
robustness to thermal fluctuations, as evidenced by a
ness according to η and η , respectively. To facilitate
s n significantly lower than chance minimum folding energy
an overview of the extent of excess mutational robust-
among miRNA precursor sequences. Defining the envi-
ness we also calculated the average of the rank scores
ronmental robustness measure η simply as minus the
E
overall miRNA precursorsequences r¯ and r¯ as well as
s n minimum folding energy we also find r¯ = 0.278, R =
E E
the fraction of miRNA precursor sequences with higher
0.796, S = 0.220 using unbiased sampling. The cor-
E
than average robustness (i.e. rank-scores < 0.5) R and
s relation between r and r across miRNA precursor se-
s E
R and the fraction of sequences with statistically sig-
n quences is, however,rather low with a Pearson’scorrela-
nificantincreasedrobustness(i.e. rank-scores<0.05,see
tioncoefficientofc(r ,r )=0.217andc(r ,r )=0.071.
s E n E
MaterialsandMethods) S andS , respectively,accord-
s n The minimum folding energy is a somewhat crude mea-
ingtoagivemeasure. Thestatisticalsignificanceofboth
sure of thermodynamic robustnessand does not evenre-
rankscoresforindividualmiRNAprecursorsequencesas
flect the excess mutational robustness according to the
well as that of the finding a given fraction of robust se-
measure η . There is no good reason to assume that
s
quences among a group of sequences was determined as
a low MFE in itself confers environmentalrobustness,as
detailed in the Materials and Methods section.
evenforrelativelyhighfreeenergiesagivensequencemay
ReexaminingthemutationalrobustnessofmiRNApre- none the less with high probability fold into structures
cursorsequencesincomparisontoanunbiased sampleof sufficiently similarto the MFE structureto remainfunc-
6
To construct an appropriate measure of thermody-
fractionrobustR fractionresolved
t
a 0 .18 0. 91 namic robustness that also reflects this observation we
00..46 00..78 would need to know the extent of similarity that is re-
0.2 0.6
0 0.5 quired to retain functionality – indeed we would require
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
threshold threshold detailed knowledge of the interaction between the RNA
averagerankscorer¯t correlationc(rt,rs) substrateandtheenzymeDicertoestablishanappropri-
0.5 1
00..34 00..68 ate measure of structure similarity. As such information
0.2 0.4
0 .01 0. 20 isnotatpresentavailablewechosetousethemostsimple
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
threshold threshold and widely employed structure similarity measure, the
base-pairdistance usedabove. Inorderto determine the
b 0.4
extentofsimilarityrequiredtoretainfunctionalitywede-
0.35
fined the threshold thermodynamic robustness measure
quences 0 .002..523 dds tt t hhr ..u c== t u 32 r 05 a l dtthhisrrteeassnhhcooell ddb atthhseeerrdmm moodudytyanntaiaommniiaccl rrrooobbbuuussstttnnneeessssss ηeqt(udatthi.n)gaistawfituhnctthioenproofbtahbeilitthyreinshtohlde edqisutialinbcreiudmtht.,hebry-
of stem-loop se 00 ..010.581~~ c(rs,rt(P3e0a)r)s=on0's. 7c1orarneadltcio(rns c,orte(f2f.5: ))=0.73 smdpiosetdcatynntcoaemstheicequeManlFsetEmosbotlrreulecostfsusrttehr,uaicn.et.uaretshrtehsahtolhdadvteh.bwasiteh-praeir-
n 0.06
fractio 00..0024 ηt(dth.)=XH(dth.−di) e−EZi/kT, (1)
i∈Ω
0
0 0.1 0.4 0.6 0.8 1
rank score wherethesumgoesoverthesetofallpossiblestructures
Ω, d denotes the base-pair distance of structure i to the
i
FIG.3: Toquantifytheextentofthermodynamicrobustness MFE structure, Z = e−Ei/kT is the partition sum
Pi∈Ω
presentinmiRNAprecursorsequencesweexaminedtherank andH(x)istheunitstepfunction,i.e.H(x)=0ifx<0
statistics ofthethresholdthermodynamicrobustnessηt(dth.) and H(x)=1 if x≥0.
whichwedefinedasthecumulativefrequencyofstructuresin
Examining the thermodynamic robustness of miRNA
the thermodynamic ensemble that are equal to or less than
precursor sequences in comparison to an unbiased sam-
a distance threshold dth.. a For threshold values lager than
dth. ≈ 20 the average rank of miRNA precursor sequences ple of sequences with identical MFE structure we found
with respect to ηt(dth.), denoted by r¯t(dth.), becomes larger thatmiRNAprecursorsequenceshavesignificantlymore
thanr¯s,whileremaininghighlycorrelatedwithit. Forthresh- structures in their equilibrium thermodynamic ensemble
olds dth. > 25 the ηt(dth.) values for an increasingly larger that are similar to the MFE structure than sample se-
fraction of random sample sequences becomes indistinguish- quences (see Fig. 2a,b and Table I). In other words
able from the ηt(dth.) valueof the original miRNA precursor miRNA precursor sequences tend to adapt more similar
sequenceduetoalackofstructureswithd>dth. inourfinite structures as a result of thermal fluctuations than ran-
resolution sample of the equilibrium ensemble. The fraction
dom sample sequences with the same structure. Calcu-
of sequences which remain distinguishable are labeled frac-
tion resolved. b Among the majority of the random sample lating the average rank score r¯t(dth.) and the fraction of
sequences that are resolved, however, the threshold robust- robust Rt(dth.) and significantly robust St(dth.) miRNA
ness of large number of the miRNA precursor sequences be- precursorsequences,withrespecttothemeasureηt(dth.)
(Fig. 3a)andexaminingthedistributionofstructuresas
comes highly significant (black and white bars) compared to
themutationrobustnessasmeasuredbyηs(greybars),whilst afunctionofthebase-pairdistanceforindividualmiRNA
retaining a high correlation among rs and rt(dth.). precursorsequences(seee.g.Fig.2a)indicatesthatabove
a threshold distance dth. ≈20 the measures start to sat-
urate, yielding an estimate of the required similarity to
retain function. The correlation between the rank-score
ofmiRNA precursorsequences accordingto the distance
tional. The largenumber ofmiRNA precursorsequences similarity based mutational robustness measure and the
that exhibit excess mutational robustness as measured thresholdthermodynamicsmeasureishighforallthresh-
by the structural similarity based measure η suggests old values. This is the direct result of the high degree
s
that a strict adherence to the MFE structure is not nec- of similarity between the distribution of structures in
essary to retain functionality – a sufficiently similar, but the thermodynamic ensemble and the mutational neigh-
not necessarily identical, secondary structure is enough borhood (Fig. 2). The average rank score r¯t(dth.) and
toguaranteetheexcisionofthepropersubsequence. This the fraction of robust Rt(dth.) and significantly robust
is further supported by the fact that folding free energy St(dth.) miRNA precursor sequences with respect to the
aloneisnotsufficienttodiscriminatemiRNAprecursors, threshold thermodynamic robustness measure indicate a
as well as recentevidence that a diverse set ofstructural markedly larger extent of excess robustness than their
features are needed for successful cleavage of a miRNA counterparts for mutation robustness, i.e. r¯ , R and S
s s s
precursor sequence [Ritchie et al. 2007]. (see Fig. 3a,b and Table I) above dth. >20.
7
Discussion evolution of environmental and genetic robustness.
The correlation between the response to heritable
The results presented above demonstrate the corre- (mutational) and nonheritable (thermodynamic) pertur-
latedpresenceofexcessenvironmental(thermodynamic) bation, and hence the congruent evolution of genetic
and genetic (mutational) robustness amongmiRNA pre- and environmental robustness may extend to other sys-
cursor sequences as measured according to, respectively, temswithgenotype-phenotypemapsdifferentfromRNA
ηs andηt(dth.). A rathergeneralcausality betweenenvi- secondary structure. In particular, Xia and Levitt
ronmentaland genetic robustnessin the contextof RNA [Xia and Levitt 2002],havefoundcompellingevidenceof
secondary structure has been suggested by Ancel and thecorrelatedevolutionofincreasedthermodynamicsta-
Fontana [Ancel and Fontana 2000], who studied the dy- bilityandthenumberofneutralneighboursinlatticepro-
namicsofaninsilico populationofRNAsequencesevolv- tein models. Understanding the relationshipbetween se-
ingtowardsapredefinedtargetshape. Theyfoundthata quence, structure, and function is, and will remainto be
correlationexistsbetweenthesetofshapesinthe plastic in the foreseeable future, a central theme in both molec-
repertoire of a sequence and the set of dominant (min- ular and evolutionary biology. A comprehensive view of
imum free energy) shapes in its genetic neighborhood. how the relationship between sequence, structure, and
They argue that this statistical property of the RNA function is shaped during the course of evolution must
genotype-phenotype map, which they call plastogenetic take into consideration both the potential correlations
congruence, traps populations in regions where most ge- that arise from the physics of the structure-sequence re-
netic variation is phenotypically neutral. In other words lationshipas wellas the relevantpopulationgenetic con-
RNA sequences explore a similar repertoire of subopti- ditions in the context of which it takes on the role of
mal structures as a result of perturbations due to muta- a genotype-phenotype map. In the context of compu-
tions and perturbations resulting from thermal fluctua- tational miRNA gene discovery our thermodynamic ro-
tions,andselectionforagiventargetstructurefavoursse- bustness measure potentially offers an improved struc-
quences with higher robustness to perturbations of both tural feature that may perform better than the free en-
type. ergyscoreof the hairpinorits ensemble diversity(which
Since, in contrast to genetic robustness, environmen- have proved uninformative [Freyhult et al. 2005]).
tal robustness does not require high values of uN , as Acknowledgments
e
it is a property of the sequence and not its mutational
neighborhood, we contend that the observedbias in mu- This work was partially supported by the Hungarian
tational robustness is in fact the result of the congruent Science Foundation (K60665).
[Ancel and Fontana 2000] Ancel LW, Fontana W (2000) [Ciliberti et al. 2007b] Ciliberti S, Martin OC, Wagner A
Plasticity, Evolvability, and Modularity in RNA J. Exp. (2007) Innovation and robustness in complex regulatory
Zool. 288: 242-283. gene networks. Proc. Natl. Acad. Sci. U.S.A. 104: 13591-
(DOI:10.1002/1097-010X(20001015)288:3¡242::AID- 13596.
JEZ5¿3.0.CO;2-O) (DOI:10.1073/pnas.0705396104)
[Azevado et al. 2007] Azevedoetal. (2007)Sexualreproduc- [Crombach and Hogeweg 2008] Crombach A, Hogeweg P
tion selects for robustness and negative epistasis in artifi- (2008) Evolution ofEvolvability in Gene Regulatory Net-
cial gene networks. Nature 440: 87-90. works PLoS Comput. Biol. 4: e1000112.
(DOI:10.1038/nature04488) (DOI:10.1371/journal.pcbi.1000112)
[Bartel 2004] Bartel DP (2004) MicroRNAs: genomics, bio- [Di Giulio 1987] Di Giulio M (1987) The extension reached
genesis, mechanism, and function. Cell 116: 281–297. by the minimization of the polarity distances during the
(DOI:doi:10.1016/S0092-8674(04)00045-5) evolution of the genetic code. J. Mol. Evol. 29: 288-293.
[Bloom et al. 2007] Bloom JD,LuZ, ChenD,RavalA,Ven- (DOI:10.1007/BF02103616)
turelliOS,ArnoldFH(2007)Evolutionfavorsproteinmu- [Fisher 1928] FisherRA(1928) Thepossiblemodificationsof
tationalrobustnessinsufficientlylargepopulations. BMC theresponseofthewildtypetorecurrentmutations.Am.
Biol. 5: 29. Nat. 62:115-126.
(DOI:10.1186/1741-7007-5-29) (DOI:10.1086/280193)
[Borenstein and Ruppin2006] Borenstein E, Ruppin E [Forester et al. 2006] Forster R, Adami C, Wilke CO (2006)
(2006) Direct evolution of genetic robustness in mi- Selection for mutational robustness in finite populations
croRNA Proc. Natl. Acad. Sci. U.S.A.103: 6593-6598. J. Theor. Biol. 243: 181-190.
(DOI:10.1073/pnas.0510600103) (DOI:10.1016/j.jtbi.2006.06.020)
[Ciliberti et al. 2007a] Ciliberti S, Martin OC, Wagner A [Freyhult et al. 2005] Freyhult E, Gardner PP, Moulton V
(2007) Robustness can evolve gradually in complex reg- (2005) A comparison of RNA folding measures. BMC
ulatory networks with varying topology. PLoS Computa- Bioinformatics 6: 241.
tional Biology 3: e15. (DOI:10.1186/1471-2105-6-241)
(DOI:10.1371/journal.pcbi.0030015) [Griffiths-Jones et al. 2006] Griffiths-Jones S, Grocock RJ,
8
van Dongen S, Bateman A, Enright AJ (2006) miRBase: [van Nimwegen et al. 1999] van Nimwegen E, Crutchfield P,
microRNAsequences,targetsandgenenomenclature.Nu- HuynenM(1999)Neutralevolutionofmutationalrobust-
cleic Acids Research 34: D140. ness Proc. Natl. Acad. Sci. U.S.A.96: 9716-9720.
(DOI:10.1093/nar/gkj112) (DOI:n.a. )
[Haigh and Hurst 1991] Haig D, Hurst LD (1991) A quanti- [Ritchie et al. 2007] Ritchie W, Legendre M, Gautheret D
tative measure of error minimization in the genetic code. (2007) RNA stem-loops: to be or not to be cleaved by
J. Mol. Evol. 33: 412-417. RNAseIII.RNA 13: 457-462.
(DOI:10.1007/BF02103132) (DOI:10.1261/rna.366507)
[Haldane 1930] Haldane, JBS (1930) A note on Fisher’s the- [S´anjuan et al. 2007] Sanju´an R, Cuevas JM, Furi´o V,
ory of dominance. Am. Nat. 64:87-90. Holmes EC, Moya A (2007) Selection for robustness in
(DOI:10.1086/280299) mutagenized RNAviruses. PLoS Genet. 3: e93.
[Hillenmeyer et al. 2008] Hillenmeyer ME et al. (2008) The (DOI:10.1371/journal.pgen.0030093.eor)
Chemical Genomic Portrait of Yeast: Uncovering a Phe- [Shu et al. 2007] ShuWetal.(2007) Insilicogeneticrobust-
notypefor AllGenes Science 320: 362. nessanalysisofmicroRNAsecondarystructures: potential
(DOI:10.1126/science.1150021) evidenceofcongruentevolutioninmicroRNA.BMC Evo-
[Hofbacker et al. 1994] HofackerIL,etal.(1994)FastFolding lutionary Biology 7: 223.
and Comparison of RNA Secondary Structures. Monat- (DOI:10.1186/1471-2148-7-223)
shefte fu¨r Chemie 125: 167-188. [Siegal and Bergman 2002] Siegal LS, Bergman A (2002)
(DOI:doi:10.1007/BF00818163) Waddington’s canalization revisited: Developmental sta-
[Kacser and Burns 1981] Kacser H, Burns JA (1981) The bilityandevolutionProc.Natl.Acad.Sci.USA99: 10528.
molecular basis of dominance. Genetics 97:6639-6666. (DOI:10.1073/pnas.102303999)
(DOI:n.a. ) [Szathm´ary and Smith 1997] Szathm´aryE,SmithJM(1997)
[Kimura 1966] Kimura M,Maruyama T (1966) The muta- The major transitions in evolution (Oxford University
tional load with epistatic gene interactions in fitness. Ge- Press, Oxford, UK)
netics 54: 1337-1351. (ISBN:019850294X)
(DOI:n.a. ) [Sz¨oll˝osi and Der´enyi2008] Sz¨oll˝osiGJ,Der´enyiI(2008)The
[Krakauer and Plotkin 2002] Krakauer DC, Plotkin JB effectofrecombinationontheneutralevolutionofgenetic
(2002) Redundancy, antiredundancy, and the robustness robustness. Math. Biosci. 214: 58-62.
of genomes. Proc. Natl. Acad. Sci. USA99: 1405. (DOI:10.1016/j.mbs.2008.03.010)
(DOI:10.1073/pnas.032668599) [de Visser et al. 2003] de Visser JAG et al. (2003) Perspec-
[Lagos-Quintana et al. 2001] Lagos-Quintana M, Rauhut R, tive: EvolutionandDetectionofGeneticRobustnessEvo-
LendeckelW,TuschlT(2001)Identificationofnovelgenes lution 57: 1959-1972
coding for small expressed RNAs.Science 294:853-858. (DOI:10.1554/02-750R)
(DOI:10.1126/science.1064921) [Wright 1934] WrightS(1934)Physiologicalandevolutionary
[Lau et al. 2001] LauNC,LimLP,WeinsteinEG,BartelDP theories of dominance. Am. Nat. 68:25-53.
(2001)AnabundantclassoftinyRNAswithprobablereg- (DOI:10.1086/280521)
ulatory roles in Caenorhabditis elegans. Science 294:858- [Waddington 1957] Waddington, CH, Kacser H (1957) The
862. Strategy of the Genes: A Discussion of Some Aspects of
(DOI:10.1126/science.1065062) Theoretical Biology (MacMillan, NewYork USA).
[Lee and Ambros2001] LeeRC,AmbrosV(2001) Anexten- (ISBN:n.a. )
sive class of small RNAs in Caenorhabditis elegans. Sci- [Wagner and Stadler 1999] Wagner A, Stadler PF (1999)
ence 294:862-864. Viral RNA and evolved mutational robustness. J Exp
(DOI:10.1126/science.1065329) Zool. 285:119-27
[Lynch and Conery 2003] Lynch M, Conery JS (2003) The (DOI:10.1002/(SICI)1097-010X(19990815)285:2¡119::AID-
Origins of Genome Complexity Science 302: 1401 - 1404. JEZ4¿3.0.CO;2-D )
(DOI:10.1126/science.1089370) [Wilke and Adami 2003] Wilke CO, Adami C (2003) Evolu-
[Mayo and Bu¨rger 1997] MayoO,Bu¨rgerR(1997)Evolution tion of mutational robustness. Mutat. Res. 522: 3-11
ofdominance: atheorywhosetimehaspassed? Biol.Rev. (DOI:10.1016/S0027-5107(02)00307-X)
72: 97-110. [Xia and Levitt 2002] XiaY,LevittM(2002) Rolesofmuta-
(DOI:10.1111/j.1469-185X.1997.tb00011.x) tionandrecombinationintheevolutionofproteinthermo-
[Monetville et al. 2005] MontvilleRetal.(2005)Evolutionof dynamics Proc. Natl. Acad. Sci. U.S.A.99:10382-10387.
Mutational Robustness in an RNA Virus PLOS Biology (DOI:10.1073/pnas.162097799)
3:e381
(DOI:10.1371/journal.pbio.0030381)