Evolutionary Potential of a Duplicated Repressor-Operator Pair Simulating Pathways Using Mutation Data

Evolutionary Potential of a Duplicated Repressor-Operator Pair Simulating Pathways Using Mutation Data EvolutionaryPotentialofaDuplicated Repressor-OperatorPair:SimulatingPathways UsingMutationData Frank J.Poelwijk[,Daniel J.Kiviet[,Sander J.Tans* FOMInstituteforAtomicandMolecularPhysics(AMOLF),Amsterdam,Netherlands Ample evidence has accumulated for the evolutionary importance of duplication events. However, little is known about the ensuing step-by-step divergence process and the selective conditions that allow it to progress. Here we presentacomputationalstudyonthedivergenceoftworepressorsafterduplication.Acentralfeatureofourapproach is that intermediate phenotypes can be quantified through the use of in vivo measured repression strengths of Escherichiacolilacmutants.Evolutionarypathwaysareconstructedbymultipleroundsofsinglebasepairsubstitutions andselectionfortightandindependentbinding.Ouranalysisindicatesthatwhenaduplicatedrepressorco-diverges together with its binding site, the fitness landscape allows funneling to a new regulatory interaction with early increasesinfitness.Wefindthatneutralmutationsdonotplayanessentialrole,whichisimportantforsubstantial divergenceprobabilities.Byvaryingtheselectivepressurewecanpinpointthenecessaryingredientsfortheobserved divergence.Ourfindingsunderscoretheimportanceofcoevolutionarymechanismsinregulatorynetworks,andshould berelevantfortheevolutionofprotein-DNAaswellasprotein-proteininteractions. Citation:PoelwijkFJ,KivietDJ,TansSJ(2006)Evolutionarypotentialofaduplicatedrepressor-operatorpair:Simulatingpathwaysusingmutationdata.PLoSComputBiol 2(5):e58.DOI:10.1371/journal.pcbi.0020058 conditions that favor the evolution toward independent Introduction regulation. Interestingly, such regulatory divergence is Initially put forward by Stevens in 1951 [1] and later inherently a coevolutionary process, where repressors and advocatedbyOhnoinhisseminalwork[2],geneduplication operatorsmustbeoptimizedinacoordinatedfashion. followed by functional divergence is now seen as a general The mere presence of a selective pressure is clearly not a mechanism for acquiring newfunctions [3]. Also,regulatory suf,cient condition to achieve a new function. Rather, the networks are thought to be shaped signi,cantly by genetic evolutionary potential and limitations can be seen as duplication [4]. For instance, sequence analysis of tran- governed by the shape of the actual ,tness landscape and scriptionfactorfamiliespointstovarioushistoricalduplica- the evolutionary search within it. Studying these intrinsic tion events [5,6]. However, very little is known about the limitations to divergence represents the overall aim of this subsequent mutational divergence pathways or about the work. Many open questions arise when considering the correspondingstepwisephenotypicalchangesthataresubject formation of a new protein-DNA interaction, which may be to selection. While these issues have not yet been explored viewed as the construction of a new lock and uniquely experimentally, related generic aspects of mutational plasti- matching key. For instance, how should the protein be city have been addressed theoretically [7–11]. However, a modi,ed step-by-step to recognize a new DNA-binding site central obstacle in studying mutational pathways through thatalsodoesnotyetexist,orviceversa?Onewouldexpect computersimulationsremainstheunknownrelationbetween thatcomplementarymutationsneedtooccurintheprotein the sequence and binding af,nity, for which, in general, a andDNA-bindingsite.Doesthismeanthattemporarylosses ratherabstractmapping hastobeassumed.Todescribethe in,tnessmustbeenduredwhentakingsingle-mutationsteps? formationofanewregulatoryinteractionafteraduplication And,howmanymutationsmustminimallyaccumulatebefore event, which is our current aim, such an abstract approach anoticeablenewrecognitionisobtainedonwhichselection wouldbeparticularlyspeculative. Herewereasonthatmanycharacteristicsoftheadaptation Editor:SarahTeichmann,MRCLaboratoryofMolecularBiology,UnitedKingdom of real protein-DNA contacts are hidden in the extensive ReceivedFebruary17,2006;AcceptedApril12,2006;PublishedMay26,2006 body of mutational data that has been accumulated over many years (e.g., [12–14] for the Escherichia coli lac system). ApreviousversionofthisarticleappearedasanEarlyOnlineReleaseonApril12, 2006(DOI:10.1371/journal.pcbi.0020058.eor). These measured repression values can be used as ,tness landscapes,inwhichpathwayscanbeexploredbycomputing DOI:10.1371/journal.pcbi.0020058 consecutive rounds of single base pair substitutions and Copyright:Ó2006Poelwijketal.Thisisanopen-accessarticledistributedunder selection. Here we develop this approach to study the thetermsoftheCreativeCommonsAttributionLicense,whichpermitsunrestricted use, distribution, and reproduction in any medium, provided the original author divergence of duplicate repressors and their binding sites. andsourcearecredited. More speci,cally, we focus on the creation of a new and *Towhomcorrespondenceshouldbeaddressed.E-mail:tans@amolf.nl uniqueprotein-DNArecognition,startingfromtwoidentical [Theseauthorscontributedequallytothiswork. repressorsandtwoidenticaloperators.Weconsiderselective PLoSComputationalBiology|www.ploscompbiol.org 0467 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData motif[16].Mostoften,selectivepressuresforanewfunction Synopsis willbeabsent,inwhichcasesilencingofoneoftheduplicates isthemostprobableoutcome[3,17].However,therarecases Theevolutionofanewtraitcriticallydependsontheexistenceofa pathofviableintermediates.Generallyspeaking,fitnessdecreasing whereaselectivepressureispresentarecrucialtodeveloping stepsinthispathhamperevolution,whereasfitnessincreasingsteps newfunctions. accelerateit.Unfortunately,intermediatesarehardtocatchinaction Weaimedtode,neatransparentselectionpressureforthe since they occur only transiently, which is why they have largely divergence of these regulatory interactions. Attributing a beenneglectedinevolutionarystudies. ,tnessvaluetoanetworkfunctionisnon-trivial:unlikeforan enzymaticfunction,network,tnesscannotbecapturedina The novelty of this study is that intermediate phenotypes can be single biochemical parameter. Here we propose to assign a predicted using published measurements of Escherichia coli ,tnessvaluebasedonthedesiredinput-outputrelationofthe mutants. Using this approach, the evolution of a small genetic networkissimulatedbycomputer.Followingtheduplicationofone network (see Figure 1A and 1C). For simplicity, only two ofitscomponents,anewprotein-DNAinteractiondevelopsviathe concentration levels (high and low) of input and output accumulationofpointmutationsandselection.Theresultingpaths protein are considered, resulting in four possible input reveal a high potential to obtain a new regulatory interaction, in conditions. For each of these input conditions, it follows which neutral drift plays an almost negligible role. This study straightforwardly which repressor-operator interactions providesamechanisticrationaleforwhysuchrapiddivergencecan should be maximized and which must be minimized. The occur andunder which minimal selective conditions. Inaddition it interaction strength between operator Oi and repressor yields a quantitative prediction for the minimum number of homo-dimerRjisexpressedbyrepressionvalues(FOiRj ).This essentialmutations. valuerepresentstheexpressionlevelofadownstreamgenein theunrepressedconditiondividedbytherepressedcondition anditisobtaineddirectlyfrommeasureddata(seebelowand can act? The latter is an important point: mutations Materials and Methods). Taking the ,tness to scale linearly conferringaselectiveadvantagespreadmorereadilythrough with the repression values, the ,tness of the complete a population [15], resulting in a drastic increase of the networkisdenotedbytheproductofalloptimizationfactors: divergence probability. These questions are addressed by exhaustivelysearchingthelandscapeforoptimalpathways,as FO1R1 FO2R2 Fitness?maxðFO1 O2 maxðFð1Þ Þ FO1R2 Þ wellasbycomplementarypopulationdynamicssimulations. FO2R1 Previously it has been shown that lac repressor mutants Inthisexpressionmax(FOi)denotestherepressionvalueof indeedexistthatcanbindexclusivelytomutantlacoperators thestrongestinteractionwithOi,eitherbyhomodimersofR1 [14]. Our simulations reveal that a duplicated repressor- orR2orthehetero-dimercomposedofR1andR2(seeFigure operator pair can readily evolve to achieve such independ- 1andMaterialsandMethods). enceofbinding,whilemonotonouslyincreasingits,tnessin The,tnessde,nitioncomesdowntoaminimumsetoftwo astep-by-stepprocess.Moreover,simplyfollowingthe,ttest demands for regulatory binding: each operator must bind a mutantsdoespredominantlyguidethesystemtothedesired repressortightly(max(FO1)andmax(FO2)shouldbelarge)but global optimum, which indicates funnel-like features in the alsoexclusively(FO1R1/FO1R2andFO2R2/FO2R1shouldbelarge). ,tness landscape. A detailed analysis of the subsequent Priortodivergence the,rstdemand isalready met, but the network changes indicates a generic sequence of events, of latterisnot.Thechallengeduringdivergenceisthereforeto which we study the underlying mechanisms by varying the improvebindingexclusivity,whilemaintainingtightbinding. appliedselectivepressure.Next,weshowthatthetrajectories Tight and exclusive binding is a core functionality of most we ,nd in the optimal pathway simulations are not rare regulatory systems, and most pairs of existing transcription exceptions, since similar trajectories are followed using a factors must therefore score well on the employed ,tness probabilistic scheme for accepting a mutation. The results de,nition. Take for instance the LacI and RafR repressors, further suggest the feasibility of studying regulatory diver- which regulate enzymes required for growth on lactose and gence in laboratory evolution experiments, and ,nally we raf,nose,respectively.Ifoperatorbindingwouldnotbetight makeacomparisontoalternativemodelsforthecreationof in the absence of lactose and raf,nose, the wasteful newregulatoryinteractions. expressionofthedownstreammetabolicenzymeswouldlead tosub-optimalgrowthspeeds[18,19].IfRafRwouldalsobind Results tothelacoperator(andthusbindnon-exclusively),theeffect on growth speed would also be negative since the mere SelectivePressureandFitnessLandscape absence of raf,nose would then lead to insuf,cient b- Weconsideranecologicalsituationwherenaturalselection galactosidaseforhighlactoseconcentrations. would favor independent regulation of two genes X and Y. Regulation is not independent in the initial symmetric One therefore expects a conservative selective pressure network with duplicated components (see Figure 1): X and that minimally includes binding tightness and exclusiveness, to keep the lac and raf regulation intact. Important here is Y have two identical upstream binding sites (O1 and O2), which bind two identical repressors (R1 and R2) equally thatthelacandrafrepressorsareinfactrelated:theirorigin strongly. Such a situation will, for instance, arise upon has been traced to duplication events from a common ancestor [6]. If a conservative pressure keeps their function duplicationofarepressorthatregulatestwoormoregenes. Note that this selective pressure, of course, is not a generalintact now, it seems a good candidate for the initial outcomeofarepressorduplication.Aduplicationeventmay divergence pressure as well. Full divergence to the current ariseinthecontextofadifferentfunctionalpressure,which lacandrafsystemsclearlyinvolvesmanyadditionaldevelop- could direct the evolution toward a different topological ments after duplication. For instance, the divergence of PLoSComputationalBiology|www.ploscompbiol.org 0468 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData Figure1.DivergenceProcess,FitnessCriterion,andMutationalDatasetofRepressionValues (A)Diagramofthestudieddivergenceprocess:afteraduplicationevent,anewregulatoryinteractioncanbeformedbymutatingthetwooperators,O1 andO2,andtworepressors,R1andR2. (B)Duplicationanddivergenceyieldsheterodimers,whichcanallbindtotheoperator.The(initiallysymmetric)operatorsandrepressorsarebasedon thelacsequence,asindicated.Basepairsthatarekeytoalteringspecificity(coloredredandblue)canbemutatedtoarbitrarysequence. (C)Theselectivepressureforindependentregulationfollowsfromfourinputconditionsthatcontributetothetotalfitness.When,e.g.,R1ishighand O1R1 and(FO2R1 R2low,thisimpliesthatXshouldbelowandYhigh.Outofallinteractionparametersofthenetwork,inthiscaseonlyFÀ1 andneedtobeoptimized.WhenR1andR2arehigh,bothXandYshouldbelow,regardlessofwhichrepressor-dimercausesrepression.Therefore ) arerelevant max(FO1) (the strongest interaction with O1 by either homodimers of R1 or R2 or by the heterodimer of R1 and R2) and max(FO2) need to be be optimized.WhenbothR1andR2arelow,noparametersneedtobeoptimized. (D)Resultingrepressionvaluelandscape,showingrepressionvaluesbasedonactualmeasurementsofmutants. DOI:10.1371/journal.pcbi.0020058.g001 ligand-binding properties[20]might haveoccurred priorto restrictednumberofbasepairsthatcanbemutatedinsilico operator-bindingdivergence.Whiletheseconsiderationsput andlinkedtoexperimentaldata.Atthesametimehowever, additionalconstraintsontheentiredivergenceprocess,they while the tightness of DNA binding is the result of the do not alter the particular operator-binding divergence integralproteinarchitecture,surprisinglyfewbasepairs(ten studiedhere. intotal)havebeenfoundtobeimportantforalteringbinding speci,city [14] (see Figure 1B). Focusing on these key base A remaining question still is how the various demands pairs is therefore reasonable for the minimal paths that we shouldbe weighedin thetotal ,tness. Thatchoice isclearly areinterestedinhere.Usingmeasurementson1,286mutants, not general: it will strongly depend on the operons in repressionvaluesofallvariantsinthesekeybasepairscould questionandonthechangingcellenvironment.Forexample, convincingly be determined [13,14,21]. These variants thus ifactiveRafRispresentmorethanhalfofthetime,thenits include all multiple mutants in both the repressor and the cross-interaction with the lac operator would be compara- operator.Repressionvaluesofheterodimersandasymmetric tivelymoreharmfulbecauseitlastslonger.Inordertogivea operatorsarecalculatedusinganadditivecontributionofthe uniform presentation we weighed the factors of the four repressor monomers to the dimer-DNA binding [22] (see input states equally, which would correspond to an equal Materials and Methods). In total, about 1 3 107 possible contributionofthese phasestotheoverall,tness.However, repressor-operatorcombinationsareobtained(seeFigure1D weighingthefactorsunequally(e.g.,byincreasingthepower forthehomodimervariants). ofthetightoperatorbinding,orthecross-interactionfactors Everymutationalpathstartswiththeduplicatedsequence from1to2)didnotalterthemainconclusions. of a tight binding repressor-operator combination (repres- MutationDataandPathwaySimulations sionvalue 100).Thesepossiblestartingsequencesobviously In our simulations, the strength of a mutant repressor- include wild-type lac, but also e.g., the gal and ebg systems, operatorinteraction(asexpressedbytherepressionvalueF), which are part of the same family of repressors. Their high isassignedusingdatafrommutationalanalysis[14].Inthese measured repression values are rather remarkable because experiments,repressionvalueshavebeendeterminedinvivo therestofthegal,ebg,andlacsequenceshaveinfactdiverged from the repressed and unrepressed expression levels of a considerably.Theseobservationsfurtherindicatethatthekey lacZgene,controlledbyamutantlacoperatorandmutantlac basepairsplaythecentralroleinspeci,crecognition. repressor (see Materials and Methods). Obviously not all The aim of the simulation method (see Materials and possible lac mutants have been constructed. Therefore, a Methods for details) has been to reveal the intrinsic potentially signi,cant limitation of our simulations is the possibilitiesforthedivergenceofrepressor-operatorbinding, PLoSComputationalBiology|www.ploscompbiol.org 0469 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData giventhemeasureddataandtheconstraintsofsinglebasepair substitutions and no ,tness decreases. For this purpose, we searchthelandscapeforoptimalpathsandstudywhattheir limitationsandpotentialare.Totracetheseoptimalpaths,all mutants with a single base pair substitution with respect to their parents are evaluated based on the ,tness described above, and the best performers are selected for the next round.ThenumberofselectedmutantsLisvariedtoassessits effect.Wealsoquestionwhethertheseoptimalpathsarenot justrarecases,bycomparingthemwithpathwaysgenerated byadifferentsimulationmethod,wherearandommutationis accepted with a probability that depends on its associated ,tnessincrease[23](seeProtocolS1). ShortCo-DivergencePathways The simulations show that paths to independent recog- nitionarereadilyfound.Evenwhenonlythebestnetworkis carried to the next round (L ? 1), which implies always following the steepest ascent in ,tness, some starting sequences can evolve to the highest ,tness in the sequence space.Inthesenetworks,bothrepressorsbindtightlytoone operator (FO1R1 ?520 and FO2R2 ?200, respectively), while notatalltotheother(F O1R2 ?1,F ?1).Weconsidered paths to be successful when the ,tness value is within an O2R1 order of magnitude of the highest ,tness in the landscape, which is a strict criterion given the fact that the ,tness parameter is a product of six factors. The diverged fraction increases for higher L (Figure 2A, solid line), which is expected since it allows alternative paths to be explored. Figure2.DivergenceSuccessRatioandPathLengthDistributions Moresurprisingisthatsuccessfultrajectoriescaneventually (A) Fraction of starting sequences (numbering 132 in total) that befoundfromallstartingpoints,butnotethatpathsthatcan successfully diverge, as a function of the number of networks carried only be followed for higher L are increasingly less probable to the next round (L). Dashed line, idem, but with the additional requirementofcontinuedtightbinding(F100)forbothrepressors. becausetheyimplymore(near)neutralmutations. (B)Distributionofpathlengthsuntildivergence.Redcolormap,optimal Mostoptimalpathsarerathershort:70%requirejust,ve co-divergencepathways.Bluecolormap,pathwayswiththeadditional requirement of F 100 for both repressors. Note that a vertical toninemutationsforL?20(Figure2B).Thesystemsalmost summationofthecolormapsyieldsthelinesin(A). exclusively,ndthenearestdivergedstateinsequencespace DOI:10.1371/journal.pcbi.0020058.g002 (Figure3B)anddosowithouttakinganydetours(Figure3A). Notably, despite the fact that the starting points lie in very ces.Eventuallyallsuccessfultrajectoriesrecoverpalindromic differentareasofthesequence space,agenericsequence of operators, even as the selective pressure does not explicitly network changes is generally observed (see Figure 4 for an specifythis.Withalldimervarietiespresent,ahomodimeris example). First of all, one repressor-operator combination available and now binds most tightly to the palindromic remains unchanged, except at the very end, as the other operator. diverges away. This is an example of asymmetric divergence In order to obtain a better insight in the essential due to positive selection, as has also been found in ingredients for the observed evolvability, various additional phylogenetic analysis of duplicate genes in eukaryotes [24]. simulationswereperformed.Forinstance,weweretriggered A striking general feature of the pathways is an early by the recurrent early knockout of one of the repressors, reductioninthebindingstrengthofthedivergingrepressor, whichisoneofthemostnoticeablefeaturesofthemutational brought about by a single base pair substitution (Figure 4B, pathways. To test for the importance of this step, both redtrace).Suchamutationwouldbeunfavorableforasingle repressor-operator pairs were required to maintain a repressor-operator pair, but here it can be ,tness neutral, signi,cant repression (FO1R1 100 and FO2R2 100). partly because the unchanged duplicate repressor ensures a Divergence is indeed signi,cantly frustrated by these con- continued repression. At this speci,c point the diverging ditions (Figure 2A, hatched line). The amount of selected repressorisfreedfromfunctional constraintsandtherefore mutants needs to be two orders of magnitude larger (L . most vulnerable to degenerative mutations resulting in 1,000) for half of the starting sequences to diverge. The silencingofthegene.Theprobabilityofsilencingisreduced saturation of the diverged fraction for very high L, where however,becausealreadyatthesecondmutationandonward, prolongedneutraldriftisallowed,indicatesthatfor22%of newanduniqueprotein-DNArecognitioncanbebuiltup.At the starting sequences no pathways exist. Moreover, in the sequence level, this phase is characterized by transient contrast to the optimal paths, the nearest diverged state in asymmetries.Theoperatormustgothroughnon-palindromic the landscape is generally not found, and the paths contain sequencesbecauseitcanonlyreceiveonemutationatatime. signi,cant detours (Figure 3). The same is seen from the Heterodimers are the best binders in this phase because of increasedpathlength:70%ofthepathstake11–21mutations theirabilitytomirrorthenon-palindromicoperatorsequen- PLoSComputationalBiology|www.ploscompbiol.org 0470 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData Figure3.AnalysisofPathwayDetoursandLocalEnvironmentofFitnessOptima (A)Histogramshowingthenumberofdetourmutationsofthedivergencepathways.TheHammingdistancedH oftwosequencesisdefinedasthe numberofpositionsatwhichtheyhavedifferentbasepairs.PathsthatarelongerthandH arriveatanoptimumafteradetour. (B)HistogramoftheHammingdistancebetweentheoptimumthatisfoundandtheclosestoptimum.Ifthismeasureiszero,apathleadstotheclosest optimum. (C)MedianfitnessvalueasafunctionoftheHammingdistancefromaglobaloptimum(solidline).Graylevelsindicatethespreadofthefitnessvalues. DOI:10.1371/journal.pcbi.0020058.g003 (Figure2B).Thesepathslackarecurringmutationpatternas early on is important, since it greatly enhances divergence observed for the optimal paths and instead show a large probabilities [15]. A lack of early selection would result in variation in the sequence of events. Both repressors and much higher probabilities of silencing of one of the operatorsaresigni,cantlymutated,andthe,tnessincreases duplicatesbytheaccumulationofmutations[3,17]. slowlyorisneutralovermultiplerounds(seeFigureS1foran While the presented systematic search for optimal path- example). ways is useful in revealing necessary conditions for diver- gence,onemaywonderwhetherpathsarenotverydifferent Another de,ning feature of duplicated transcription in a probability-based ,xation process that typi,es natural factors is the heterodimerization of transcription factor evolution. However, we found that population genetics monomers.Itisnotapriorievidentwhetherthisconstraint simulations reveal the same pathway characteristics: a on the network topology either promotes or hampers divergence. To assess its effect, simulations were performed signi,cantfractionofpathsaresuccessfulwithmonotonous whereheterodimersarenotabletoform(unpublisheddata). ,tnessincreases,onerepressiondipearlyon,andfewneutral The results indicated a surprisingly limited effect on the mutationsarepresent(seeProtocolS1andFigureS2). divergence. The paths do initially show a slower ,tness The coevolutionary search for a new and independent increase,butthepathlengthdoesnotappearmuchaffected, recognition, which is relevant for both protein-DNA and nor the success rate of divergence. The other simulation protein–protein interactions, comprises fundamental differ- enceswithoften-consideredevolutionofligand-bindingand variationsweconducted(withunequally weightedfactorsin enzymatic activity [26–28]. While in the latter cases the new the ,tness de,nition), did not qualitatively alter the main evolutionary target is ,xed, here it is open-ended: as with divergence features, such as substantial divergence success locks and keys, many possible combinations are unique without,tnessdecreases,shortpaths,andanearlyrepression matches, and each of those is equally suitable. This large dip,indicatingtherobustnessofourresults. degree of freedom allows the system to choose the solution that is most accessible. Another difference with ,xed-target Discussion evolution lies in the selective pressure. Binding is already DuplicationandCoevolutionaryDivergence tight to both operators at the start of the coevolutionary Weobtaina,rstviewona,tnesslandscapeforregulatory scenario,sotheinitialpressuretochange,infactcomesfrom divergence that is based on actual measured data. We show bene,ts of not binding another operator. This pressure for that the landscape allows evolutionary paths toward inde- unique recognition is characteristic for regulatory interac- pendent repressor-operator interactions, exhibiting a step- tionsbutplaysmuchlessaroleindevelopingotherfunctions by-step increasing ,tness, starting as early as the ,rst or suchasenzymaticactivity.Thesecharacteristicsofacoevolu- secondmutation.Sincethepossibilityoffollowingsuchpaths tionarymechanism,togetherwiththeremarkableplasticityof critically depends on molecular properties, the use of protein-DNAinteractionsresultinahighlyevolvablesystem. empirical data is essential for such claims. One could also have imagined ,tness landscapes where paths to diverged FitnessLandscapeFunnels networksdonotexist,orwheretheyareverylong,involving The diversity of molecular architectures is not only large detours. Our results contrast with the notion that a constrained by their inherent physico-chemical limitations, number of neutral or even deleterious mutations have to but also by the existence of viable evolutionary routes that accumulate before a new function can develop (see for a shapethem.Forinstance,inapopulationofbacteriathereis discussion e.g., [25]). Having bene,cial mutations available only a small probability that an advantageous function PLoSComputationalBiology|www.ploscompbiol.org 0471 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData Figure4.TypicalDivergencePathway:NetworkChanges,Fitness,andSequence (A)Evolvinginteractionnetwork,wherelinethicknessdenotesbindingstrengthbetweenrepressormonomerandoperator-half.Dottedlinesdenote negligiblerepression.Yellowcrossesindicaterepressorandoperatormutations,whicharepositionedatthetopandbottomoftheinteractionlines respectively. (B)Fitnesstrajectory(black)andcorrespondingrepressionofeachrepressoronitsoperator(redandblue).Fitnessisnormalizedtothemaximumvalue 10 (;1310 ). (C)Sequencesforeachround.Mutatedpositionsarecoloredwhite. DOI:10.1371/journal.pcbi.0020058.g004 mutationdecreasesaninteractionthatshouldbemaximized, emergesifatemporary,tnessdecreaseisrequired,rst.Put differently,theshapeofthe,tnesslandscapeiscritical,and this negative effect on the ,tness is partly balanced by the onecanreadilyimagine,tnesslandscapeswheretheoptima decrease of an unwanted cross-interaction. A second mech- are very dif,cult to reach. Upon ,rst inspection, the anism is a coevolutionary twist on Ohno’s original idea, in measured landscape we consider indeed contains many which one repressor-operator pair can search for a new potential frustration sources: over 99% of all optima in therecognition, while the other repressor maintains repression landscape are in fact below our divergence criterion. Such onbothoperatorsintheveryearlystages.Aswehaveobserved local optima represent traps in which the system gets thatadropinthebindingstrengthisnecessaryforef,cient permanently stuck once it encounters one. However, the divergence, the ability to compensate for its negative results show that the system is still guided in the right contributiontothe,tnessiscrucialforfunneling. direction to (near) global optima, which indicates that the The evolutionary fate of redundant genes has previously ,tnesslandscapecontainsfunnel-likefeatures.Moreover,the been studied primarily using sequence analysis [3,29]. By optimalpathscontainnegligibledetours(Figure3A)andlead using a different dataset and approach, our simulations to the nearest optimum (Figure 3B), showing that the strengthen recent evidence for a more rapid ,xation of funneling is ef,cient and not constrained by ruggedness. A mutations in redundant genes [29] (termed ‘‘accelerated funnel-likeorganizationofthelandscapeisalsosupportedby evolution’’). Our analysis enables a next step in our under- the monotonous ,tness increases of the probabilistic path- standingofthisimportantprocess:Itprovidesamechanistic ways (Figure S2C), as well as by the smooth ,tness decrease rationaleforwhysucharapiddivergencecanindeedoccur, whensteppingawayfromaglobaloptimum(Figure3C). in terms of minimal selective conditions bacteria must Theunderlyingcausesforfunnelsinthe,tnesslandscape experience, in combination with independently measured maybefoundattwolevels. The,rstlevel isthatofasingle plasticityofprotein-DNAinteractions.Furthermoreityields repressor-operatorinteraction.Thesurfacesmoothnessthat a quantitative prediction for the minimum number of isneededforthefunnelsmaybepartlyunderstoodfromthe essentialmutationstoachievedivergence. reportedadditivecontributionsofthelacaminoacidstothe SuggestedExperiments binding energy. In mathematical models, additive interac- tions have been shown to yield smoother ,tness surfaces Ourresultsshowdivergencetobepossiblewithmonotonic increasing,tness,whichhintsatthefeasibilityofmonitoring becausetheycanbeoptimizedindependently[7]. similarprocessesinexperiments.Ithasrecentlybeenshown Atahigherlevel,featuresofthenetworktopologyshapethe that the serial dilution assay, as pioneered by Lenski and landscapesurfaceanddivergencepotential.Wefoundthatthe coworkers[30],canbeemployedtoadaptbacterialstrainsto tightlyinterconnectedtopology,aspresentaftertheduplica- a new condition within weeks [19,31]. Similarly, one could tion,doesnotfrustratedivergencebutinsteadpromotesit.In attempt to evolve a duplicate lac repressor/operator copy contrasttoanisolatedrepressor-operatorpair,whereadrop towards the independent regulation of a second operon. in the binding strength decreases the ,tness, the same mutation can be neutral in the interconnected topology. However, this more complex assay does require key mod- Compensation for the decrease in binding strength can be i,cations: (1) growth and selection of the mutants should attributedtotwofeaturesofthetopology.First,thereisthe occur in alternating media, in analogy to our discussion of characteristicpressuretonotbindtherivaloperator:whena multipleinputconditions,and(2)astartingnetworkmustbe PLoSComputationalBiology|www.ploscompbiol.org 0472 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData engineered that satis,es the conditions for DNA-binding pressureforindependentregulationseemstobeadominant divergence: a duplicate repressor/operator and a selective one, as evidenced by the many independent transcription pressurefortightandindependentbinding. factorsthatareparalogs,duplicationsalsohaveyieldedother network motifs. An interesting example is the UxuR/ExuR In practice, one could place the lac operator upstream of theraf,noseutilizationoperon,andconstructalacIduplicate pair of repressors. Like the case studied in this paper, they that is sensitive to raf,nose. This initial situation is now haveoriginatedbyduplicationandsharetwooperators(see similartooursimulations:twolacrepressorsbindtothetwo Protocol S3). However, they seem to have diverged under a lacoperators.Theemployed,tnessde,nitionisalsosuitable: differentselectivepressure,sincetheircrossinteractionwas (1) in media where the two metabolites are both low noteliminated,butinsteadhasbeenretained,formingaso- (supplementede.g.,byanothercarbonsource),themetabolic calledbi-fanmotif[16]. enzymes should not be expressed. The resulting optimal This work describes how regulatory network connections growth is well represented by positive contributions to the can be formed and broken after a duplication event. Our overall,tnessbyhighvaluesfortightbinding.(2)Whenjust quantitative approach takes the selective conditions and onemetaboliteispresent,onescreensforexclusivebinding. molecularadaptabilityexplicitlyintoaccount,andopensupa In a medium without lactose the lactose-sensitive repressor new angle on the duplication-divergence question that is shuts both operons down if binding is still non-exclusive. complementarytoexistingapproaches.Evolutionofnetwork Uponmutationsthatallowthisrepressortobindexclusively connectionsistreatedmoreabstractlyinnumericalstudiesof to the operator of the lactose operon, raf,nose metabolic biological network growth, which have recently received enzymeswouldbeexpressed.Theresultingfastergrowthdue muchattention[10,36,37].Theuseofexperimentaldatawill toraf,noseutilizationthuscorrelateswellwithhighervalues help to perform such studies on a more realistic footing. for exclusive binding. The pressure for a correct behavior Finally, the promising new ,eld of experimental network under multiple conditions prevents the ,xation of trivial engineering [38–40] and evolution (see e.g., [41]) will also solutionsthatwouldjustworkunderonecondition. bene,tfromthequanti,cationofnetworkadaptability. OtherNetworkGrowthScenarios MaterialsandMethods For biological regulatory networks to grow, not only new components are required, but also new and independent Mutational dataset. In thispaper we used anextensive dataset of bindingaf,nitiesoflacrepressorandoperatormutants,obtainedby interactions. Next to the coevolutionary duplication-diver- B.Mu?ller-Hillandcoworkers.Intheseexperiments,repressionvalues gencescenariofornetworkgrowth,alternativemodelsforthe F have been determined in vivo as the ratio of the unrepressed OiRj creationofnewregulatoryinteractionshavebeenproposed. and repressed expression levels of a b-galactosidase (lacZ) reporter gene,controlledbyamutantlacoperatorO andmutantlacrepressor i Inthe,rstalternative,anewoperatormustemergeupstream R.ThiswasdoneusingthestandardassaybyMiller[42].Sincetheb- j oftheregulatedgeneinaneffectivelyrandomDNAsequence galactosidase synthesis is proportional to the fraction of free [32].Thisscenariohasmainlybeenconsideredforeukaryotes operator (seee.g., [43]), we,nd forthe repressionvalue F OiRj ?1þ [R]/K ,whereK istheequilibriumdissociationconstantand[R]is with large upstream regulatory regions and short binding j D D j the concentration of active repressor R. The dataset contains j sites. For longer operators in prokaryotes, this scenario repression values of base pair substitutions leading to changes in requires many neutral mutations before improvements can amino acid residues 1 and 2 of the recognition helix of the lac be selected for (see Protocol S2), which represents a major repressor(Y17andQ18)andbasepairs4and5ofthesymmetriclac evolutionaryobstacle. operator [44].Theseresiduesandbasepairswerefoundto bemost importantforalteringrepressoroperator-bindingaf,nities[14].The Anotherpossiblesourcefornewregulatoryinteractionsis dataset covers a considerable fraction of all possible substitutions lateral gene transfer, which is thought to be the source of involvingahomodimericrepressorandasymmetricoperator(1,286 many paralogs found in prokaryotes [33]. In this scenario outofatotalof6,400).PartofthisrawdataispublishedinLehming divergence would occur while two genes each reside in etal.[14];thefulldatasetisfoundin [21].Thecontributionsofthe different organismal lineages (essentially being orthologs at tworepressoraminoacidstotherepressionvaluewerefoundtobe additive. With this knowledge, repression values could convincingly that stage) and each experiencing different selective con- be assigned to all mutants, including those that were not measured straints. Lateral gene transfer unites the diverged genes, [14].Inthepresentstudyweusetheseassignedrepressionvalues,all resulting in immediate contributions to ,tness by both ofwhicharegivenin[14].Moreover,weextendthedatasettoinclude homologous genes. Although examples of this scenario have heterodimericrepressorsandnon-palindromicoperators(seebelow), been found for enzymes [34], transcription factor-operator to obtain the complete mapping between sequence and repression valuesforallpossiblemutants(13107)inthekeyrepressorresidues interactionsareaspecialcase,asthereisnoobviousinternal Repressionvaluesofheterodimersandnon-palindromicoperators. andoperatorbasepairs. orexternalselectionpressurefortheirinteractiontodiverge Weconsidertherepressorstoactasdimers.Aftertheirduplication, byitself.Ourresultsillustratethefeasibilityofcoevolutionary oncetherepressorsgenesaremutated,thisleadstoheterodimeriza- tionofdistinctmonomers.Whileheterodimerbindingstrengths(F ) He divergence of two transcription factors within a single havenotbeen directlymeasured, theycan be derivedfromthe two organismallineage.These,ndingsaresupportedbythelack corresponding homodimer repression values (F and F ), meas- Ho1 Ho2 ofevidenceforhorizontaltransferofthelacsysteminE.coli ured on a palindromic operator. The heterodimer binding energy DG isthesumofthemonomer-monomerandthedimer-operator He [35].However,thisisnottosaythatlateralgenetransferand binding energy. Simple equilibrium considerations lead to the duplication-divergencearemutuallyexclusive.Summarizing, followingexpression,where[R]inthiscaseisthetotalconcentration the coevolutionary divergence studied here differs from ofrepressorsubunits: alternative models of network growth by providing both a FHe ?1þ?R?2eÀDGHe=kT ?1þ ðF Ho1 Ho2 À1ÞðF À1Þ ð2Þ high probability of selective advantageous point mutations p,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, With this equation, repression values involving non-palindromic andarationaleforadivergencepressure. operators are also automatically taken into account: each dimer- Finally, it is of interest to consider different selective operatorinteractionisbuiltupadditively[22]fromtwointeractions pressures within the same duplication scenario. While the between a monomer and an operator-half. In this derivation the PLoSComputationalBiology|www.ploscompbiol.org 0473 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData dimerization free energy was assumed to be ,xed, since it does not (C)Sequencesforeachround.Mutatedpositionsarecoloredwhite. directlyaffectthespeci,citybywhichtherepressorsrecognizetheir FoundatDOI:10.1371/journal.pcbi.0020058.sg001(349KBPDF). operators. The heterodimer repression value then becomes inde- FigureS2.QualitativeFeaturesofSuccessfullyDivergingPathsinthe pendentofthedimerizationenergy. Optimal pathway simulations. Each repressor monomer is repre- ProbabilisticPathwaySimulations sentedbysixbasepairs(twoaminoacidresidues),andeachoperator Simulations were performed with a 5% growth advantage of a by four base pairs, which are key to speci,c binding. The complete diverged network over the initial duplicate network, and a 5 network with duplicates is thus represented by 20 base pairs. Each populationsizeof10 .Ofalltracedpaths,17%successfullydiverged, simulation run starts with the duplication of a tight binding despitethestrictrequirementsthatpromotetrappinginlocaloptima repressor-operatorpair,havingarepressionvalueof100orhigher. (,tness cannot decrease). Relaxing these conditions would lead to Out of all possible repressor-operator combinations (homodimers largerdivergenceprobabilities. and palindromic operators), there are 132 ful,lling this condition. (A) Histogram showing the number of base mutations until Changingthisthresholddidnotsigni,cantlyaltertheoutcomeofthe divergenceforthesuccessfulpathways. simulations. In order to avoid any bias due to codon usage of the (B)Histogramshowingthelowestrepressionvaluesofeachrepressor startingrepressor,separatesimulationswererunstartingfromeach onitsoperatorduringthesuccessfuldivergencepathways. of its synonymous codon versions. Thesesimulationswereaveraged (C)Histogramshowingthenumberofneutralmutationsthatoccur toproducethepresentedresults.In order to determine the optimal mu tational pathways in the untilthepathwayssuccessfullydiverged. ,tness landscape, an evolutionary algorithm was employed. Begin- FoundatDOI:10.1371/journal.pcbi.0020058.sg002(39KBPDF). ningwithoneofthestartingsequences,eachroundwegeneratedall Protocol S1. Simulation of Mutational Pathways Incorporating mutants that differ by one base pair (60 in total). Of each mutant ProbabilisticPopulationDynamics network, the strength of all eight possible interactions was deter- FoundatDOI:10.1371/journal.pcbi.0020058.sd001(28KBDOC). mined(seeFigure1Bwherefourpossibleinteractionsareschemati- cally shown between the repressor dimers and one of the two Protocol S2. Comment on Neutral Mutations Required for the operators). Interactions between repressor homodimers and palin- EmergenceofaNewOperator dromic operators were directly assigned from the published FoundatDOI:10.1371/journal.pcbi.0020058.sd002(20KBDOC). repression values [14]. Other interactions were calculated from the measured data as described above. Next, we selected the best L Protocol S3. Alternative Selective Pressures and the Escherichia coli networks to the next round based on a ,tness parameter that is RegulatoryNetwork directly calculated from the interaction strengths (see equation 1). FoundatDOI:10.1371/journal.pcbi.0020058.sd003(333KBDOC). The next round started by again generating all single base pair 5 mutants of the L selected networks. The effect of L was assessedby varying it between 1 and 10 . Decreasing ,tness steps were not allowed,andincaseofequal,tness,parentswererankedabovetheir Acknowledgments offspring.Theserulesmakedivergenceharderbecausetheyconstrain thespacethatcanbeexplored.Theevolutionarycyclewasrepeated We thank Benno Mu?ller-Hill for kindly providing additional data, until the ,tness could not be further improved. Pathways were James Shapiro for stimulating initial discussions, and Dennis Bray, consideredtobesuccessfulwhenthe,tnesscamewithinafactor10 Harald Tepper, Marileen Dogterom, Pieter-Rein ten Wolde, Kobus ofthehighest,tnessinthelandscape. Kuipers, and Sorin Tanase-Nicola for critical reading of the manu- script. Author contributions. FJP conceived and designed the model, SupportingInformation analyzedthedata,andwrotethepaper.DJKdesignedandperformed FigureS1.TypicalDivergencePathway,withtheAdditionalRequire- the computational experiments, analyzed the data, and wrote the mentofContinuedTightBindingofBothRepressors(F100) paper. SJT conceived and designed the model, designed and performed the computational experiments, analyzed the data, and (A) Evolving interaction network, where line thickness denotes wrotethepaper. binding strength between repressor monomer and operator-half. Funding.ThisworkispartoftheresearchprogramoftheStichting Dotted lines denote negligible repression. Yellow crosses indicate voor Fundamenteel Onderzoek der Materie, which is ,nancially repressor and operator mutations, which are positioned at the top andbottomoftheinteractionlinesrespectively. supported by the Nederlandse Organisatie voor Wetenschappelijke (B) Fitness trajectory (black) and corresponding repression of each Onderzoek. repressoronitsoperator(redandblue).Fitnessisnormalizedtothe Competinginterests.Theauthorshavedeclaredthatnocompeting & maximumvalue(;131010). interestsexist. References of the lac repressor. XIV. Analysis of 4,000 altered Escherichia coli lac repressorsrevealsessentialandnon-essentialresidues,aswellas‘‘spacers’’ 1. Stephens SG (1951) Possible signi,cance of duplication in evolution. Adv Genet4:247–265. whichdonotrequireaspeci,csequence.JMolBiol240:421–433. 2. OhnoS(1970)Evolutionbygeneduplication.NewYork:Springer-Verlag. 13. Mu?ller-HillB(1996)Thelacoperon:Ashorthistoryofageneticparadigm. 160p. Berlin:DeGruyter.207p. 3. Lynch M, Conery JS (2000) The evolutionary fate and consequences of 14. LehmingN,SartoriusJ,Kisters-WoikeB,vonWilcken-BergmannB,Mu?ller- duplicategenes.Science290:1151–1155. Hill B(1990) Mutantlac repressors withnewspeci,citieshint atrulesfor 4. TeichmannSA,MadanBabuM(2004)Generegulatorynetworkgrowthby protein-DNArecognition.EMBOJ9:615–621. duplication.NatGenet36:492–496. 15. Kimura M (1962) On the probability of ,xation of mutant genes in a 5. Madan Babu M, Teichmann SA (2003) Evolution of transcription factors population.Genetics47:713–719. and the gene regulatory network in Escherichia coli. Nucleic Acids Res 31: 16. Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, et al. (2002) 1234–1244. Networkmotifs:Simplebuildingblocksofcomplexnetworks.Science298: 6. Nguyen CC, Saier JMH (1995) Phylogenetic, structural, and functional 824–827. analyses of the LacI-GalR family of bacterial transcription factors. FEBS 17. HughesAL(1994)Theevolutionoffunctionallynovelproteinsaftergene Lett377:98–102. duplication.ProcRoySocLondonB256:119–124. 7. Kauffman SA(1993) Theoriginsoforder:Self-organizationandselection 18. KochAL(1983)Theproteinburdenoflacoperonproducts.JMolEvol19: inevolution.Oxford:OxfordUniversityPress.734p. 455–462. 8. Bray D, Lay S (1994) Computer simulated evolution of a network of cell- 19. Dekel E, Alon U (2005) Optimality and evolutionary tuning of the signalingmolecules.BiophysJ66:972–977. expressionlevelofaprotein.Nature436:588–592. 9. Franc?ois P, Hakim V (2004) Design of genetic networks with speci,ed 20. Fukami-Kobayashi K, Yoshio Tateno Y, Nishikawa K (2003) Parallel functionsbyevolutioninsilico.ProcNatlAcadSciUSA101:580–585. evolution of ligand speci,city between lacI/galR family repressors and 10. Barabasi AL, Oltvai ZN (2004) Network biology: Understanding the cell’s periplasmicsugar-bindingproteins.MolBiolEvol20:267–277. functionalorganization.NatRevGenet5:101–113. 21. Lehming N (1990) Regeln fu?r Protein/DNA-Erkennung [PhD thesis]. 11. SenguptaAM,DjordjevicM,ShraimanBI(2002)Speci,cityandrobustness intranscriptioncontrolnetworks.ProcNatlAcadSciUSA99:2072–2077. Universita?tzuKo?ln.99p. 12. MarkiewiczP,KleinaLG,CruzC,EhretS,MillerJH(1994)Geneticstudies 22. HollisM,ValenzuelaD,PioliD,WhartonR,PtashneM(1988)Arepressor PLoSComputationalBiology|www.ploscompbiol.org 0474 May2006|Volume2|Issue5|e58 EvolutionaryPathwaysfromMutationData heterodimer binds to a chimeric operator. Proc Natl Acad Sci U S A 85: genomic repertoires in bacteria. PLoS Biol 3: e130. DOI: 10.1371/journal. 5834–5838. pbio.0030130 34. InagakiY,DoolittleWF,BaldaufSL,RogerAJ(2002)Lateraltransferofan 23. MacArthurS,Brook,eldJFY(2004)Expectedratesandmodesofevolution EF-1agene:OriginandevolutionofthelargesubunitofATPsulfurylasein ofenhancersequences.MolBiolEvol21:1064–1073. eubacteria.CurrBiol12:772–776. 24. ConantGC,WagnerA(2003)Asymmetricsequencedivergenceofduplicate 35. Stoebel DM (2005) Lack of evidence for horizontal transfer of the lac genes.GenomeRes13:2052–2058. operonintoEscherichiacoli.MolBiolEvol22:683–690. 25. LynchM(2005)Simpleevolutionarypathwaystocomplexproteins.Protein 36. Bhan A, Galas DJ, Dewey TG (2002) A duplication growth model of gene Sci14:2217–2225. expressionnetworks.Bioinformatics18:1486–1493. 26. PilarFrancinoM(2005)Anadaptiveradiationmodelfortheoriginofnew 37. Wagner A (2002) How the global structure of protein networks evolves. genefunctions.NatGenet37:573–577. ProcRoySocLondonB270:457–466. 27. Ju?rgens C, Strom A, Wegener D, Hettwer S, Wilmanns M, et al. (2000) 38. Kobayashi H, Kaern M, Araki M, Chung K, Gardner TS, et al. (2004) Directedevolutionofa(ba)8-barrelenzymetocatalyzerelatedreactionsin Programmable cells: Interfacing natural and engineered gene networks. twodifferentmetabolicpathways.ProcNatlAcadSciUSA97:9925–9930. ProcNatlAcadSciUSA101:8414–8419. 28. O’BrienPJ,HerschlagD(1999)Catalyticpromiscuityandtheevolutionof 39. Weber W, Fussenegger M (2002) Arti,cial mammalian gene regulation newenzymaticactivities.ChemBiol6:R91–R105. networks: Novel approaches for gene therapy and bioengineering. J 29. Kondrashov FA, Rogozin IB, Wolf YI, Koonin EV (2002) Selection in the Biotech98:161–187. evolutionofgeneduplications.GenomeBiol3:research0008. 40. Farmer WR, Liao JC (2000) Improving lycopene production in Escherichia 30. LenskiRE,TravisanoM(1994)Dynamicsofadaptationanddiversi,cation: colibyengineeringmetaboliccontrol.NatBiotech18:533–537. A 10,000-generation experiment with bacterial populations. Proc Natl 41. YokobayashiY,WeissR,ArnoldFH(2002)Directedevolutionofagenetic AcadSciUSA91:6808–6814. circuit.ProcNatlAcadSciUSA99:16587–16591. 31. Ibarra RU, Edwards JS, Palsson BØ (2002) Escherichia coli K-12 undergoes 42. MillerJH(1972)Experimentsinmoleculargenetics.NewYork:ColdSpring adaptive evolution to achieve in silico predicted optimal growth. Nature HarborLaboratoryPress.466p. 420:186–189. 43. SadlerJR,NovickR(1965)Thepropertiesofrepressorandthekineticsof 32. Berg J, Willmann S, La?ssig M (2004) Adaptive evolution of transcription itsaction.JMolBiol12:305–327. factorbindingsites.BMCEvolBiol4:42. 44. Sadler JR, Sasmor H, Betz JL (1983) A perfectly symmetric lac operator 33. Lerat E, Daubin V, Ochman H, Moran NA (2005) Evolutionary origins of bindsthelacrepressorverytightly.ProcNatlAcadSciUSA80:6785–6789. PLoSComputationalBiology|www.ploscompbiol.org 0475 May2006|Volume2|Issue5|e58