Table Of ContentPurifying Selection Determines the Short-Term Time
Dependency of Evolutionary Rates in SARS-CoV-2 and pH1N1
Influenza
Mahan Ghafari,1 Louis du Plessis,1 Jayna Raghwani ,1 Samir Bhatt,2 Bo Xu,3 Oliver G. Pybus ,1 and
Aris Katzourakis*,1
1DepartmentofZoology,UniversityofOxford,Oxford,UnitedKingdom
2MRCCentreforGlobalInfectiousDiseaseAnalysis,JameelInstituteforDiseaseandEmergencyAnalytics,ImperialCollegeLondon,
London,UnitedKingdom
3DepartmentofEarthSystemScience,TsinghuaUniversity,Beijing,China
*Correspondingauthor:E-mail:[email protected]. D
o
Associateeditor:HarmitMalik w
n
lo
a
d
e
Abstract d
fro
High-throughput sequencing enables rapid genome sequencing during infectious disease outbreaks and provides an m
h
opportunity to quantify the evolutionary dynamics of pathogens in near real-time. One difficulty of undertaking evo- ttp
s
lutionaryanalysesovershorttimescalesisthedependencyoftheinferredevolutionaryparametersonthetimespanof ://a
observation.Crucially,thereareanincreasingnumberofmolecularclockanalysesusingexternalevolutionaryratepriors ca
d
to infer evolutionary parameters. However, it is not clear which rate prior is appropriate for a given time window of e
m
observation due to the time-dependent nature of evolutionary rate estimates. Here, we characterize the molecular ic
.o
evolutionarydynamicsofSARS-CoV-2and2009pandemicH1N1(pH1N1)influenzaduringthefirst12monthsoftheir u
p
respective pandemics. We use Bayesian phylogenetic methods to estimate the dates of emergence, evolutionary rates, .co
m
andgrowthratesofSARS-CoV-2andpH1N1overtimeandinvestigatehowvaryingsamplingwindowanddatasetsizes /m
affect the accuracy of parameter estimation. We further use a generalized McDonald–Kreitman test to estimate the be
numberofsegregatingnonneutralsitesovertime.Wefindthattheinferredevolutionaryparametersforbothpandemics /artic
aretimedependent,andthattheinferredratesofSARS-CoV-2andpH1N1declineby(cid:2)50%and(cid:2)100%,respectively, le
/3
over the course of 1 year. After at least 4months since the start of sequence sampling, inferred growth rates and 9
/2
emergence dates remain relatively stable and can be inferred reliably using a logistic growth coalescent model. We /m
showthatthetimedependencyofthemeansubstitutionrateisduetoelevatedsubstitutionratesatterminalbranches sa
c
whichare2–4timeshigherthanthoseofinternalbranchesforbothviruses.Theelevatedrateatterminalbranchesis 00
9
strongly correlated with an increasing number of segregating nonneutral sites, demonstrating the role of purifying /6
A5
selectioningeneratingthetimedependencyofevolutionaryparametersduringpandemics. 09
r5
Keywords:substitutionrate,molecularclock,clockrate,purifyingselection. t23
i b
cy
l g
eue
s
Introduction t o
n
ofcoalescentpopulationmodels,changesinnaturalselection, 1
4
Rapidwhole-genomesequencinghasbecomepartofpatho- differing replication rates in various host reservoirs, or se- F
e
gensurveillance systems and is important to both infection quencingerrors),misspecificationofclockmodelsandsatu- b
ru
controlandenablesadetailedinvestigationoftheepidemi- ration of nucleotide changes can result in rate a
ry
ological and evolutionary characteristics of pathogens. underestimation (Clark and Whittam 1992; Sullivan and 2
0
Quantifyinginfectiousdiseaseevolutionenablestheinference Joyce 2005; Duch^ene et al. 2014; Ghafari et al. 2021). This 22
ofparameterssuchastimesoforigin,epidemicgrowthrates, nonstationarity can render findings that are confusing or
andevolutionaryrates(Fraseretal.2009;Smithetal.2009;Lu conflicting such that an estimated substitution rate over
etal.2020). onetimewindowisnottransferabletotheanalysisofanother
Oneoftheperilsofmakingsuchinferencesovershorttime (Ho et al. 2011). For example, during the 2014–2016 West
periods(i.e.,afewmonthsoryears)isthattheinferredparam- Africa Ebolavirus epidemic, there was controversy and con-
eters of interest may vary over the timespan of observation cernregardingthevirus’“mutationrate,” becauseearlyesti-
(Meyeretal.2015).Althoughseveralfactorscanleadtothe matesofthesubstitutionratefromtheepidemicappearedto
misestimationofsubstitutionrates(suchasmisidentification be approximately twice the average rate across multiple
(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(2):msac009 doi:10.1093/molbev/msac009 AdvanceAccesspublicationJanuary17,2022 1
MBE
Ghafarietal. . doi:10.1093/molbev/msac009
outbreaks (Gire et al. 2014; Holmes et al. 2016). Confusion molecularclockrate,timeoforigin,growthrate,andnumber
suchasthisisalsocommonintheSARS-CoV-2literatureand ofnonneutralsitesfordifferentdatasetsthatrepresentdif-
arises in part due to an incomplete understanding of esti- ferenttimescalesofgenomicobservation.Weshowthattime
mated substitution rates by users of phylogenetic analysis dependency in estimates of the mean substitution rate is
software. Given the importance of phylogenetic dating and dominated by elevated rates at terminal branches, whereas
clock estimation for all SARS-CoV-2 genomic epidemiology theestimatedrateofsubstitutionatinternalbranchesdoes
worldwide,understandingthepatternofinferredsubstitution notexhibitatime-dependentdecay.Wethenuseageneral-
rate over time and investigating the potential underlying ized McDonald–Kreitman test based on nucleotide site fre-
mechanismsinvolvedincreatingtime-dependentrateeffects quenciesthatallowspurifyingselectiontobequantifiedover
in viruses can shed light on the molecular evolutionary dy- shorttimescalesanddemonstratethatthereisastrongcor-
namicsofviruses. relationbetweentheelevatedratesatterminalbranchesand
Another major obstacle is that during an ongoing out- the high number of low frequency nonneutral sites in both
break, the level of sequence diversity may be so low that SARS-CoV-2andpH1N1genomes. D
o
theamountofaccruedevolutionarychangeisnotsufficient w
n
to make informative phylogenetic inference possible Results lo
a
d
(Duch^eneetal.2020).Several statisticaltestshavebeende- e
We use Bayesian phylogenetic methods implemented in d
velopedtoensureenoughtemporalsignalispresentinaset BEAST v.1.10 to estimate the molecular clock rate, times of fro
m
of temporally sampled genome sequences to reliably infer
origins, and growth rates of the SARS-CoV-2 and pH1N1 h
evolutionary parameters(Duch^eneet al.2015; Murrayetal. pandemics. The data sets varied in size and temporal sam- ttps
2016;Duch^eneSandDuch^eneDA2020). plingrangeandwerecreatedtoreflectthewayinwhichreal ://a
Morespecifically,thereisextensiveevidencethatpurifying datasetsaccumulateinsizeanddiversityduringthecourseof cad
selectioninvirusescanresultinvariedestimatesforratesof e
anepidemic.Specifically,weaggregateallavailablesamplesup m
substitutionandevolutionaryrateratio,dN/dS,acrossdiffer- ic
toeachrespectivemonthandinfertheparametersofinterest .o
ent timescales (Sharpet al. 2001; Holmes2003; Hughesand u
Hughes 2007). It is likely that a greater proportion of poly- (fig. 1 and supplementary fig. S1, Supplementary Material p.co
online). We first compare the results for SARS-CoV-2 using m
morphisms observed among genomic sequences sampled twocoalescentgrowthpriors:exponentialgrowthandlogistic /m
early in an epidemic are segregating deleterious mutations, be
growth. Although this is not a comprehensive comparison /a
whichwillonlypersistforalimitedtimebeforebeingelimi- betweenallphylogeneticmodelsthatcanbeusedtoexplain rtic
natedbypurifyingselection(Lametal.2016).Thedurationof le
the evolutionary dynamics, it allows us to find the better- /3
thistime-dependenteffectmaybeprolongedduetoincom- 9
fittingmodelthatintroduceslesserrorinparameterestima- /2
pletepurifyingselectioninpopulationswithverylargeeffec- /m
tion—the demographic model is effectively a nuisance vari- s
tivepopulationsizes(Woodhams2006). a
Althoughevidenceofstrongpurifyingselectionhasmostly able.Ouranalysissuggests,exceptforthemonthofJanuary, c00
thelogisticgrowthcoalescenttreepriorisabetterfittothe 9
been reported at the species level and over very long time- /6
scales (i.e., thousands to millions of years), using inference data (table 1). As reported in previous studies, we find that 509
there is not enough temporal signal in the data to reliably 5
methods based on the dN/dS ratio (Ho et al. 2005; Pybus 2
3
etal.2007;WertheimandPond2011),therehavebeenfew infertheevolutionaryparametersduringthefirst2monthsof by
theSARS-CoV-2pandemic(Duch^eneetal.2020).Thisresults g
studies at the intra-population level and over short time- u
sscuaffliecsi,emntalyinllayrbgeecnauumsebtheresreoifssoafmtepnlensotohpropuogrthuntiimtyetotocotrllaecckt isntatthiseticuanlduenrceesrttiaminattyiofnorotfhesupbasrtaitmuteitoenrersattiemsaatsesw.Tehllealsachkigohf est on
low-frequencyvariants(Hedgeetal.2013;Meyeretal.2015). temporal signal in the SARS-CoV-2 samples is also evident 14 F
Crucially, althoughpurifyingselectionhasbeenput forward from the number of new singletons (i.e., single nucleotide eb
asapossibleexplanationforthetimedependencyofsubsti- variations in the data set compared with the ancestral se- rua
tution rates over such timescales, there have not been any quence) added the data set per month during the first 2– ry 2
systematicstudiestoinvestigatetheroleofpurifyingselection 3months(supplementaryfig.S2,SupplementaryMaterialon- 02
2
and quantify its contributionto altering the inferredsubsti- line). Ontheother hand, forpandemicH1N1(pH1N1),the
tutionrateofvirusesovertime.Furthermore,usingstandard molecularclockrateisuptotwotimeshigherforthefirst3
phylogeneticmethodstocomputethedN/dSratioforcon- months than for the following months (fig. 1A and B).
specific sequences sampled from a single population over AlthoughtheinferredsubstitutionrateofSARS-CoV-2tends
shorttimescalesmaybeinappropriateasthedifferencesbe- todecreaseasweincreasethetimespanofmeasurement,the
tweensequencesoversuchtimescalesrepresentsegregating rateforpH1N1doesnotchangeconsiderablyafterthefirst3–
polymorphismsasopposedtofixedsubstitutionsalonginde- 4monthsofmeasurement,inagreementwithpreviousfind-
pendent lineages. The former has been shown to produce ings(Meyeretal.2015).
verydifferentestimatesofthedN/dSratioovertime(Rocha BoththeinferredtimesoforiginandgrowthratesofSARS-
etal.2006;KryazhimskiyandPlotkin2008). CoV-2andpH1N1inearlymonthshavewidecredibleinter-
Inthisstudy,weaimtoidentifyandquantifythesourceof valsduetothehighuncertaintyassociatedwithsmallsample
time dependency of virus substitution rate estimates over sizes and narrow sampling windows. The precision of these
short time periods and characterize estimates of the twoinferredparametersincreaseswiththeadditionofdata
2
MBE
EvolutionaryRatesinSARS-CoV-2andpH1N1Influenza . doi:10.1093/molbev/msac009
A B
D
o
w
n
lo
a
C D 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
E F le
/3
9
/2
/m
s
a
c
0
0
9
/6
5
0
9
5
2
3
b
y
g
u
e
s
t o
n
1
4
F
e
b
ru
a
ry
2
FIG.1.(A,B)Inferredrates,(C,D)timesoforigin,and(E,F)growthratesofSARS-CoV-2(leftcolumn)andpH1N1influenza(rightcolumn)usingan 0
2
exponential(red)andalogistic(black)growthmodel.OpencirclesrepresentnonconvergenceforatleastoneparameterintheBayesiananalysis. 2
Notethaty-axesarenotonthesamescaleforSARS-CoV-2andpH1N1.
fromlatermonthsandremainsroughlyconsistentafterthe expectation after 6months of sampling and yield unreliable
first4–5monthsofmeasurement(fig.1C–F).Usingthelogis- estimates(fig.1CandD).Thisislikelybecauseafterthefirst
tic growth model, we find the estimated times of origin for few months, the growth rate declines and the population
SARS-CoV-2andpH1N1tobeOctober28,2019(95%HPD: dynamics deviate from the exponential growth model
September 30, 2019, November 24, 2019) and January 18, (fig. 1E and F). This, in turn, results in the underestimation
2009(95% HPD: December 14, 2008, February 22, 2009), re- ofgrowthrateofSARS-CoV-2andinferringanoldertimeof
spectively, which is also in agreement with previous studies origin. We also note that, from a molecular epidemiology
(Smithetal.2009;Hedgeetal.2013;Luetal.2020).However, perspective,usingtheexponentialgrowthcoalescentmodel
under the exponential growth model, the inferred time of wouldnotbeappropriatenorrealisticparticularlyoverlonger
origin ofSARS-CoV-2 samples significantly diverge from the timescales (i.e., >6months) as the pandemic dynamics did
3
MBE
Ghafarietal. . doi:10.1093/molbev/msac009
Table1. Log-MarginalLikelihoodsofExponentialandLogisticGrowthModelswithIncreasingTemporalRangesofSamplingDates.
MonthofSampling NumberofSamples Log-MarginalLikelihood Log-MarginalLikelihood BayesFactor
ExponentialGrowthModel LogisticGrowthModel
Jan 41 241,022.77 241,025.56 22.79
Feb 121 243,062.00 243,047.11 114.88
Mar 181 244,627.77 244,607.95 119.78
Apr 241 246,524.75 246,506.25 118.53
May 301 248,587.39 248,548.15 139.24
Jun 361 252,310.61 252,248.23 162.38
Jul 421 255,138.43 255,036.60 1101.83
Aug 481 258,580.94 258,430.09 1150.85
Sep 541 262,101.22 261,805.93 1295.29
Oct 601 265,433.20 265,005.09 1428.11
Nov 661 269,121.03 268,921.98 1199.05 D
Dec 721 273,094.26 272,389.69 1704.57 ow
n
NOTE.—Takingexponentialgrowthasthenullmodel,weselectthelogisticgrowthmodelforanydatasetwithapositiveBayesfactor. loa
d
e
d
fro
not continue to grow exponentially (i.e., the growth rate replacement sites is always greater than silent sites in m
started to drop after the first few months due to various ORF1abwhereas,forHA,silentsitesareinmajorityinmost http
behavioral changes in the population and the implementa- months(fig.2CandD).Thetotalnumberofnonneutralsites s://a
tionofnonpharmaceuticalinterventions).Despitethediffer- in ORF1ab is higher because it is a much longer gene c
a
ence in the estimated time of origin, both the exponential ((cid:2)21,000nt)comparedwithHA((cid:2)1,800nt). de
m
and logistic growth models infer very similar clock rates for Oneoftheimpactsofpurifyingselectionovershorttime- ic
SARS-CoV-2overtime. scales is that the number of replacements on terminal .ou
p
Purifyingselectionhasbeenoftencitedasoneofthemain branches should be higher than that on internal branches. .c
o
evolutionary processes contributing to the elevation of esti- To measure this effect during the SARS-CoV-2 and pH1N1 m
/m
mated molecular evolutionary rates over short timescales epidemics,weusedabranch-specifictwo-parametermolec- b
e
(Hedgeetal.2013).Theargumentisthatlow-frequencydel- ularclockmodelinBEASTandinferseparateratesofsubsti- /a
eterious mutations can segregate in a population for some tutionforterminalandinternalbranches(seefig.2EandF). rticle
time before being purged because of purifying selection. Ourresultsshowthatforbothviruses,theratesofsubstitu- /3
9
Therefore,theproportionofallchangesthataredeleterious tiononterminalbranchesare2–4timeshigherthanonin- /2/m
ishighwhenthephylogenetictreeisshort,andlowerwhen ternal branches. The difference is more dramatic in SARS- s
a
c
thetreeislongerasittakesmoretimeforsharedneutralor CoV-2wheretheaverageratioofsubstitutionratesattermi- 0
0
9
advantageouschangestoaccruewithrespecttotheancestral nal branches are four times higher than that of internal /6
5
state.Toevaluatethishypothesis,weinvestigatethecorrela- branches whereas it is only two times higher in HA during 0
9
5
tionbetweenthenumberofnonneutralpolymorphismsand thefirst12monthofobservation.Figure2EandFalsoshows 2
3
estimatedsubstitutionratesovertime.Ourresultsshowthat thatthemeanandcredibleintervalsoftheinferredsubstitu- b
y
thereisahigherproportionofsegregatingmutationsatlow tion rate at terminal and internal branches are nonoverlap- g
u
e
frequencies(<15%)duringthefirstfewmonthsofobserva- ping during several months of rate measurement which s
t o
tion, in both the ORF1ab of SARS-CoV-2 and the HA of demonstrates that there is a significant difference between n
1
pH1N1influenza(fig.2AandB).Inparticular,theproportion therateofevolutionatthesetwobranchcategories.Thisalso 4
F
oflow-frequencynonneutralsitesinORF1abwashigherdur- agrees with the observation of a higher proportion of non- e
b
ingthefirst4–5monthsofobservationandhadaroughly4- neutralsitesinSARS-CoV-2samplescomparedwithpH1N1 ru
a
folddropduringthatperiodcomparedwitha2-folddropin because there is going to be more deleterious mutations at ry 2
0
HA gene. After the first 5months, the proportion of low terminal branches of SARS-CoV-2. Further, there is a clear 2
2
frequency nonneutral sites remains roughly constant at decline in the estimated rate at terminal branches as the
around 15% in ORF1ab, whereas it shows an uptick from sampling window increases. In contrast, the estimated rate
15%to17%inHAtowardthebeginningof2010.Theslight at internal branches gradually increases through time (with
increasein theproportion ofnonneutral sites inHA during theexceptionofthefirst2–3monthsforpH1N1).
this period is also in agreement with a rise in the relative We also investigate the role of sequencing error as a po-
genetic diversity of pH1N1 around the world in late 2009/ tentialconfounderfortheelevatedratesatterminalbranches.
early2010whichmayalsohaveresultedinanincreaseinthe Bymaskingmorethan110sitesofSARS-CoV-2,inadditionto
numberofsegregatingdeleteriousmutationsinthepopula- our standard filtering step (see Materials and Methods sec-
tion(Suetal.2015).Wenotethatalthoughitisunlikelyfor tion),whicharesuggestedtobepronetorecurrentsequenc-
mutationsinthelowfrequencyclasstohavestronglypositive ingerrorandappeartobehighlyhomoplastic,werecalculate
fitness effects, they may also contain neutral or adaptive therateofevolutionatterminalandinternalbranches(sup-
mutationsthathavenotreachedsufficientlyhighfrequencies plementaryfig.S3,SupplementaryMaterialonline).Wefind
yet. We also find that the number of low frequency that although the overall inferred rates drop for the
4
MBE
EvolutionaryRatesinSARS-CoV-2andpH1N1Influenza . doi:10.1093/molbev/msac009
A B
non-neutralsites non-neutralsites
fonoitcarF Frac(cid:2)on of
D
C D ow
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
E F o
m
)YSS SSY) /mbe
(etarnoitu u(cid:2)on rate( /article/39
ttsbi bs(cid:2)t /2/m
u u
s s s
d d a
e e c
rre err 00
fnI Inf 9/6
5
0
9
5
2
3
b
G H y g
u
e
eutralsites utralsites st on 14 F
n e e
ononfnoitc- n(cid:2)on ofnon- bruary 202
arF Frac 2
Subs(cid:2)tu(cid:2)on rateatterminalbranches Subs(cid:2)tu(cid:2)on rateatterminalbranches
FIG.2. (A,B)NumberofnonneutralsitesovertimeforSARS-CoV-2(leftcolumn)andpH1N1(rightcolumn).(C,D)Numberofreplacement
(dashedblueline),silent(solidblueline)sites,andtheirratio(redline)overtime.(E,F)Meanclockrate(black),andtheratesattheterminal(red)
andinternal(blue)branches.TheMCMCchainsforthefirstmonthofsamplingSARS-CoV-2andthefirstthreemonthsofpH1N1donotconverge
usingthelogisticgrowthcoalescentmodel.Instead,theexponentialgrowthcoalescentmodelwasused.(G,H)Correlationcoefficientbetweenthe
substitution rate at terminal branches and number of nonneutral sites—excluding the estimates from the first month of sampling due to
inadequatetemporalsignalandsignificantuncertaintyintheinferredparameters.
5
MBE
Ghafarietal. . doi:10.1093/molbev/msac009
alignments with masked sites (due to their lower genetic Ourresultsprovidedfurtherevidencethatshort-termrate
variation compared with unmasked alignments), the time- estimates are subject to time-dependent rate effects largely
dependent rate drop at terminal branches is retained (sup- due to transient polymorphisms (Hedgeet al. 2013; Meyer
plementary fig. S3B, Supplementary Material online). et al. 2015). In addition, our estimated mean substitution
Therefore,sequencingerrorcannotbetheunderlyingsource rates, times of origin, and growth rates agree with previous
fortimedependencyofevolutionaryrates.Infigure2GandH, studies of SARS-CoV-2 (Lu et al. 2020; Volz et al. 2020;
wefurthershowthatthereisanoverallsignificantcorrelation Duch^eneetal.2020;duPlessisetal.2021)andpH1N1influ-
betweentherateofsubstitutionatterminalbranchesforthe enza (Rambaut and Holmes 2009; Smith et al. 2009; Hedge
unmasked alignments and number of nonneutral sites in et al. 2013; Meyer et al. 2015). We found that although the
bothSARS-Cov-2(R2¼0.74andP<0.001)andpH1N1(R2 time-dependent rate effect for pH1N1 is less pronounced
¼ 0.53 and P<0.01) (fig. 2G and H). By examining the site after the first 3months of sampling, the same effect in
frequencyspectrumforbothdatasets,wecanalsoseethat SARS-CoV-2 is more pronounced and continues even after
mostofthevariationcomesfromthelow-frequencyregime 12months of sampling. This may be the result of a much D
o
w
(<15%) with very limited variation in the mid- to high- higher mean substitution rate in pH1N1 compared with n
lo
frequency regimes (supplementary fig. S2, Supplementary SARS-CoV-2. Wealso foundahigherproportionoflowfre- a
d
Materialonline). quencynonneutralsitesfortheSARS-CoV-2datasetduring ed
Furthermore,inthepH1N1dataset, wecanseeadipin the first 8months of observations and a higher ratio of re- fro
m
the inferred rate at terminal branches in December 2009 placementtosilentsitesduringtheentire1yearofobserva- h
(fig.2F).Thereisalsoasuddenlossoflowfrequencygenetic tioncomparedwithpH1N1.Thesamemethodusedinthis ttp
s
diversity in the data set during the same period which is studytoinferthenumberoflowfrequencynonneutralsites ://a
c
partiallyrecoveredinFebruaryandMarch(seethetrajectory hasalsobeenusedtocharacterizethemolecularevolutionary ad
e
ofnewlyaddedsingletonsandlowfrequencyvariantsinsup- dynamicsofotherRNAvirusesincludingseasonalcoronavi- m
ic
plementary fig. S2D and F, Supplementary Material online). rusesandemergingavianinfluenzaviruses(Bhattetal.2013; .o
u
Ontheotherhand,thepatternofsubstitutionrateatinternal Luetal.2018;KistlerandBedford2021). p.c
o
branchesisnonmonotonicovertime.Inparticular,weseean By dividing the data into successively longer temporal m
increaseintherateatinternalbranchesofSARS-CoV-2from intervals,weshowedthattheoverabundanceofdeleterious /m
b
e
Augustonward.Similarly,thereisasteadyriseintherateat mutationsatterminalbranchesisthemainreasonbehindthe /a
internalbranchesofpH1N1fromFebruaryonward. gradualdecay inmeansubstitutionrateovertimeinSARS- rtic
le
CoV-2 and pH1N1 influenza. Further work can be done to /3
Discussion 9
quantifythelonger-termeffectofpurifyingselectiononthe /2
/m
Wefoundthattheoverabundanceofdeleteriousmutations timedependencyofsubstitutionrates. s
a
duringtheearlystagesofboththeSARS-CoV-2andpH1N1 c0
Materials and Methods 0
influenzapandemicsstronglycorrelateswithhighersubstitu- 9
/6
tionratesatterminalbranchesrelativetointernalbranchesof WedownloadedallSARS-CoV-2sequencesfromGISAIDand 50
9
inferred phylogenetic trees. Once there is enough temporal pH1N1influenzasequencesfromGenBankandalignedthem 5
2
3
signalinthedatatoreliablyestimaterelevantepidemiological usingMUSCLEv3.8.425(Edgar2004)—acompletemetadata b
y
and evolutionary parameters, the mean substitution rates tableacknowledgingtheauthors,originatingandsubmitting g
u
decline over the course of 1 year of rate measurement, laboratoriesoftheSARS-CoV-2sequencedataisavailablein es
whereas the estimated time of origin and growth rate are supplementarytableS1,SupplementaryMaterialonline.For t on
more stable and remain roughly the same after the first 4– the SARS-CoV-2 data set, we mask the first 54 and last 240 14
F
5monthsofmeasurement.Wefoundthatthisdecliningpat- sitesofthealignmentsandonlyincludecompletesequences e
b
terninmeansubstitutionrateiscausedbytheratedecayat thataremorethan29,000ntlongwithhighcoverageasde- ru
a
terminalbrancheswhereastherateatinternalbranchesdoes terminedbyGISAID’sdefaultsearchoption(i.e.,entrieswith ry
2
notexhibitsuchapatternovertime.Inparticular,thegradual <1% Ns and <0.05% unique amino acid substitutions). 02
2
increaseinthesubstitutionrateatinternalbranchesoverthe Furthermore, for part of the analysis where we investigate
span of several months for both SARS-CoV-2 and pH1N1 the role of sequencing error in time-dependent rate effects
influenzacouldberelatedtosignaturesofadaptiveevolution. atterminalbranchesforSARS-CoV-2,wemaskanadditional
Wealsofoundthatsitesthatarepronetorecurrentsequenc- 114 sites that appear to be highly homoplastic or prone to
ingerrorandappeartobehighlyhomoplasticdonotmakea recurrentsequencingerror(DeMaioetal.2020,https://viro-
significant contribution to the time-dependent pattern of logical.org/t/issues-with-sars-cov-2-sequencing-data/473).We
substitutionrates.Althoughbothsequencingerroranddel- specificallyselectpH1N1sequencesfromGenBankforthefull
eterious mutations appear as singletons or low-frequency coding region of the hemagglutinin (HA) and ensure that
variationsinthedata,deleteriousmutationsarestronglybi- none of the SARS-CoV-2 and pH1N1 samples have under-
ased toward changing amino acids, whereas sequencing goneserialpassaging.
errors are equally likely to be synonymous or nonsynony- To investigate the effects of the temporal range of sam-
mous.Therefore,theydonotmakesimilarcontributionsto plingdatesontheaccuracyofparameterestimation,wein-
alteringtheinferredsubstitutionratesovertime. crementally increase the size of each data set by adding 60
6
MBE
EvolutionaryRatesinSARS-CoV-2andpH1N1Influenza . doi:10.1093/molbev/msac009
genomes, chosen randomly, for every additional month of EstimatingtheSiteFrequencyandNumberof
sampling. Thus, after 12months of sampling since the first NonneutralSites
sequencewasuploadedonGISAIDand/orGenBank,wehave We use the adapt-a-rate package (Bhatt et al. 2010, 2011;
720samples.Wenotethatduetoalackoftemporalsignalin Raghwani et al. 2016), a generalized McDonald–Kreitman
theearlySARS-CoV-2samples(i.e.,mostoftheearlysamples test,toestimatethenumberofnonneutralsitesbyassuming
werealmostcompletelyidentical)andfailureoftheMarkov that deleterious mutations are mostly confined to the low
chainMonteCarlo(MCMC)chainsforthephylogeneticanal- frequencyrange(0–15%),neutralmutationstothemidfre-
ysestoconverge,weonlyused41samplescollectedbetween quencyrange(15–75%),andadaptivemutationstothehigh
December24,2019andJanuary31,2020(labeledas“January frequency range (75–100%). We then estimate the site fre-
sequences” in our analysis) and took 20 additional samples quency spectra by comparing the main alignments to an
ancestralsequence.ForSARS-CoV-2,theancestralsequence
duringthenextmonthtomatchwiththe120samplesused
is the earliest sample collected from Wuhan, Wuhan/
forpH1N1bytheendofthesecondmonthofsampling(see
IPBCAMS-WH-01/2019,andforpH1N1itistheearliestsam- D
supplementaryfig.S1,SupplementaryMaterialonline). o
ple from Mexico, ACQ99614jA/Mexico/4108/2009. The wn
PhylogeneticAnalyses choice of the cut-off frequencies is based on the diffusion loa
d
approximationofallelefrequencieswhereby,atequilibrium, e
We use BEAST v1.10 (Suchard et al. 2018) for the Bayesian d
phylogeneticanalysisoftheentiredatasetusinganHKYþC most deleterious and adaptive mutations are confined to fro
m
substitution model with a Laplace prior (mean¼0 and frequencies <15% and >75%, respectively (Bhatt et al. h
2011). It has further been shown that the exact choice of ttp
scale¼100)onthecoalescentgrowthrate,aLognormalprior s
(mean¼1 and SD¼2) on the coalescent population size, the cut-off frequency for the three frequency classes does ://ac
notsignificantlychangetheestimatednumberofnonneutral a
andacontinuoustimeMarkovchainpriorontheevolution- d
sites(Bhattetal.2011). em
aclroyccklomckodraetlea.nFodretxhpeofniresnttpiaalrtanodftlhoegiasntiaclygsrios,wwtheucsoealaesscterinctt We use the ORF1ab of SARS-CoV-2 and HA of pH1N1 ic.ou
influenzaforthisanalysis.Wealsocarriedoutasimilaranalysis p
demographicmodelsofSARS-CoV-2evolutionandonlythe .c
for the other genes of SARS-CoV-2, including the S gene. o
m
logistic growth model for the HA segment of pH1N1. To
However,becauseofthelimitedgeneticdiversitypresentin /m
quantify the relative fit of the two coalescent models for b
thesequencesandtheirrelativelyshortsize,ourgeneralized e
SARS-CoV-2,wecomputetheirlogmarginallikelihoodsusing /a
the generalized stepping-stone sampling method and com- MtimcDatoentahlde–nKurmeibtmeraonfnteosntnuesuintgraaldsiatpest.-aT-hraetreeaissounnafobrlecthoooess-- rticle
pare their Bayes factors (Fan et al. 2011; Baele and Lemey ingtoanalyzetheHAisprimarilybecauseoftheavailabilityof /39/2
2014; Baele et al. 2016). We note that although the growth alargenumberofsequencesfromGenBankforthefullcoding /m
s
coefficient of the two coalescent models is not expected to regionofHAwhichallowsustorandomlysample60align- ac
0
converge to the same value (as they correspond to intrinsi- mentspermonthfromApril2009toMarch2010. 0
9
callydifferentpopulationdynamics),thetwomodelsarefor- /6
5
Supplementary Material 0
mallynested,andthelikelihoodfunctionoftheexponential 9
5
2
modelexistsasalimitofthelikelihoodfunctionofthelogistic Supplementary data are available at Molecular Biology and 3
b
model.Therefore,eventhoughthegrowthrateiseffectivelya Evolutiononline. y g
nuisancevariableforourpurposes,thelikelihoodsofthetwo ue
Acknowledgments s
modelscanbecomparedviaalikelihoodratiotesttoselect t o
n
thebetter-fittingone(PybusandRambaut2002). This work was supported by the Biotechnology and 1
4
Foreachdataset,weperformMCMCrunsfor100million Biological Science Research Council (BB/M011224/1 to F
e
steps, sample trees every 10,000 steps, and remove the first M.G) and the European Research Council (101001623- bru
10% of the steps as burn-in. We ensure that the effective PALVIREVOLtoA.K.). ary
sample size for every parameter of interest is >200 using 20
Data Availability 2
Tracer v1.7 (Rambaut et al. 2018). For the second part of 2
theanalysis,weuseatwo-parametermolecularclockmodel Datasetsrequiredtoreproducetheanalysesareavailableat
withonestrictclockrateforterminalbranchesandonestrict https://github.com/mg878/twoclock_rate. The adapt-a-rate
clock for internal branches, using the same priors as before. package is also available at https://github.com/jnarag/
Thismolecularclockmodelisaversionofafixedlocalclock teaspoon.
whereby rather than having a single global rate on all
branches,theterminalandinternalbranchesareallowedto References
evolveaccordingtodifferentevolutionaryrateswhereasrate
BaeleG,LemeyP.2014.Bayesianmodelselectioninphylogeneticsand
constancy is assumed along the respective branches (Yoder genealogy-based population genetics. In: Bayesian phylogenetics,
and Yang 2000; Drummond and Suchard 2010). We use an methods, algorithms, and applications. Boca Raton (FL): CPC
Press.p.59–93.
MCMCchainoflength50millionsteps,samplingevery1,000
BaeleG,LemeyP,SuchardMA.2016.Genealogicalworkingdistributions
steps and evaluate sampling of the parameters of interest
forBayesianmodeltestingwithphylogeneticuncertainty.SystBiol.
usingTracerv1.7. 65(2):250–264.
7
MBE
Ghafarietal. . doi:10.1093/molbev/msac009
BhattS,HolmesEC,PybusOG.2011.Thegenomicrateofmolecular Kistler KE, Bedford T. 2021. Evidence for adaptive evolution in the
adaptation of the human influenza A virus. Mol Biol Evol. receptor-binding domain of seasonal coronaviruses OC43 and
28(9):2443–2451. 229e.eLife10:e64509.
BhattS,KatzourakisA,PybusOG.2010.Detectingnaturalselectionin KryazhimskiyS,PlotkinJB.2008.ThepopulationgeneticsofdN/dS.PLoS
RNA virus populations using sequence summary statistics. Infect Genet.4(12):e1000304.
GenetEvol.10(3):421–430. LamTTY,ZhuHC,GuanY,HolmesEC.2016.Genomicanalysisofthe
Bhatt S, Lam TT,Lycett SJ, Leigh Brown AJ, Bowden TA,Holmes EC, emergence,evolution,andspreadofhumanrespiratoryRNAviruses.
GuanY,WoodJLN,BrownIH,KellamP,etal.2013.Theevolutionary AnnuRevGenomicsHumGenet.17:193–218.
dynamicsofinfluenzaAvirusadaptationtomammalianhosts.Phil LuJ,duPlessisL,LiuZ,HillV,KangM,LinHF,SunJF,FrancoisS,Kraemer
TransRSocB.368(1614):20120382. MUG,FariaNR,etal.2020.GenomicepidemiologyofSARS-CoV-2in
ClarkAG,WhittamTS.1992.Sequencingerrorsandmolecularevolu- GuangdongProvince,China.Cell181(5):997–1003.e9.
tionaryanalysis.MolBiolEvol.9(4):744–752. LuJ,RaghwaniJ,PryceR,BowdenTA,ThezeJ,HuangSQ,SongYC,Zou
DeMaioN,WalkerC,BorgesR,WeilgunyL,SlodkowiczG,GoldmanN. LR,Liang LJ, Bai R, et al. 2018.Molecular evolution,diversity,and
2020. Issues with SARS-CoV-2 sequencing data.Available from: adaptationofinfluenzaA(H7N9)virusesinChina.EmergInfectDis.
https://virological.org/t/issues-with-sars-cov-2-sequencing-data/473. 24(10):1795–1805. D
DrummondAJ,SuchardMA.2010.Bayesianrandomlocalclocks,orone MeyerAG,SpielmanSJ,BedfordT,WilkeCO.2015.Timedependenceof ow
du rPaletesstisoLru,lMetchCermonaell.JTB,MZCarBeibols.k8i:1A1E4,. Hill V, Ruis C, Gutierrez B, ebvreoalukt.iVoniraursyEmvoelt.r1ic(1s)d:vuervin0g06t.he2009pandemicinfluenzavirusout- nloa
d
Raghwani J, Ashworth J, Colquhoun R, Connor TR, et al. 2021. MurrayGGR,WangF,HarrisonEM,PatersonGK,MatherAE,HarrisSR, ed
EstablishmentandlineagedynamicsoftheSARS-CoV-2epidemic HolmesMA,RambautA,WelchJJ.2016.Theeffectofgeneticstruc- fro
intheUK.Science371(6530):708–712. tureonmoleculardatingandtestsfortemporalsignal.MethodsEcol m
Duch^eneS,Duch^eneDA.2020.Estimatingevolutionaryratesandtime- PybuEsvoOl.G7(,1R):a8m0–b8a9u.t A.2002.GENIE:estimatingdemographic history http
scalesfromtime-stampeddata.In:Themolecularevolutionaryclock. s
Cham(Switzerland):Springer.p.157–174. frommolecularphylogenies.Bioinformatics18(10):1404–1405. ://a
Duch^eneS,Duch^eneD,HolmesEC,HoSYW.2015.Theperformanceof Pybus OG, Rambaut A, Belshaw R, Freckleton RP, Drummond AJ, cad
the date-randomization test in phylogenetic analyses of time- Holmes EC. 2007. Phylogenetic evidence for deleterious mutation em
structuredvirusdata.MolBiolEvol.32(7):1895–1906. loadinRNAvirusesanditscontributiontoviralevolution.MolBiol ic
Evol.24(3):845–852. .o
Duch^ene S, Featherstone L, Haritopoulou-Sinanidou M, Rambaut A, u
RaghwaniJ,BhattS,PybusOG.2016.Fasteradaptationinsmallerpop- p
Lemey P, Baele G. 2020. Temporal signal and the phylodynamic .c
ulations:counterintuitiveevolutionofHIVduringchildhoodinfec- o
thresholdofSARS-CoV-2.VirusEvol.6(2):veaa061. m
Duch^eneS,HolmesEC,HoSYW.2014.Analysesofevolutionarydynam- tion.PLoSComputBiol.12(1):e1004694. /m
RambautA,DrummondAJ,XieD,BaeleG,SuchardMA.2018.Posterior b
ics in viruses are hindered by a time-dependent bias in rate esti- e
summarizationinBayesianphylogeneticsusingTracer1.7.SystBiol. /a
mates.ProcRSocB.281(1786):20140732. 67(5):901–904. rtic
EdgarracRyCa.n2d00h4i.gMhUthSrCoLuEg:hmpuutlt.iNpluecsleeiqcuAecnidcesRaleigs.n3m2(e5n)t:1w79it2h–h1i7g9h7a.ccu- RambautA,HolmesE.2009.Theearlymolecularepidemiologyofthe le/3
swine-origin A/H1N1 human influenza pandemic. PLoS Curr. 9
FanY,WuR,ChenMH,KuoL,LewisPO.2011.Choosingamongpar- /2
1:RRN1003. /m
titionmodelsinBayesianphylogenetics.MolBiolEvol.28(1):523–532.
FraserC,DonnellyCA,CauchemezS,HanageWP,VanKerkhoveMD, RochEaJ.E2P0C06,.SmCoitmhpJMar,isHounrsstofLDdN,H/doSldaerneMtimTGe,dCeopoepnedreJnE,tSfmoritchloNseHly,Freei-l sac0
HollingsworthTD,GriffinJ,BaggaleyRF,JenkinsHE,LyonsEJ,etal. 0
latedbacterialgenomes.JTheorBiol.239(2):226–235. 9
2009.Pandemicpotential ofastrainofinfluenzaA(H1N1):early /6
SharpPM,BailesE,ChaudhuriRR,RodenburgCM,SantiagoMO,Hahn 5
findings.Science324(5934):1557–1561. 0
BH.2001.Theoriginsofacquiredimmunedeficiencysyndromevi- 9
GhafariM,SimmondsP,PybusOG,KatzourakisA.2021.Amechanistic ruses:whereandwhen?PhilTransRSocLondB.356(1410):867–876. 52
evolutionarymodelexplainsthetime-dependentpatternofsubsti- SmithGJD,VijaykrishnaD,BahlJ,LycettSJ,WorobeyM,PybusOG,Ma 3 b
tutionratesinviruses.CurrBiol.31(21):4689–96.e5. SK,CheungCL,RaghwaniJ,BhattS,etal.2009.Originsandevolu- y g
GireSK,GobaA,AndersenKG,SealfonRSG,ParkDJ,KannehL,JallohS, tionarygenomicsofthe2009swine-originH1N1influenzaAepi- ue
MomohM,FullahM,DudasG,etal.2014.Genomicsurveillance demic.Nature459(7250):1122–1125. st o
elucidatesEbolavirusoriginandtransmissionduringthe2014out- SuYCF,BahlJ,JosephU,ButtKM,PeckHA,KoayESC,OonLLE,BarrIG, n 1
HedbgreeaJ,kL.yScceietntcSeJ,3R4a5m(6b2a0u2t):1A3.6290–1133.7R2e.al-timecharacterizationofthe VfluijeanykzraishrenvaeaDls,StmheithtraGnJsDit.io2n01f5r.oPmhyhloodsytnaadmaipctsatoifonH1tNo1i/m20m09unine-- 4 Fe
b
molecular epidemiology of an influenza pandemic. Biol Lett. drivenselection.NatCommun.6(1): ru
9(5):20130331. SuchardMA,LemeyP,BaeleG,AyresDL,DrummondAJ,RambautA. ary
Ho SYW, Lanfear R, Bromham L, Phillips MJ, Soubrier J, Rodrigo AG, 2018.Bayesianphylogeneticandphylodynamicdataintegrationus- 2
0
CooperA.2011.Time-dependentratesofmolecularevolution.Mol ingBEAST1.10.VirusEvol.4:vey016. 22
Ecol.20(15):3087–3101. SullivanJ,JoyceP.2005.Modelselectioninphylogenetics.AnnuRevEcol
HoSYW,PhillipsMJ,CooperA,DrummondAJ.2005.Timedependency EvolSyst.36(1):445–466.
ofmolecularrateestimatesandsystematicoverestimationofrecent VolzE,BaguelinM,BhatiaS,BoonyasiriA,CoriA,Cucunub(cid:2)aZ,Cuomo-
divergencetimes.MolBiolEvol.22(7):1561–1568. Dannenburg G, Donnelly CA, Dorigatti I, FitzJohn R, et al. 2020.
HolmesEC.2003.Patternsofintra-andinterhostnonsynonymousvar- Report 5: phylogenetic analysis of SARS-CoV-2. MRC Centre for
iation reveal strong purifying selection in dengue virus. J Virol. GlobalInfectiousDiseaseAnalysis.ImperialCollegeLondon.
77(20):11296–11298. WertheimJO,PondSLK.2011.Purifyingselectioncanobscurethean-
HolmesEC,DudasG,RambautA,AndersenKG.2016.Theevolutionof cientageofvirallineages.MolBiolEvol.28(12):3355–3365.
Ebola virus: insights from the 2013–2016 epidemic. Nature WoodhamsM.2006.Candeleteriousmutationsexplainthetimedepen-
538(7624):193–200. dencyofmolecularrateestimates?MolBiolEvol.23(12):2271–2273.
HughesAL,HughesMAK.2007.Moreeffectivepurifyingselectionon YoderAD,YangZH.2000.Estimationofprimatespeciationdatesusing
RNAvirusesthaninDNAviruses.Gene404(1–2):117–125. localmolecularclocks.MolBiolEvol.17(7):1081–1090.
8
