Dynamic RNA Fitness Landscapes of a Group I Ribozyme during Changes to the Experimental Environment

1.1 MB·English

by  Gianluca Peri, Clémentine Gibard, Nicholas H Shults, Kent Crossin, Eric J Hayden

Dynamic RNA Fitness Landscapes of a Group I Ribozyme during
Changes to the Experimental Environment
Gianluca Peri,*,1 Cl(cid:2)ementine Gibard,†,2 Nicholas H. Shults,†,2 Kent Crossin,2 and Eric J. Hayden*,1,2
1BiomolecularSciencesGraduatePrograms,BoiseStateUniversity,Boise,ID,USA
2DepartmentofBiologicalScience,BoiseStateUniversity,Boise,ID,USA
†Presentaddress:FACibleBiodiagnostics,366SW5thAve.,Suite104,Meridian,ID,USA
*Correspondingauthors:E-mails:[email protected];[email protected].
Associateeditor:JianzhiZhang
Abstract
Fitness landscapes of protein and RNA molecules can be studied experimentally using high-throughput techniques to D
o
measure the functional effects of numerous combinations of mutations. The rugged topography of these molecular wn
fitnesslandscapesisimportantforunderstandingandpredictingnaturalandexperimentalevolution.Mutationaleffects loa
d
are also dependent upon environmental conditions, but the effects of environmental changes on fitness landscapes ed
remainspoorlyunderstood.Here,weinvestigatethechangestothefitnesslandscapeofacatalyticRNAmoleculewhile fro
m
changingasingleenvironmentalvariablethatiscriticalforRNAstructureandfunction.Usinghigh-throughputsequenc- h
ing of in vitro selections, we mapped a fitness landscape of the Azoarcus group I ribozyme under eight different ttp
s
concentrations of magnesium ions (1–48mM MgCl2). The data revealed the magnesium dependence of 16,384 muta- ://a
c
tional neighbors, and from this, we investigated the magnesium induced changes to the topography of the fitness a
d
landscape. The results showed that increasing magnesium concentration improved the relative fitness of sequences em
at higher mutational distances while also reducing the ruggedness of the mutational trajectories on the landscape. ic.o
u
As a result, as magnesium concentration was increased, simulated populations evolved toward higher fitness faster. p
.c
Curve-fitting of the magnesium dependence of individual ribozymes demonstrated that deep sequencing of in vitro o
m
reactions can be used to evaluate the structural stability of thousands of sequences in parallel. Overall, the results /m
b
highlight how environmental changes that stabilize structures can also alter the ruggedness of fitness landscapes and e
/a
alterevolutionaryprocesses. rtic
le
Keywords:fitnesslandscape,gene–environmentinteractions,molecularevolution,ribozyme,noncodingRNA. /3
9
/3
/m
s
a
b
Introduction 3
helped explain the evolutionary constraints to optimizing 7
3
Fitness landscapes are maps of how changes to biological gene functions, which have important implications for un- A/65
r0
sequences(genotypes)resultinphenotypic orfitnessdiffer- derstandingandpredictingevolutionaryoutcomes(Wagner t22
ences at the genome or individual gene level. The fitness 2011; Ferretti et al. 2018). The rich data sets are providing ic89 b
landscapes of entire genomes are too vast to explore com- opportunities to interpret past evolution (Yang et al. 2019; ly
e g
prehensively,andexperimentalapproacheshavefocusedon Bendixsen et al. 2021), predict future evolution (de Visser ue
s
theinvestigationofsmallerlocalfitnesslandscapes(deVisser andKrug2014),andmoreeffectivelyharnessmolecularevo- t o
n
and Krug 2014). The genotypes in these local fitness land- lutioninthelaboratory(Ogdenetal.2019;Townshendetal. 0
5
scapesarechosenfromcombinationsofmutationsthathave 2021). A
p
beenidentifiedasimportantbasedontheoutcomesofnat- Epistasisisanimportantfeatureofmolecularfitnessland- ril 2
ural or experimental evolution, or from bioinformatics scapes.Epistasisinthiscontextreferstoanonadditiveeffect 0
2
approaches. Molecular fitness landscapes result from the when two or more mutations are combined, such that the 2
studyofnumerouscombinationsofmutationswithinasingle fitnesseffectofamutationdependsonthepresenceorab-
geneproduct,suchasaproteinoranoncodingRNAmole- senceofothermutations.Inextremecases(signepistasis)the
cule. For example, the fitness landscape of beta lactamase effectofamutationcanbebeneficialordeleteriousdepend-
proteins has been used to study the evolution of antibiotic ingonthepresenceorabsenceofothermutations.Interms
resistance(Weinreichetal.2006),andRNAfitnesslandscapes ofthetopographyoffitnesslandscapes,signepistasiscauses
havebeenusedtostudytheemergenceofnewRNAstruc- ruggedness—multiple peaks and valleys in the landscape—
tures (Jim(cid:2)enez et al. 2013; Bendixsen et al. 2019; Pressman where peaks are commonly defined as genotypes where
et al. 2019). The study of molecular fitness landscapes has all the analyzed mutational neighbors have lower fitness
(cid:2)TheAuthor(s)2022.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionLicense(https://creativecommons.
org/licenses/by/4.0/),whichpermitsunrestrictedreuse,distribution,andreproductioninanymedium,providedtheoriginalworkis Open Access
properlycited.
Mol.Biol.Evol.39(3):msab373 doi:10.1093/molbev/msab373 AdvanceAccesspublicationJanuary10,2022 1
MBE
Perietal. . doi:10.1093/molbev/msab373
(Poelwijk et al. 2011). Landscape ruggedness is thought to selection can be detected by high-throughput sequencing
constrain evolution by making some sequential orders of (Hayden et al. 2011, 2015). This experimental system allows
mutationsunlikelytooccur,evenifthefinalmutationalcom- the study of numerous genetic differences while systemati-
bination is beneficial (Ferretti et al. 2018). To evolve away callyalteringasingleisolatedenvironmentalvariable.
from a peak genotype toward higher fitness may require Thecombinationsofmutationsstudiedherewerechosen
thepreservationofadeleteriousmutationinthepopulation based on mutational effects observed in a previous study
longenoughforasecondmutationtoimprovefunction,or where random mutations were introduced into the
thesimultaneousoccurrenceoftwoormorespecificmuta- Azoarcus ribozyme followed by in vitro selection for
tions. If a genotype is a peak in the local fitness landscape, reverse-splicing activity under low (2mM) or high (25mM)
randommutationsareunlikelytoimprovethefunction,and MgCl concentrations (Hayden et al. 2015). The selection
2
thereforeevolutionislesslikelytopreservethemultiplemu- enriched for ribozymes that were active under each condi-
tationalchangesrequiredtogenerateahigh-fitnessgenotype. tion.Followingthisselection,severalsequenceswereisolated
Numerous experimental data sets have demonstrated that andstudiedfortheirsensitivitytoMgCl bymeasuringtheir
2
D
molecularfitnesslandscapesareinfactruggedwithmultiple activityrelativetothewild-typeribozymeatboth2mMand o
w
suboptimal peaks caused by the epistatic interactions of 25mMMgCl .ThisanalysiscanbefoundintheSupporting n
2 lo
mutations(Szendroetal.2013). Informationofthepreviouspublication(Haydenetal.2015). a
d
e
Anotherimportantaspectoffitnesslandscapesisthein- Toidentifythenucleotidepositionsstudiedherein,wechose d
teractionbetweengeneticdifferencesandtheenvironment. from this previous analysis one sequence with high Mg2þ fro
m
ProteinandRNAfunctionsaresensitivetonumerouscom- sensitivity,onesequencewithlowMg2þ sensitivityandone h
ponentsofthecellular environment,suchasionic strength, sequencewithintermediateMg2þsensitivityfromtheprevi- ttps
pH, temperature, the expression level of other genes, and ous analysis. The high sensitivity sequence (E27-1) showed ://a
c
metabolite concentrations. Experimental investigations of very low activity at 2mM MgCl2 and high activity at ade
RNAandproteinlandscapesinyeasthaveshownthatalter- 25mM MgCl and contained the mutations C90G, A128C, m
2 ic
ingthegrowthconditionscanchangeepistaticinteractions, andC177U.ThesequencewithlowMgCl2sensitivity(D45-5) .ou
andthetopographyofthemolecularfitnesslandscape(Liand showed high activity at both 2mM and 25mM MgCl and p
2 .c
Zhang 2018; Flynn et al. 2020). The constant fluctuation in contained the mutations A153G and U184C. The third se- om
natural environments means that fitness landscapes are in quence (E27-12) showed intermediate activity at 2mM Mg /m
b
factdynamic.Measuringchangestotheeffectsofnumerous andhighactivityat25mMMgandcontainedthemutations e/a
mutationsindifferentenvironmentsremainsachallengefor A119G and C158U. To make the library used here, we ran- rtic
characterizing fitness landscapes and for using the informa- domizedallsevennucleotidepositions(90,119,128,153,158, le/3
9
tionforvariousapplications. 177,184)wheremutationswerefoundinthesethreesequen- /3
Here,westudyfitnesslandscapesofanoncodingRNAina ces(fig.1A).Thisresultedinalibraryof47¼16,384different /m
s
a
cell-freebiochemicalreactionwhileexperimentallychanging sequences(fourpossiblenucleotidesatsevenpositions)and b
3
themagnesiumionconcentration.RNAfoldingandfunction included numerous mutations of unknown function. We 73
arestronglydependentupontheconcentrationofionsinthe expectedthatthemutationswouldshowarangeofmagne- /65
0
environment.Magnesiumisparticularlyeffectiveatstabilizing siumdependenceandprovideopportunitiesformagnesium 22
8
thenativefoldsofRNAmolecules,andtheamountofmag- inducedchangestothetopographyofthefitnesslandscape. 9
b
nesiumrequiredisrelatedtothethermodynamicstabilityof Sevenrandomnucleotidepositionswerechosentoachievea y g
the native state (Draper 2008). Group I self-splicing introns numberofsequencevariantsthatcouldberealisticallystud- ue
s
havebeenafoundationalmodelofRNAstructureandfolding ied under multiple magnesium concentrations, in triplicate, t o
n
inpartbecausethecatalytic activity oftheseribozymesis a and still have full sequencing coverage at adequate depth 0
5
very sensitive readout for native tertiary structure (Rangan usingcurrenthigh-throughputsequencingplatforms. A
p
etal.2004;Sinanetal.2011).Forourexperiments,wechoseto ril 2
study mutational variants of the group I self-splicing intron 0
Results 22
from the bacterium Azoarcus (Tanner and Cech 1996;
GenBank accession number DQ103524.1). This ribozyme is Todeterminethefitnessofall16,384RNAsequencesineight
well characterized in terms of reaction kinetics (Gleitsman different magnesium concentrations, the RNA pool was
andHerschlag2014),folding(Ranganetal.2004),andstruc- allowed to react with a reverse-splicing substrate under a
ture(StahleyandStrobel2005),providingguidanceforreac- givenmagnesiumconcentration,andthereactedribozymes
tion conditions and for interpreting results. In addition, the were amplified by reverse-transcription PCR (RT-PCR). The
Azoarcusribozymehasbeenengineeredtostudytheoriginof frequency ofeach sequenceinthepoolsafterselection was
life(Vaidyaetal.2012),andtodelivermRNArepairinggene compared with their frequency before selection as deter-
therapies(DolanandMu¨ller2014).Importantly,thereverse- mined by high-throughput sequencing (see Materials and
splicingreactionofthisribozymeallowsfortheenrichmentof Methods). High-throughput sequencing yielded sufficient
molecular variants based on how much reacted during a readcoverageforeachgenotypeinthepools(supplementary
biochemicalassay(BeaudryandJoyce1992).Numerousribo- fig. 2, Supplementary Material online). We define ribozyme
zymesequencevariantscanbecosynthesizedandstudiedin fitnessbasedonthechangeinpopulationfrequencyofeach
parallel. The enrichment of each sequence from in vitro sequence because the relative activity of each ribozyme
2
MBE
DynamicRNAFitnessLandscapesofaGroupIRibozyme . doi:10.1093/molbev/msab373
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
/3
9
/3
/m
s
a
b
FIG.1.Invitroselectionofribozymeactivityunderdifferentmagnesiumconditions.(A)ThesecondarystructureoftheAzoarcusribozymewith 37
3
nucleotidesnumberedaccordingtointronicnucleotides.Randomizednucleotides(N)arehighlightedinred.TheRNAsubstrateusedtoselectfor /6
5
reverse-splicingactivityisshownaslowercase.Thecurvedarrowindicatesthereactionusedforselection,whichresultsinaportionofthesubstrate 0
2
appendedtothe30-endofactiveribozymes.(B)Experimentalconcept.Theinvitrotranscribedribozymelibrary(bluetube,crystalstructure1ZZN) 2
8
9
isreactedindifferentmagnesiumconcentrations(redtubes).High-throughputsequencedataareusedtogenerateRNAfitnesslandscapesby b
determiningtheRNAFitness(enrichment/depletion)ofeachribozymevariantateachmagnesiumconcentration[Mgþþ]. y g
u
e
s
sequence largely determines how much of that sequence dataanalysistoolthatwasusedtoexplorethedata(available t o
n
reacts and then gets amplified by RT-PCR (supplementary at https://haydenlab.gitlab.io/azo). The results showed that 0
5
fig.1,SupplementaryMaterialonline).Sequenceswithalarger relative fitness of many sequences in the library increased A
p
increaseinfrequencyhavehigherfitness,andthosethatde- with increasing magnesium concentration (supplementary ril 2
crease in frequency have low fitness. The frequency of each figs. 6, 7, and 10, Supplementary Material online). Both the 0
2
2
sequencewas normalized to the frequency of thewild-type maximumfitnessobservedandtheaveragefitnessobserved
Azoarcus ribozyme at that magnesium concentration. The increased with increasing magnesium concentration (fig. 3).
activity of the wild-type ribozyme at each magnesium con- Thenumberofsequenceswithfitnesshigherthanwild-type
centrationwasalsodeterminedbyPAGE(supplementaryfig. (fitness>1)alsoincreasedwithhighermagnesium.However,
4,SupplementaryMaterialonline).Fitnessvalueswerehighly severalsequenceswithhigherthanwild-typefitnessathigher
correlated between replicates (supplementary fig. 3, concentrationsofmagnesiumhadmuchlowerfitnessatlow
Supplementary Material online) and with low coefficients magnesium(oneofthesegenotypes ishighlighted infig. 2).
ofvariationthatincreasedforthesequenceswiththelowest Specifically, several sequences with three or four mutations
activity(fitness<0.15)(supplementaryfig.11,Supplementary were more active than wild-type only at the two highest
Materialonline). magnesium concentrations (16 and 48mM magnesium).
Tovisuallyevaluatechangesintheribozymefitnessland- Theseresultssuggestthatcertaincombinationsofmutations
scapecausedbydifferentmagnesiumconcentrations,fitness canbebothbeneficialforribozymeactivityyetdestabilizing
was plotted as a function of mutational distance from the for the ribozyme structure, requiring more magnesium to
wild-type sequence (fig. 2). We also created an interactive appearasbeneficial.
3
MBE
Perietal. . doi:10.1093/molbev/msab373
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
/3
9
/3
/m
s
a
b
3
7
3
/6
5
0
2
2
8
9
b
y
g
u
e
s
t o
n
0
5
A
p
ril 2
0
2
2
FIG.2.TheRNAfitnesslandscapesunderdifferentmagnesiumionconcentrations.Foreachgraph,nodesrepresentindividualsequencesplotted
withmutationaldistance(xaxis)andfitness(yaxis)relativetothewild-typeribozyme.Eachnodeshowstheaveragefitnessofthreeexperimental
replicates.Nodesareconnectedbyablueedgeifthetwosequencesdifferbyasingle-nucleotidedifference.Violinplots(gray)representthe
distributionoffitnessvaluesateachmutationaldistance.Themagnesiumchlorideconcentration(MgCl )usedforeachdatasetisshownabove
2
eachgraph.Asequenceishighlighted(yellowedges)toillustratemagnesiuminducedchangestofitness.
Next,weaskedwhethertheruggednessofthelandscapes set (Szendro et al. 2013). Sign epistasis occurs when one
changed with magnesium concentration. As a first ap- individualmutationisdeleteriousbuthasabeneficialeffect
proach to measure ruggedness, we determined the fre- in the presence of a second already beneficial mutation.
quency of sign and reciprocal sign epistasis in each data Reciprocal sign epistasis occurs when both individual
4
MBE
DynamicRNAFitnessLandscapesofaGroupIRibozyme . doi:10.1093/molbev/msab373
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
/3
9
/3
FIG.3. Magnesiuminducedchangestothefitnesslandscapes.(A)Themeanfitnessofallvariantsinthelibrary,relativetowild-type,ateach /m
s
magnesiumionconcentration([MgCl ]).Threedatapointsateachconcentrationrepresentexperimentalreplicates.Thebluelinetracesthemean a
2 b
valueofthetriplicatesateachmagnesiumionconcentration.(B)Landscaperuggednesscalculatedasruggedness¼2(fRS)þfS,wherefRSisthe 373
fractionofreciprocalsignepistasisandfSisthefractionofsignepistasisforallpairsofmutations.(C)Meanfitnessplottedasafunctionof /65
mutationaldistancefromwild-type.Datapointsrepresentthemeanofthefitnessvaluesofallgenotypeswiththatmanymutationsatthat 0
2
magnesiumconcentration.Linesarebest-fitcurvestotheequationw(n)¼exp((cid:2)anb).(D)Effectofmagnesiumconcentrationontheexponential 28
9
decayparameteraanddirectionalepistasisparameterbfromthecurvefitting. b
y
g
mutationsaredeleterious,butcombinedtheyarebeneficial. deviation from a purely additive model by fitting a multi- ue
s
Thissevereformofepistasisisthecauseofpeaksandvalleys dimensionallinearmodeltothedata.Wefoundthatr/salso t o
n
infitnesslandscapes.Eachlandscapewasdividedintofour- decreasedwithincreasingmagnesium(supplementaryfig.9, 0
5
genotype subgraphs (squares) containing a starting geno- Supplementary Material online). Bothmeasurements indi- A
p
type,twosinglemutantsofthisgenotypeandthegenotype cate that increasing magnesium decreased ruggedness in ril 2
with both individual mutations combined. The fitness of thislocalribozymefitnesslandscape. 0
2
2
thefourgenotypesinthesquarewasusedtoidentifyinstan- Tofurthercharacterizethefitnesslandscapes,weanalyzed
cesofepistasis.Specifically,thefractionofsquareswithsign howtheinteractionsbetweenmostlydeleteriousmutational
epistasis and reciprocal sign epistasis was determined. We effects changed with magnesium concentration on average.
used the formula ruggedness¼2(f )þf , where f is the Forthisanalysis,wefittheequationW(n)¼exp((cid:2)anb)tothe
RS S RS
fractionofreciprocalsignepistasis andf isthe fractionof fitness (W) by mutational distance(n)dataateachmagne-
S
signepistasis(Wuetal.2016).Thedatashowedasignificant siumconcentration.Thisequationhasbeenpreviouslyused
negative correlation between ruggedness and magnesium todetecttheaveragedeleteriouseffectofmutations(a)and
concentration (Spearman correlation coefficient r¼(cid:2)0.94, any predominant direction of epistasis (b) (Bendixsen et al.
P¼2.15(cid:3)10(cid:2)11)(fig.3B).Withincreasingmagnesiumthere 2017).Wefoundthattheseparametersshowaninversere-
was a decrease in the fraction of squares that exhibit lationshipinresponsetomagnesiumconcentration(fig.3D).
signepistasisandreciprocalsignepistasis.Wealsoassessed The deleterious effect of individual mutations a decreased
changes to ruggedness using the roughness to slope withincreasingmagnesium,whichisexpectedasmagnesium
ratio, r/s (Carneiro and Hartl 2010; Szendro et al. 2013). shieldsthedestabilizingeffectofmostmutations.Valuesforb
This approach measures the level of epistasis as the were positive at all magnesium concentrations, indicating a
5
MBE
Perietal. . doi:10.1093/molbev/msab373
predominanceofnegativeepistasiswheremutationstendto genotypes emerge and replace lower fitness genotypes.
have a more deleterious effect in combination than what Thesesimulationshavetheadvantageoftakingintoaccount
would be expected if the individual effects were additive. both the frequency and the magnitude of epistatic interac-
Valuesforbincreasedwithincreasingmagnesium,indicating tionswhilealsoallowingstochasticprocessesthatmimicnat-
that whilehigher magnesiumreduces individualmutational uralevolution.
effects,thehighermagnesiumconcentrationsstudiedcould To characterize how magnesium-induced changes in
not completely shield the destabilizing effect of increasing ruggednesscanalterevolutionarydynamics,wemeasured
numbers of mutations. The inverse relationship between theratesofadaptationduringevolutionarysimulationson
theaverageeffectofindividualmutationsandtheircombined the different landscapes. The simulations started with
effects has been observed previously in antibiotic resistance 1,000individualsallwiththesamegenotype.Thisstarting
proteinswhereantibioticconcentrationwasvaried(Bershtein genotype was chosen because it had low fitness and was
etal.2006).Ourresults,observedundervaryingmagnesium seven mutations away from all the highest fitness geno-
concentrations, further support the generality of this typesfromeachmagnesiumconcentration.Theevolution
D
phenomenon. ofthispopulationwassimulatedfor2,000generationswith o
w
Magnesium titrations have been used extensively to un- a mutation rate of 0.01 and the population size was kept n
lo
derstand RNA folding because the magnesium required to constant at 1,000 individuals. The average fitness of the a
d
e
achieve the native folded state is related to the thermody- population was recorded at each generation. For each d
namicsoftheunfoldedtofoldedtransition. Inourdata,in- landscape, the 2,000-generation simulation was repeated fro
m
dividual active ribozymes showed a transition from lower 100times.Whenthe100simulationsfromaspecificmag- h
tohigherfitnesswithmagnesium,andwefittheHillequation nesium concentration are combined, we observed a sig- ttps
tothesedata(fig.4).Wewereabletofitacurveandextract moidalincreaseinfitnessovertime.Asasimplifiedmethod ://a
c
best fit parameters for all active ribozymes that achieved to extrapolate the growth rate of each fitted curve, we ad
e
a maximum fitness greater than fitness¼0.15. We found arbitrarily set a threshold at generation 50 and recorded m
ic
thatonaverage,themidpointofthetransitiontendstoin- themeanvalueofpopulationfitnessatthatgenerationfor .o
u
crease as a function of mutational distance from wild-type. each magnesium concentration. This simplified growth p
.c
Thisfurtherindicatesthatmostmutationshaveadestabiliz- ratetendedtoincreaseastheconcentrationofmagnesium om
ing effect, and on average sequences with more mutations is increased (fig. 5I), indicating that populations were /m
b
require more magnesium to stabilize the native, active fold. evolving faster at higher magnesium concentration. e/a
Theresultsindicatethattheeffectofmagnesiumtitrationson Furthermore, we fit the results from the 100 replicates at rtic
RNAfoldingcanbeevaluatedbydeepsequencingofinvitro each magnesium concentration to the logistic growth le/3
9
selections,whichcanenablethestudyofRNAstructureand equation: f(x)¼L/(1þe^((cid:2)k(x(cid:2)x0))) (fig. 5A–H), where /3
foldingofmanydifferentRNAsequencesinahighlyparallel L is the maximum value reached by the curve, k is the /m
s
a
fashion. “growth rate,” and x is the midpoint of the sigmoid. We b
0 3
The observed changes to the ruggedness of the fitness found that theLvalues obtainedfor eachlandscape indi- 73
landscapes suggest that changes in magnesium could alter cated that the populations were reaching the highest fit- /65
0
evolutionary processes. For example, the general smoothing ness genotype of the landscape by the end of the 22
8
of thelandscapes at higher magnesium indicates that more simulations (fig. 5J). However, the k values for each land- 9
b
mutational pathways to higher fitness become “accessible,” scape increased with magnesium indicating a more rapid y g
which could promote more rapid evolutionary adaptation. increase in fitness with higher magnesium (fig. 5K). The ue
s
However,evolvingpopulationsaresubjecttorandommuta- midpoints (x0) also indicated more rapid adaptation at t on
tionsandthereforeencounterbothoptimalandsuboptimal higher magnesium. The strong decrease in x0 with in- 05
intermediates in a stochastic fashion. To better understand creased magnesium indicates that less time (generations) A
p
the evolutionary consequences of the changing RNA land- was required to reach half the maximum fitness with ril 2
scapes,weusedcomputationalsimulationsofpopulationsof highermagnesium.Wefoundthatthepopulationsneeded 0
2
2
ribozymescompetingforsurvival.Eachsimulatedevolutionis more than twice as long to reach half maximal fitness at
runforasetnumberof“generations.”Foreachgeneration,a 1mM Mgþþ (x >70) as compared with 48mM Mgþþ
0
populationof1000individualsischosenfromthegenotypes (x <30)(fig.5L).Thesesimulationresultssupportthecon-
0
presentinthepreviousgeneration.Individualsarechosenat clusion that the general smoothing of the landscape
randomtoreproduce,buttheprobabilityofreproductionis causedbyhighermagnesiumcanincreasetherateatwhich
relative to the fitness of that genotype in the given magne- populationsevolvetoward higherfitness.
siumenvironment(seeMaterialsandMethods).High-fitness These simulations were also repeated using randomized
genotypesaremorelikelytoreproduce.Eachindividualthat startinggenotypesandpopulationsized,butnonotabledif-
reproduces also has a random chance of acquiring a single ferences were observed from the simulations with a fixed
mutation that converts it to a neighboring genotype. If this starting genotype and population size, except for growth
genotypeisofhigherfitness,thehigherfitnessgenotypehasa ratesbeingindependentfromthemagnesiumconcentration,
higherprobabilityofsurvivalinthesubsequentgenerations. likely caused by a considerable amount of simulations
Through this process of reproduction with mutation, the starting from high-fitness genotypes (supplementary fig. 12,
population tends to increase in fitness as higher fitness SupplementaryMaterialonline).
6
MBE
DynamicRNAFitnessLandscapesofaGroupIRibozyme . doi:10.1093/molbev/msab373
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
FIG.4.Magnesiuminducedchangestoindividualribozymevariants.(A)Exampleribozymeswithdifferentnumbersofmutations.Curvesarethe ic.o
Hillequationfittothefitnessdata.Thecolorofthelinesrepresentsthegeneticdistancefromwild-type,whichisindicatednexttoeachcurve.(B) u
p
Curvefitstoallsequencesinthedata,plottedbymutationaldistance.Thelinecolorsrepresentthemidpointofeachcurve(gradientfromblueto .c
o
m
redasthemidpointmovesawayfrom1mM).Onlythecurveswithmaxfitnessbetween0.15and3,slopelowerthan5,andmidpointlowerthan16
/m
wereincluded.Distributionsofmaximumfitness(C),hillslopes(D),andmidpointsfromtheHillequationcurvefitting,withsequencescategorized b
e
bymutationaldistancefromwild-type. /a
rtic
le
Discussion 2007; Wagner 2008). Although under a constant selection /39/3
We have presented analysis of how a single environmental pressure, it is likely that populations would be attracted to /m
s
a
variablecanchangethetopographyofanRNAfitnessland- local fitness peaksthatdisfavor sequence spaceexploration. b
3
scapeandalterevolutionaryprocesses.Becausetheenviron- Ourresultsillustratehowchangestoasingleenvironmental 73
mental variable studied impacts the stability of the active variable, such as Mg2þ ions, can alter the future adaptive /65
0
structure of these molecules, we anticipate that the general potentialofanevolvingpopulation. 2
2
8
observationsmaybesimilartoothervariablesthatalterthe Theresultspresentedheredemonstrateamethodforhigh- 9
b
foldingandstabilityofmacromolecules,suchastemperature, throughput biophysicalcharacterizationof RNA. Our results y
g
ionicstrength,crowding,thepresenceofdenaturantsormo- showthatmagnesiumtitrationsofribozymescanbestudied ue
s
lecular chaperones (Rutherford and Lindquist 1998; Desai by deep sequencing, enabling advancements in thestudy of t o
n
et al. 2014; Gracia and Russell 2014; Paudel and Rueda howRNAstabilityorfoldingchangeswithmutations.Assuch, 0
5
2014). In addition, because decreasing magnesium typically ourapproachjoinsothersinusingsequencingtechnologyto A
p
lowersRNAfitness,thedatapresentedherecanbecompared turn low-throughput biochemistry approaches into high- ril 2
withparametersinotherexperimentalsystemsthatincrease throughput approaches for characterizing RNA sequence to 0
2
selectionpressure,suchasincreasedantibioticsintheevolu- functionrelationships(PittandFerr(cid:2)e-D’Amar(cid:2)e2010).Forex- 2
tionofbeta-lactamasegenes(Bershteinetal.2006). ample,thedeepsequencingofribozymereactionsundervary-
Theobservationthatincreasingmagnesiumreducedrug- ing substrate concentrations has been used to measure
gednesswhilealsofavoringsequenceswithmoremutations enzyme kinetics for vast sequence libraries (Pressman et al.
hasimportantevolutionaryimplications.Theresultssuggest 2019;Andreassonetal.2020).Similarapproacheshavebeen
thatselectionunderhighermagnesiumconcentrationwould used to investigate RNA sequence requirements for RNA-
preservesequencesinthisdatasetathighermutationaldis- binding proteins (Guenther et al. 2013). One advantage of
tances,andthereducedruggednessmeansthatsomeofthe these approaches is that all sequence variants are studied
intermediatesonthemutationaltrajectoriestohigherfitness simultaneously in the same reaction, ensuring identical con-
would also be preserved. Together, periods of selection at ditions, which may improve the resolution of subtle differ-
higher magnesium would allow a population to explore ences between variants by reducing experimental noise. In
more sequence space in this local fitness landscape. Such addition, our approach used the enrichment of sequences
sequence space exploration is thought to be important for during in vitro selection, suggesting that the analysis can be
futureadaptationandevolutionaryinnovations(Amitaietal. adapted to other types of selections, whether in vitro or
7
MBE
Perietal. . doi:10.1093/molbev/msab373
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/m
b
e
/a
rtic
le
/3
9
/3
/m
s
a
b
3
7
3
/6
5
0
2
2
8
9
b
y
g
u
e
s
t o
n
0
5
A
FIG.5. Computationalsimulationsofpopulationsevolvingonthefitnesslandscapes.Evolutionarysimulationswerecarriedoutwithconstant p
populationsize(N¼1,000)andmutationalprobability(m¼0.01).About2,000generations-longsimulationswererepeated100timesoneach ril 2
0
landscape.(A–H)Populationmeanfitness(orangedots)ofthe100replicatesimulations(lightbluelines)oneachlandscapewasrecordedateach 2
2
generation, and mean values were fit to a logistic equation: f(x)¼L/(1þe^((cid:2)k(x(cid:2)x ))) (red lines). The best fit curves are shown for each
0
magnesium concentration. (I) The average fitness of the simulated populations at generation 50 (vertical dashed lines in A–H). This “G50
meanfitness”valueisplottedforeachmagnesiumconcentration.(J)Changeinthemaximumfitnessfromthelogisticfittedcurvescompared
withthechangeinmagnesium.(K)Increaseinthelogisticgrowthrate(k)ofthebest-fitsigmoidalcurveateachmagnesiumconcentration.(L)
Changeinthemidpoint(x )ofthebest-fitcurvesateachmagnesiumconcentration.Midpointsindicatethegenerationsrequiredtoreachhalfthe
0
maximumfitness.Lowermidpointsindicatefasteradaptation.
cellularbased(Townshendetal.2015;LiandZhang2018).We effects of mutations in a yeast tRNA were highly correlated
anticipate that future experiments will combine approaches between four different environments that altered tempera-
(i.e.,kineticsandmagnesiumtitrations)tobetterunderstand tureandRNAstability(3%DMSO).Thecomparisonoffitness
theinterplaybetweentheevolutionofstabilityandcatalysis. measurements in pairs of environments typically showed a
Itisinformativetocompareourresultswithrecentfindings Pearson’scorrelation>0.90,enablingthepredictionoffitness
oftheeffectofenvironmentalchangetothefitnesslandscapes inoneenvironmentbasedonfitnessmeasuredinadifferent
oftRNAinyeast.LiandZhang(2018)foundthatthefitness environment.IntermsofRNAstructure,alteringmagnesium
8
MBE
DynamicRNAFitnessLandscapesofaGroupIRibozyme . doi:10.1093/molbev/msab373
ions should have a similar effect to altering temperature, in Darwin,anditiswellacceptedthatchangesintheenviron-
thatunstablemoleculeswillnotfoldproperlyatlowmagne- mentarecriticalfortheemergenceandoptimizationofgene
sium or high temperature. However, in our data, we only functions. Experimental fitness landscapes are now offering
observesuchhighcorrelationswhencomparingsimilarmag- quantitative, molecular level insight into the processes of
nesiumconcentrations(supplementaryfig.5,Supplementary evolution.Thesequencespaceofgenomes,andmostgenes,
Materialonline).Someoftheribozymevariantsstudiedhere issovastthatitwillnotbeexhaustivelyexploredexperimen-
change from very low fitness to very high fitness over the tally,especiallywhenthisiscombinedwiththecombinatorial
magnesium concentrations studied, which contributes to complexity of environmental variables. The goal remains to
low correlations. This further suggests that the magnesium extrapolategeneraltrendsfromexperimentaldata,requiring
concentrations used in our experiments may cover a larger more significant collaborations between theory and experi-
rangeofthermodynamic“stress”thanthetemperaturediffer- ment.Coordinatedadvancementsonbothfrontsareneeded
encesusedintheyeastexperiments(23(cid:4)C,30(cid:4)C,or37(cid:4)C), tomoreefficientlyandaccuratelyexploreandpredictfunc-
whichisoneimportantdifference.Inaddition,themutational tionsingenotypespace.
D
librarieswerequitedifferent.ThetRNAlibraryusedrandom o
Materials and Methods w
mutagenesistoexhaustivelycoversingleanddoublemutants n
lo
fromthewild-typesequence. The Azoarcus ribozymelibrary LibrarySynthesis ad
e
hfreorme ftohceusweidldo-tnypsee,qwueitnhcesesvwenithrahnidgohmernmuuctleaotitoidnealpdoisstitainonces A DNA template was chemically synthesized with the d fro
sequence: m
andallcombinationsofmutationsatthesesevenpositions.If h
mostmutationshaveadestabilizingeffect,thenitwouldbe GTGCCTTGCGCCGGGAAACCACGCA ttps
expected that higher mutational distances would also have AGGGATGGTGTCAAATTCGGCGAAACC ://a
c
morevariedresponsetoenvironmentalchangesthatimpact TAAGCGCCCGCCCGGGCGTATGGCA ad
e
thermodynamicstability.Clearly,yeastcellscontainnumerous ACGCNGAGCCAAGCT m
ic
geneproductsthatchangewithtemperatureandthatcould TCGGCGCCTGCGCCGATGNAGGT .o
u
havetheeffectofbufferingtheeffectsoftemperatureonRNA p
GTAGNGACTAGACGGCACCCACC .c
structureandfunction.However,thedifferencesinthelibrary TAAGGCNAACGNTATGGTGAAG om
designandexperimentalconditionsareaninterestingpossi- /m
GCATAGTCNAGGGAGNGGCGAAAGTCACAC b
biliDtyetshpaittectahnehinigshp-irtehrfouutughrepuextpneartiumreenotfso.urexperiments,the AAACCGG, e/artic
data only explores a very small region of the total genotype whereNindicatesabasesynthesizedwithequalamountsof le/3
9
spacenearthewild-typeAzoarcusgroupIribozyme.Thislocal all four phosphoramidites (“machine mixed,” 4nmol /3
landscape data may be most informative for modeling and Ultramer, desalted, IDT). The template begins with the first /m
s
a
predicting the optimization of a known molecular function. three nucleotides of the ribozyme internal guide sequence b
3
Inrelatedworkbyothers,thesequencespaceofsmallerRNAs (IGS)andismissingthefirsteightnucleotidesofthenatural 73
was comprehensively evaluated for their ability to catalyze a intron sequence. This ssDNA template was PCR amplified /65
0
specificreaction.Suchcomprehensivelandscapesmaybemore with primers (table 1: Tas2.1a and T20a) that extend the 22
8
informative for understanding the origins of RNA functions templateby28basepairsthatincludethepromoterforT7 9
b
(Pressmanetal.2019).Thispriorworkfoundseverallocalpeaks RNA polymerase. PCR reactions contained 200ml 2(cid:3) Kapa y g
thatwereseparatedbyfitnessvalleysofseveralmutationsthat HifiMaster Mix(Roche),12mlofprimer1 (10mM),12mlof ue
s
causedundetectableactivity.Therewerenumerouspathways primer 2 (10mM), 16ml of library DNA template (100nM), t o
n
tohigherfitnessnearalocalpeak,butthepathwaysbetween and water to reach a reaction volume of 400ml (8(cid:3)50ml 0
5
peaks were not obviously accessible by evolutionary events, reactions). The library DNA was amplified for 18 cycles of A
p
suggestingthatevolutionwouldrarelyfindthemostoptimal 98 (cid:4)C for 20s and 74 (cid:4)C for 30s, and a final extension at ril 2
sequence(globalpeak).Ourresultssuggestthatenvironmental 72(cid:4)Cfor2min.Amplificationwasconfirmedona2%agarose 0
2
2
change could help bridge fitness valleys in RNA fitness land- gel. PCR products were purified on silica-based columns
scapes, similar to what has been found in the evolution of (DCC,Zymo),elutedinRNasefreewater(IDT),andquanti-
antibiotic degrading protein enzymes studied in Escherichia fiedbyA260(NanoDrop).TheRNAwastranscribedintrip-
coli during oscillating levels of antibiotics (Steinberg and licate (4h, 37 (cid:4)C) using T7 RNA polymerase (Thermo
Ostermeier 2016). Although both local and comprehensive Scientific). Transcriptions were purified with TRI reagent
fitnesslandscapesareinformative,moreworkwillberequired (Direct-zol, Zymo), separated on denaturing PAGE (10%,
tomergethesetwoapproachesinordertomorecomprehen- 8M urea), excised by UV shadowing, and eluted in 0.3 M
sively cover larger sequence spaces and form general or pre- sodium acetate. The elution was filtered (0.2mm), ethanol
dictivemodelsoftheeffectofenvironmentalchangesonthe precipitatedandresuspendedinRNAsefreewaterandquan-
originsandevolutionofmolecularfunctions. tifiedbyA260(Nanodrop).
Ourresultsjoinagrowingofbodyofdatathatexplorehow
fitness landscapes can be used to understand and predict InVitroSelectionsandNGSLibraryPreparation
evolution(deVisserandKrug2014).Theimportanceofthe A 30-nt RNA oligo was chemically synthesized with the se-
environment has been central to evolutionary theory since quence GGCAUAACUUCAAAUACUCCUCCCAACUCA
9
MBE
Perietal. . doi:10.1093/molbev/msab373
Table1. PCRPrimersUsedinThisStudy.
Name Sequence(50–30)
Tas2.1a CTGCAGAATTCTAATACGACTCACTATAGTGCCTTGCGCCGGGAA
T20a CCGGTTTGTGTGACTTTCGCC
Tas1.6a_illumina GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGTGAGTTGGGAGGAGTATTTGAAGTT
Illumina-reverse GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Tas2.1a_phased1 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGCATGCATGCATGCATGCGTGCCTTGCGCCGGGAA
Tas2.1a_phased2 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGTGCATGCATGCATGCGTGCCTTGCGCCGGGAA
Tas2.1a_phased3 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGATGCATGCATGCGTGCCTTGCGCCGGGAA
Tas2.1a_phased4 TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCATGCATGCGTGCCTTGCGCCGGGAA
(IDT). This substrate was designed to contain the compli- beforeselection. This valuewas normalized to wild-type¼1.
mentoftheIGS(CAU),havenopredictedsecondarystruc- The interactive data analysis tool was built with JavaScript
D
ture(NUPACK),andhaveasuitable primerbinding Tmfor usingacustomlibrarytoenablefastrenderingontheweb. o
w
reversetranscription.Ribozymereactionswerecarriedoutin n
lo
triplicatewith 10pmoloflibrary RNA, buffer(EPPSpH 7.5) CalculatingLandscapeRuggedness a
d
e
andvariousconcentrationsofaddedMgCl (1,2,3,4,6,8,16, Pairwiseepistasiswasdeterminedforeverysetoffourgeno- d
and 48mM). RNA was heat cooled prior2to adding magne- typesthatcomprisemutational“squares” inthelandscapes. fro
m
sium. The reactions were initiated with the addition of 50 One genotype was defined as a reference, with a starting h
wpmitholeEsDsTuAbsteraqtuei,minoclaurbattoedthaet 3m7a(cid:4)gCnefsoiurm1hc,oanncdensttoraptpioend. nfitonteyspseWwa0.sTcohme fpitanreedsswdiitfhfetrhenecfietnbeestswoefeangethniostryepfeerwenitchetwgeo- ttps://a
c
Reacted RNA was heated to 65 (cid:4)C for 5min and cooled to mutationaldifferences(WAB)andeachoftheindividualmu- ade
room temperature with a primer (table 1: Tas1.6a-illumina) tationaldifferences(W andW ).Asquarewasconsideredto m
A B ic
that was complementary to the substrate and added se- have no epistasis if WABþW0(cid:2)WA(cid:2)WB¼0. For all the epi- .ou
quence compatible with NGS Unique dual index primers staticsquares,twoconditionswereverified: p.c
o
(IDT) to the first strand RT product. Reverse transcription m
(ProtoscriptII,NEB)reactionswereincubatedat42(cid:4)Cfor1h, Condition 1: jWA– W0þWAB– WBj /m
then 65 (cid:4)C for 25min. Second-strand synthesis was carried ¼ðjWA– W0jþjWAB– WBjÞ be/a
outby lowcyclePCR(10 cycles,KAPAHiFi, Roche) using a rtic
le
mixture of four different forward primers (table 1: Tas2.1a- /3
9
pInhdaesxedco1m-4p)attoibaldedsesqeuqeunecnecetoditvheerssiteycoannddsNtrGaSndUn(BiqeunediDxsueanl Condition 2:<jðWjWB–B–WW0þ0jþWjAWB–ABW–AWj AjÞ: /3/ms
a
et al. 2020). This second-strand synthesis PCR product was b
3
purified with silica-based columns (DCC, Zymo) and ampli- Thesquareswherecondition1wasTrueandcondition2 73
fied with Illumina compatible NGS Unique Dual Index wasFalsewereclassifiedasmagnitudeepistasis.Thesquares /65
0
Primers (IDT). A unique index combination was used for wherecondition1isFalseandcondition2isTruewereclas- 22
each replicate. PCR products were quantified by qPCR and sified as reciprocal sign epistasis. All the remaining squares 89
b
pooled equimolar. The pooled DNA was assessed with a wereclassifiedassignepistasis.Theruggednessofalandscape y g
Fragment Analyzerandsequenced on aHiSeq4000 using a at a given magnesium concentration was calculated as ue
s
pairedend150reads(GC3F,UniversityofOregon). ruggedness¼2(fRS)þfS, where fRS is the fraction of squares t o
withreciprocalsignepistasisandfSisthefractionofsquares n 0
5
CalculatingFitnessfromSequenceDataAnalysis withsignepistasis.Roughnesstosloperatio(r/s)wascalcu- A
p
PTahiereadbeunnddarenacdesowfeeraechjosineqeduewnictehvFaLrAiaSnhtainllotwheinrgeafodrd“aotuatiaensd.” lqauteadntaifsiedseshcoriwbewdeplrleevaiochusglyen(Aoittyapeetafilt.n2e0s0s1c).aBnribefelyp,trhedeircatteido ril 20
2
curve fitting was carried out using custom python scripts fromtheeffectsofmutationsonthebackgroundofthemost 2
available on GitLab (https://gitlab.com/haydenlab/azo7). fitgenotype,assuminganadditivemodeloffitnesseffects.
Briefly,readswerecountedthatmatchedaregularexpression
Supplementary Material
of the wild-type ribozyme sequence with wild-cards at the
variable positions. For each matching read, the nucleotide Supplementary data are available at Molecular Biology and
identityateachvariablepositionwasextractedtodistinguish Evolution online. Raw sequence data in fastq format can be
betweenthesequencevariants.Allsequenceswithmutations foundattheEuropeanNucleotideArchive(accessionnum-
outside the expected positions were excluded from further berPRJEB43384).
analysis. For each replicate, the counts of each unique se-
Acknowledgments
quence variant were divided by the total count of all reads
that matched to any variant. This resulted in proportion of WeacknowledgesupportfromBoiseStateUniversity(E.J.H.,
readsrepresentedbyeachsequence.Fitnesswasdetermined BiomolecularSciencesGraduatePrograms),NationalScience
bydividingtheproportionobservedafterselectioninaspe- FoundationOfficeofIntegrativeActivities(E.J.H.,GrantNos.
cific magnesium concentration by the proportion observed OIA-1738865 and OIA-1826801), National Aeronautics and
10