Time scale evolution of avipoxviruses.

Avipoxviruses are divided into three clades: canarypox-like viruses, fowlpox-like viruses, and psittacinepox-like viruses. Several molecular clock and demographic models available in the BEAST package were compared on three avipoxvirus genes (P4b, cnpv186 and DNA polymerase genes), which enabled to determine that avipoxviruses evolved at a rate of 2-8×10(-5)substitution/site/year, in the range of poxviruses previously reported evolution rates. In addition, the date of mean time of divergence of avipoxviruses from a common ancestor was extrapolated to be about 10,000-30,000years ago, at the same period as modern poxvirus species. Our findings will facilitate epidemiological investigations on avipoxviruses' spread, origin and circulation.

Introduction

The Poxviridae virus family is divided into two subfamilies: the Chordopoxvirinae (containing vertebrate viruses) and the Entomopoxvirinae (containing insect viruses). The Chordopoxvirinae is further divided into ten genera: nine mammalian genera and the Avipoxvirus genus, subject of the present study (Haller et al., 2014, Virus Taxonomy, 2013). Based on basic genomic analyses (restriction fragment length polymorphism), avipoxviruses (APV) are currently divided into ten virus species (http://ictvonline.org/) but more recent phylogenetic studies identified three major clades: canarypox-like viruses, fowlpox-like viruses, and psittacinepox-like viruses (Jarmin et al., 2006).

Classically, APV were considered to be host species- or order- specific. Genus taxonomy had been based on this concept until recent studies showed that Otididae, Columbidae and Accipitridae can be infected by a large diversity of strains (Abdallah and Hassanin, 2013, Gyuranecz et al., 2013, Jarmin et al., 2006, Le Loc’h et al., 2014). The age of divergence between the three main APV clades has not been estimated. So far, only three genes have been widely studied: fpv094/cnpv121 (locus P4b, hereafter “P4b”), fpv140/cnpv186 (hereafter “cnpv186”) and fpv167/cnpv240 (hereafter “DNA polymerase”) (Gyuranecz et al., 2013, Jarmin et al., 2006, Le Loc’h et al., 2014). Besides, only four APV have been fully sequenced: one Fowlpox virus (Afonso et al., 2000), one Pigeonpox virus, one Penguinpox virus (Offerman et al., 2014), and one Canarypox virus (Tulman et al., 2004). Their analyses have shown large genomic rearrangements and suggest significant genomic diversity among APV.

Codivergence between DNA (and in some cases RNA) viruses and their hosts is considered as a main evolution path and should thus be associated with low mutation rates (Holmes, 2004, Holmes and Drummond, 2007, Villarreal and DeFilippis, 2000). Although this codivergence seems more likely in the case of papillomaviruses (millions of years of association, (Bernard et al., 2006)), it is less probable or at least more recent in the case of newly emerging viruses (e.g. the Severe Acute Respiratory Syndrome Coronavirus evolution (Peiris et al., 2004)). The literature suggests that a short sampling time interval (as compared to the real evolution timescale) biases the calculations (Ho et al., 2007).

The aim of the present preliminary study was to assess substitution rate and time to most recent common ancestor (TMRCA) of APV using the currently available data. Our results were compared with data available for other poxviruses and a hypothesis on the pathogen’s spread and circulation is given.

Materials and methods

Virus data sets

All available DNA sequences of APV fpv094/cnpv121 (locus P4b, hereafter “P4b”), fpv140/cnpv186 (hereafter “cnpv186”) and fpv167/cnpv240 (hereafter “DNA polymerase”) genes were retrieved from GenBank/EMBL databases. The sampling date (at least year, month and day when available) was collected for each sequence. In case no information was available, the 1st of July of the year of sequence publication was then considered as sample collection date. Only the oldest sequence was kept when identical sequences were available. Table 1 summarizes the details of included APV sequences. For each gene, DNA sequences were translated to amino acid (aa) sequences using BioEdit 7.2.5 (Hall, 1999), the aa sequences were aligned with ClustalX 2.1 (Larkin et al., 2007), and the original BioEdit format aa sequences were aligned with the ClustalX aa aligned sequences before they were reverse-translated in an aligned nucleotide file using BioEdit. The sequence of the APV pigeon/CAPM V-101/CZE/1965 P4b gene (accession number: GU982736) was removed due to many deletions resulting in non-functional short proteins; the P4b sequences of APV isolates from Hawaii EL1.1H1900, APAP16.2H03, and HAAMB.07H02 (accession numbers EF568402 to EF568404) were also removed because they were too short (116 bp). In order to exclude recombinant sequences (as done in Firth et al., 2010) which can bias the calculation, possible recombinant sequences were searched with the Recombinant Identification Program (RIP) (www.hiv.lanl.gov). None were identified. Overall, we analyzed 81 P4b gene sequences (partial open reading frame, ORF: 426 nucleotides) from specimens collected between 1865 and 2014; 31 cnpv186 genes (complete ORF: 933 to 1020 nucleotides) from specimens collected between 1965 and 2014; and 47 DNA polymerase genes (partial ORF: 555 nucleotides) from specimens collected between 1980 and 2014 (Table 1).

Table 1

Avipoxvirus sequences selected for the molecular clock analyses.

Gene	Virus clade⁎	Strain/sequence accession number	Sampling date	Estimated date
P4b	CNPV	Leptotila rufaxilla/PE/1865/HM627226	1865	1-Jul-1865
		Passer domesticus/ES/1911/HM627228	1911	1-Jul-1911
		Parus major/NO/1972/AY453173	20-Apr-1972	20-Apr-1972
		Sturnus vulgaris/US/1984/KC018063	1984	1-Jul-1984
		Serinus canaria/IT/1985/GQ180203	Dec 1985	15-Dec-1985
		Junco hyemalis hyemalis/US/1986/KC018022	1986	1-Jul-1986
		Quiscalus major/US/1989/KC018062	1989	1-Jul-1989
		Turdus philomelos/IT/1992/GQ180206	Feb 1992	15-Feb-1992
		Colinus virginianus ridgwayi/US/1993/KC018035	1993	1-Jul-1993
		Columba livia/NO/1996/AY453177	1996	1-Jul-1996
		Himatione sanguinea/US-HI/1996/KC018059	1996	1-Jul-1996
		Hemignathus virens/US-HI/1998/EF568393	1998	1-Jul-1998
		Hemignathus virens/US-HI/1998/EF568396	1998	1-Jul-1998
		Serinus canaria/DE/2001/AY530309	2001	1-Jul-2001
		Loxioides bailleui/US-HI/2002/EF568381	2002	1-Jul-2002
		Himatione sanguinea/US-HI/2002/EF568390	2002	1-Jul-2002
		Corvus brachyrhynchos/US/2003/DQ131891	2003	1-Jul-2003
		Corvus brachyrhynchos/US/2003/DQ131893	2003	1-Jul-2003
		Mimus polyglottos/US/2003/DQ131895	2003	1-Jul-2003
		Carpodacus mexicanus/US/2003/DQ131896	2003	1-Jul-2003
		Cardinalis cardinalis/US/2003/DQ131899	2003	1-Jul-2003
		Himatione sanguinea/US-HI/2003/EF568378	2003	1-Jul-2003
		Corvus hawaiiensis/US-HI/2003/EF568385	2003	1-Jul-2003
		Carduelis atrata/NL/2003/KC018038	2003	1-Jul-2003
		Ardea herodias/US/2004/DQ131898	2004	1-Jul-2004
		Hemignathus virens/US-HI/2004/EF568382	2004	1-Jul-2004
		Diomedea immutabilis/US-HI/2004/EF568383	2004	1-Jul-2004
		Serinus canaria/GB/2005/AM050375	<2006	1-Jul-2005
		Serinus canaria/GB/2005/AM050384	<2006	1-Jul-2005
		Columba livia/GB/2005/AM050386	<2006	1-Jul-2005
		Passer domesticus/GB/2005/AM050390	<2006	1-Jul-2005
		Hemignathus virens/US-HI/2005/EF568401	2005	1-Jul-2005
		Turdus migratorius/US/2005/KC018068	2005	1-Jul-2005
		Anthus berthelotii/PT/2006/EU883532	24-Sep-2006	24-Sep-2006
		Cyanistes caeruleus/ES-IC/2006/EU883533	10-Oct-2006	10-Oct-2006
		Corvus corone/IT/2007/GQ180211	Jun 2007	15-Jun-2007
		Parus major/HU/2007/EF634351	2007	1-Jul-2007
		Parus major/CZ/2008/EU798995	<2009	1-Jul-2008
		Serinus canaria/CL/2008/KC018060	2008	1-Jul-2008
		Chlamydotis undulata/MA/2009/LK021658	30-Jul-2009	30-Jul-2009
		Chlamydotis macqueenii/MA/2010/LK021653	6-Jul-2010	6-Jul-2010
		Serinus canaria/IR/2010/KC193679	10-Oct-2010	10-Oct-2010
		Erithacus rubecula/ES/2011/KF385499	20-Jul-2011	20-Jul-2011
		Chlamydotis macqueenii/MA/2011/LK021654	15-Oct-2011	15-Oct-2011
		Chlamydotis macqueenii/MA/2012/LK021655	1-Mar-2012	1-Mar-2012
		Chlamydotis macqueenii/UZ/2012/LK021657	10-Jun-2012	10-Jun-2012
		Spheniscus magellanicus/BR/2012/KC588955	1-Jul-2012	1-Jul-2012
		Spheniscus magellanicus/BR/2012/KC588956	1-Jul-2012	1-Jul-2012
		Spheniscus magellanicus/BR/2012/KC588958	1-Nov-2012	1-Nov-2012
		Spheniscus magellanicus/BR/2012/KC588959	1-Nov-2012	1-Nov-2012
		Spheniscus magellanicus/BR/2012/KC588960	1-Nov-2012	1-Nov-2012
		Burhinus oedicnemus/MA/2013/LN795883	5-Feb-2013	5-Feb-2013
		Spheniscus magellanicus/BR/2013/KF516679	15-Feb-2013	15-Feb-2013
		Chlamydotis macqueenii/AE/2014/LN795884	4-May-2014	4-May-2014
	FWPV	Fringilla coelebs/ES/1936/HM627225	1936	1-Jul-1936
		Perdrix perdrix/IT/1978/GQ180201	Dec 1978	15-Dec-1978
		Columba livia/US/1980/KC018001	1980	1-Jul-1980
		Phoebastria immutabilis/US/1983/KC017986	1983	1-Jul-1983
		Buteo jamaicensis/US/1985/KC018006	1985	1-Jul-1985
		Phalacrocorax pelagicus/US-AK/1989/KC017982	1989	1-Jul-1989
		Cygnus buccinator/US/1989/KC017995	1989	1-Jul-1989
		Uria aalge/US/1991/KC017985	1991	1-Jul-1991
		Gallus gallus domesticus/US-HI/1993/EF568397	1993	1-Jul-1993
		Gallus gallus domesticus/US-HI/1993/EF568398	1993	1-Jul-1993
Milvus milvus/ES/2000/KC018010		2000	1-Jul-2000
Falco peregrinus/DE/2001/KF956002		26-Feb-2001	26-Feb-2001
Falco/AE/2002/AY530306		2002	1-Jul-2002
Haliaeetus albicilla/DE/2002/KF956003		19-Aug-2002	19-Aug-2002
Otis tarda/HU/2003/KC017970		2003	1-Jul-2003
Otis tarda/ES/2003/KC017974		2003	1-Jul-2003
Pavo cristatus/HU/2003/KC017975		2003	1-Jul-2003
Haliaeetus albicilla/IT/2004/AB576861		2004	1-Jul-2004
Spilornis cheela/TW/2008/HQ441566		13-Jan-2008	13-Jan-2008
Spilopelia senegalensis/IN/2009/HM481408		Jan 2009	15-Jan-2009
Calypte anna/US/2009/JX418296		2009	1-Jul-2009
Gallus gallus domesticus/CN/2010/HM623675		<2011	1-Jul-2010
Spheniscus magellanicus/BR/2013/KC588962		1-Jan-2013	1-Jan-2013
Anas platyrhynchos domesticus/CN/2013/KJ192189		23-Dec-2013	23-Dec-2013
PSPV (?)	Coturnix japonica/IT/1963/GQ180200	1963	1-Jul-1963
PSPV	Amazona ochrocephala/US/1980/KC018069	1980	1-Jul-1980
	Agapornis/DE/2003/AY530311	<2004	1-Jul-2003
cnpv186	CNPV	Serinus canaria/GB/2005/AM071512	<2006	1-Jul-2005
		Chlamydotis undulata/MA/2009/LK021677	30-Jul-2009	30-Jul-2009
		Chlamydotis undulata/MA/2009/LN795880	13-Nov-2009	13-Nov-2009
		Chlamydotis undulata/MA/2009/LN795881	14-Nov-2009	14-Nov-2009
		Chlamydotis undulata/MA/2009/LK021680	17-Nov-2009	17-Nov-2009
		Chlamydotis undulata/MA/2010/LK021678	21-Jun-2010	21-Jun-2010
		Chlamydotis macqueenii/MA/2010/LK021671	28-Jun-2010	28-Jun-2010
		Chlamydotis undulata/MA/2010/LN795879	1-Jul-2010	1-Jul-2010
		Chlamydotis macqueenii/AE/2011/LK021667	14-May-2011	14-May-2011
		Chlamydotis undulata/MA/2011/LN795878	24-Jun-2011	24-Jun-2011
		Chlamydotis macqueenii/MA/2011/LK021673	15-Oct-2011	15-Oct-2011
		Chlamydotis macqueenii/MA/2012/LK021674	1-Mar-2012	1-Mar-2012
		Chlamydotis macqueenii/UZ/2012/LN795877	10-Jun-2012	10-Jun-2012
		Parus major/DE/2013/KF955997	<2014	1-Jul-2013
		Chlamydotis macqueenii/AE/2014/LN795875	4-May-2014	4-May-2014
		Chlamydotis macqueenii/KZ/2014/LN795876	27-Jul-2014	27-Jul-2014
	FWPV	1965/LN795882	1965	1-Jul-1965
		Falco peregrinus/DE/2001/KF955998	26-Feb-2001	26-Feb-2001
		Haliaeetus albicilla/DE/2002/KF955999	19-Aug-2002	19-Aug-2002
		Grus grus/DE/2004/KF955996	24-Mar-2004	24-Mar-2004
		Pigeon/2005/AM071389	<2006	1-Jul-2005
		Meleagris gallopavo/GB/2005/AM071390	<2006	1-Jul-2005
		Gallus gallus domesticus GB/2005/AM071393	<2006	1-Jul-2005
		Gallus gallus domesticus/2005/AM071394	<2006	1-Jul-2005
		Buteo buteo/IT/2005/EF133691	Sep 2005	15-Sep-2005
		Spilopelia senegalensis/IN/2009/HM481416	Jan 2009	15-Jan-2009
		Gallus gallus domesticus/CN/2010/HM623676	<2011	1-Jul-2010
		Gallus gallus domesticus/EG/2011/JX464822	Jun 2011	15-Jun-2011
		Columba livia/ZA/2011/KJ801920	2011	1-Jul-2011
		Chlamydotis macqueenii/AE/2013/LK021684	10-Feb-2013	10-Feb-2013
		Chlamydotis macqueenii/AE/2013/LK021683	31-Mar-2013	31-Mar-2013
DNApol	CNPV	Uria aalge/US/1980/KC017931	1980	1-Jul-1980
		Sturnus vulgaris/US/1984/KC017954	1984	1-Jul-1984
		Leucosticte tephrocotis/US/1985/KC017959	1985	1-Jul-1985
		Junco hyemalis hyemalis/US/1986/KC017936	1986	1-Jul-1986
		Junco hyemalis hyemalis/US/1986/KC017927	1986	1-Jul-1986
		Carpodacus mexicanus/US-HI/1987/KC017877	1987	1-Jul-1987
		Quiscalus major/US/1989/KC017922	1989	1-Jul-1989
		Grus canadensis/US/1992/KC017958	1992	1-Jul-1992
		Colinus virginianus ridgwayi/US/1993/KC017941	1993	1-Jul-1993
		Pica hudsonia/US/1997/KC017951	1997	1-Jul-1997
		Corvus brachyrhynchos/US/1999/KC017956	1999	1-Jul-1999
		Corvus brachyrhynchos/US/1999/KC017945	1999	1-Jul-1999
		Circus cyaneus/ES/2000/KC017897	2000	1-Jul-2000
		Carduelis atrata/NL/2003/KC017952	2003	1-Jul-2003
		Turdus migratorius/US/2005/KC017923	2005	1-Jul-2005
		Otis tarda/HU/2005/KC017864	2005	1-Jul-2005
		Parus major/HU/2007/KC017862	2007	1-Jul-2007
		Pyrrhula pyrrhula/BE/2008/KC017957	2008	1-Jul-2008
		Camarhynchus pallidus/EC/2008/KC017949	2008	1-Jul-2008
		Mimus parvulus/EC/2008/KC017947	2008	1-Jul-2008
		Geospiza fortis/EC/2008/KC017946	2008	1-Jul-2008
		Serinus canaria/CL/2008/KC017884	2008	1-Jul-2008
		Chlamydotis undulata/MA/2009/LN795874	30-Jul-2009	30-Jul-2009
	FWPV	Columba livia/US/1980/KC017929	1980	1-Jul-1980
		Haliaeetus leucocephalus/US-AK/1981/KC017900	1981	1-Jul-1981
		Phoebastria immutabilis/US/1983/KC017904	1983	1-Jul-1983
		Buteo jamaicensis/US/1985/KC017914	1985	1-Jul-1985
		1987/M31638	1987	1-Jul-1987
		Zenaida macroura/US/1987/KC017912	1987	1-Jul-1987
		Spheniscus demersus/SA/1988/KJ859677	1988	1-Jul-1988
		Cygnus buccinator/US/1989/KC017909	1989	1-Jul-1989
		Phalacrocorax pelagicus/US-AK/1989/KC017899	1989	1-Jul-1989
		Aix sponsa/US/1991/KC017910	1991	1-Jul-1991
		Uria aalge/US/1991/KC017902	1991	1-Jul-1991
		Branta canadensis/US/1992/KC017913	1992	1-Jul-1992
		Haliaeetus leucocephalus/US/1993/KC017932	1993	1-Jul-1993
		Columba livia/US-HI/1994/KC017868	1994	1-Jul-1994
		Milvus milvus/ES/2000/KC017896	2000	1-Jul-2000
		Buteo buteo/HU/2000/KC017861	2000	1-Jul-2000
		Otis tarda/ES/2003/KC017890	2003	1-Jul-2003
		Polytelis swainsonii/CL/2004/KC017883	2004	1-Jul-2004
		Falco peregrinus/HU/2005/KC017858	2005	1-Jul-2005
		Columba livia/HU/2005/KC017856	2005	1-Jul-2005
		Spheniscus magellanicus/AR/2007/KC017905	2007	1-Jul-2007
		Anas platyrhynchos domesticus/CN/2013/KM281933	23-Dec-2013	23-Dec-2013
		Domestic mallard duck/CN/2014/KM281932	2-Jan-2014	2-Jan-2014
	PSPV	Amazona ochrocephala/US/1980/KC017925	1980	1-Jul-1980

CNPV: canarypox-like viruses, FWPV: fowlpox-like viruses, PSPV: psittacinepox-like viruses, (?): non confirmed clade classification; <: before.

Evolutionary analyses

Phylogenies taking sampling time into account were estimated using the Bayesian Markov chain Monte Carlo (MCMC) inference methods available in the BEAST package (v.1.4.8, (Drummond and Rambaut, 2007)). The analyses presented were run using the Tamura and Nei 1993 model (TN93) with gamma-distributed rate: the closest substitution model available in BEAST to the best fitting nucleotide substitution model estimated by means of hierarchical likelihood ratio approach using Mega v 6.06 (Tamura et al., 2013); and the second best nucleotide substitution model as per MrModeltest (Nylander, 2004). We used empirical base frequencies and three partitions into codon positions. Different combinations of demographic models and clock models were compared as suggested (Drummond and Rambaut, 2007): constant size or extended Bayesian skyline plot models, strict or relaxed (lognormal or exponential) clocks. At least 200 million MCMC iterations were run for each gene and each demographic model/clock model combination. Subsampling was performed every 10,000 generations to decrease autocorrelation between model parameter samples. The estimation of parameters and divergence time was carried out using Tracer (http://beast.bio.ed.ac.uk/Tracer). A correct mixing of MCMC was verified by effective sampling size (ESS) calculations available in Tracer. The best fitting model was also assessed using “Model Comparison” available in Tracer and choosing the lowest Akaike’s information criterion (AICM) value. The dN/dS ratios (ω) were calculated using the single-likelihood ancestor (SLAC) method available in the HyPhy package available on the http://www.datamonkey.org/ server (Kosakovsky Pond, 2005, Kosakovsky Pond and Frost, 2005).

Results

Phylogenetic trees based on three APV genes (P4b, cnpv186, and DNA polymerase) DNA sequences were constructed and all confirmed the common classification in the genus: canarypox-like viruses, fowlpox-like viruses, and psittacinepox-like viruses (Fig. 1 ). To assess both substitution rate per site per year and TMRCA, different combinations of demographic models and clock models were compared: constant size or extended Bayesian skyline plot models, strict or relaxed (lognormal or exponential) clocks. The lowest AICM values were obtained with a relaxed exponential clock and the extended Bayesian skyline plot model for both cnpv186 and the DNA polymerase genes (AICM values of 11,673.74 ± 0.11 standard error (SE) and 7,316.26 ± 0.12 SE for cnpv186 and DNA polymerase genes, respectively, Table 2 ). For the P4b gene, the AICM value was slightly lower with a strict clock and the Extended Bayesian Skyline Plot model than with a relaxed exponential clock and the extended Bayesian skyline plot model (AICM values of 8,492.42 ± 0.21 SE and 8,493.49 ± 0.44 SE, respectively). We however decided to use the same combination for our three genes for consistency: a relaxed exponential clock and the extended Bayesian skyline plot model (Table 2). Analyses were also performed using the General Time Reversible substitution model with gamma-distributed rate and invariant sites (GTR + G + I, best fitted model as per MrModeltest), which gave similar results (data not shown).

Fig. 1

Phylogenetic trees analyzing all available avipoxvirus sequences of the (A) P4b, (B) cnpv186, and (C) DNA polymerase genes. Numbers at main tree nodes represent the age of the nodes in years with their [95% HPD]. “Accession number_estimated collection date” are indicated on the tree and can be linked to full strains names on Table 1. Sub-clades as described in (Le Loc’h et al., 2014) are indicated on (A).

Table 2

Mean and 95% HPD of the Bayesian posterior estimates of substitution rate and TMRCA for avipoxviruses.

Gene	Clock model1	Demographic model2	Substitution rate (subst/site/year)	Substitution rate, 95% HPD	TMRCA3 (ybp)	TMRCA3, 95% HPD	AICM4 ± SE
P4b	SC	CS	3.76E-05	2.99E-06–7.45E-05	10,133	2,849–34,025	8,516.62 ± 0.09
		EBSP	1.77E-05	1.59E-08–3.07E-05	20,246	4,754–105,117	8,492.42 ± 0.21
	LRC	CS	4.16E-05	1.85E-06–9.05E-05	11,087	2,086–48,467	8,540.54 ± 0.30
		EBSP	1.81E-05	3.82E-07–3.38E-05	23,382	3,551–110,916	8,561.01 ± 0.51
	ERC	CS	4.57E-05	2.27E-08–1.04E-04	11,315	1,761–59,744	8,502.61 ± 0.28
		EBSP	1.83E-05	2.96E-08–4.17E-05	29,258	1,861–197,714	8,493.49 ± 0.49
cnpv186	SC	CS	2.76E-05	5.99E-09–6.46E-05	20,107	3,705–163,817	11,706.96 ± 0.07
		EBSP	2.69E-05	1.51E-10–6.53E-05	22,008	2,824–155,499	11,698.32 ± 0.04
	LRC	CS	5.31E-05	2.44E-09–1.41E-04	11,446	783–96,743	11,680.38 ± 0.08
		EBSP	4.02E-05	3.32E-08–1.06E-04	15,247	841–110,873	11,676.37 ± 0.07
	ERC	CS	7.94E-05	4.46E-09–2.22E-04	7,792	595–75,721	11,677.74 ± 0.06
		EBSP	5.99E-05	4.85E-08–1.67E-04	10,891	963–93,707	11,673.74 ± 0.11
DNApol	SC	CS	4.90E-05	2.93E-08–1.06E-04	9,962	2,052–65,073	7,338.44 ± 0.10
		EBSP	5.34E-05	3.64E-09–1.07E-04	9,332	2,503–44,110	7,327.55 ± 0.08
	LRC	CS	7.12E-05	9.90E-08–1.69E-04	8,369	911–58,310	7,340.27 ± 0.14
		EBSP	6.04E-05	4.05E-08–1.35E-04	9,711	1,438–59,555	7,346.87 ± 0.13
	ERC	CS	8.34E-05	1.48E-08–2.00E-04	7,331	772–56,547	7,323.73 ± 0.11
		EBSP	7.59E-05	2.62E-09–1.85E-04	8,682	950–60,718	7,316.26 ± 0.12

SC: strict clock, LRC: lognormal relaxed clock, ERC: exponential relaxed clock.

CS: constant size, EBSP: extended Bayesian skyline plot.

TMRCA: time to most recent common ancestor (ybp: years before present).

AICM: Akaike’s information criterion, SE: standard error; in bold italic font: selected models.

Our molecular clock analyses showed that APV diverged from a common ancestor approximately 10–30 thousand years ago: 8,682 (95% Highest Posterior Density (HPD): 950–60,718), 10,891 (95% HPD: 963–93,707), and 29,259 (95% HPD: 1871–197,714) years ago for DNA polymerase, cnpv186, and P4b, respectively (Fig. 1). Interestingly, fowlpox-like and psittacinepox-like viruses shared a more recent common ancestor than canarypox-like viruses: they diverged 5–16 thousand years ago (5,034 (95% HPD: 625–34,636) years ago for DNA polymerase, 16,012 (95% HPD: 1208–117,230) years ago for P4b, no psittacinepox-like viruses sequence are available for cnpv186, Fig. 1).

When considering APV collected the same year in the same location such as Chlamydotis undulata/MA/2009/LK021680, C. undulata/MA/2009/LN795880, and C. undulata/MA/2009/LN795881, all three isolated from houbara bustards in Morocco in 2009, the cnpv186 phylogeny and molecular clock analysis (Fig. 1B) enabled to conclude to clear independent introductions of viruses and rule out the hypothesis of virus spread and evolution within the site. C. undulata/MA/2009/LN795881 interestingly clustered in sub-clade B2 (sequence identical to C. undulata/MA/2011/LK021679 described in (Le Loc’h et al., 2014)) whereas C. undulata/MA/2009/LN795880 (sequence identical to C. undulata/MA/2011/LK021670 described in (Le Loc’h et al., 2014)) and C. undulata/MA/2009/LK021680 clustered in a distinct sub-clade (B1). The two sub-clades showed a divergence 3,455 (95% HPD: 222–27,041) years ago which is not compatible with an on-site evolution of the virus.

Selection pressure analyses showed that P4b, cnpv186 and DNA polymerase genes were under purifying selection with ω (dN/dS) values of 0.067 (with a 95% confidence interval (CI) of 0.057–0.078), 0.200 (95% CI: 0.183–0.219), and 0.065 (95% CI: 0.055–0.076) for P4b, cnpv186 and DNA polymerase genes, respectively. No positively selected site was found in any of the three genes studied.

Discussion

The present study aimed at assessing the substitution rate of APV and at estimating their TMRCA. We compared several clock and demographic models available in the BEAST package and found that APV evolved at a rate of 2–8 × 10−5  substitution/site/year. The date of mean time of divergence of APV from a common ancestor was extrapolated to be about 10,000–30,000 years ago.

We observed a threefold difference between substitution rates and TMRCA of P4b on the one hand and cnpv186 and DNA polymerase on the other hand (Table 2). DNA polymerase encodes a polymerase protein, cnpv186 an immunodominant protein, and P4b a core protein (Tulman et al., 2004). The slower evolution rate of P4b is therefore somewhat surprising: while viable DNA polymerase and cnpv186 mutants may be more difficult to generate, one could expect more flexibility on a core protein. The selection pressure may however be lower on P4b than on its counterparts. Esposito et al. showed that the genes coding for host range and virulence proteins vary the most for Smallpox Virus compared to genes involved in virus replication, suggesting that these genes are most targeted by selection pressure (Esposito, 2006). But the large confidence intervals (95% HPD) associated with the values we found (within the range of what had been calculated in previous studies) does not show a statistically significant difference (Table 2). Therefore, it is not possible to draw any conclusion. In addition, our analysis has a bias as it was carried out on three genes only: a very small fragment of avipoxvirus genomes (260–306 kbp). Cnpv186, P4b and DNA polymerase genes were selected on the basis of available sequence data but might have a substitution rate different from the rest of the avipoxvirus genome. The molecular clock analyses of avipoxviruses will gain in precision when more sequence data (and more full genomes) become available. To run the analyses with as much data as possible, we extrapolated sampling years for some sequences. This also represents a bias. When excluding the sequence data for which no collection year was available, the substitution rates and TMRCA results for P4b and DNA polymerase genes were almost unchanged, but cnpv186 gene had a tenfold lower calculated substitution rate and a TMRCA of 2,274 years only.

McLysaght et al. (2003) tested poxvirus genes for selection acting on the evolution of genes by looking at synonymous and nonsynonymous mutations. They observed positive selection for genes involved in virulence or for gene candidates for host–pathogen interaction (McLysaght et al., 2003). Indeed, positive selection of viruses could also contribute to an increase of substitution rate, which can be observed either for viruses inducing a strong host immune response or for those which have been the object of vaccination or intervention campaigns (Firth et al., 2010). While the houbara bustards we focus on in our projects in Morocco, the United Arab Emirates, Kazakhstan and Uzbekistan are vaccinated annually (Le Loc’h et al., 2015), APV vaccination is not systematic worldwide, thus making selection pressure assessment difficult. Our selection pressure analyses demonstrated a purifying selection of the three APV genes we studied with ω values ranging between 0.065 and 0.200, suggesting that a non-synonymous mutation has only 7–20% as much chance as a synonymous mutation of being fixed in the population. This finding might correlate with relatively weak vaccination/intervention pressure and may be linked with the low host specificity of avipoxviruses. Vaccination is likely not the sole factor driving potential selection pressure and ecological and epidemiological factors should be further investigated for a proper understanding of the observed selection phenomenon.

Li et al. (2007) worked on Smallpox Virus evolution and linked sequence data and the first description of smallpox in humans (4th century AD in China; suspicions of smallpox as early as 1122 BC). Their dating was based on time-structured sequence data (207–231 years before present (ybp) using strict and relaxed clocks, respectively; these values were discarded for obvious discrepancy with historical and epidemiological data), or with calibration with historical records of smallpox infection (1400–6300 ybp) (Li et al., 2007). However, the latter method proved erroneous as epidemiological records may not correlate with real origin of a pathogen (Hughes et al., 2010, Shchelkunov, 2009). The absence of ancient epidemiological data on APV prevented us from including any historical calibration in our analyses.

Several research teams have in the past warned evolutionary biologists studying DNA viruses about the biases associated with the use of heterochronous phylogenetic modeling designed for faster evolving RNA pathogens (Firth et al., 2010). In addition, it was shown that molecular evolution is artificially accelerated on short timescales (Ho et al., 2007). Indeed our datasets contain very recent sequences only (from 1865 to now) compared to the TMCRA of APV. Hence our results should be considered with caution. But Babkin and Shchelkunov estimated that Poxviridae evolve at a rate of 0.9–1.2 × 10–6  substitutions/site/year. They calculated the ancestors of poxviruses, orthopoxviruses, and modern poxvirus species dated approximately 500,000, 300,000, and 14,000 years, respectively (Babkin and Shchelkunov, 2006). The evolution rate of APV estimated in the present study is therefore within the range of what was observed for poxviruses. The divergence between APV clades (canarypox-like viruses, fowlpox-like viruses, and psittacinepox-like viruses) is as old as the divergence among modern poxvirus species. Our results therefore suggest that APV may keep their “Genus” classification and not become a sub-family as proposed in the past (Jarmin et al., 2006), even if taxonomic classifications could not be based solely on the time of divergence. In addition, when Firth et al. tried to account for the specificity of double stranded DNA genomes in evolution calculations, their most robust virus model (with the different clock and demographic models tested) was Variola virus, suggesting that heterochronous phylogenetic modeling may be used for poxviruses evolution calculations (Firth et al., 2010). They however emphasized that the TMRCA calculated for Variola virus based on the high substitution rates observed were likely erroneous. Our TMRCA values should really be considered as first estimates and not overinterpretated.

In conclusion, the present study based on three genes showed that APV evolve at a rate of 2–8 × 10−5  substitution/site/year and that the date of mean time of divergence of APV from a common ancestor is about 10,000–30,000 years ago. This new finding should facilitate future epidemiological investigations on virus spread, origin and circulation.

Author contributions

G.L.L., M.F.D. and S.B. designed the study and drafted the manuscript. G.L.L. and M.F.D. performed and interpreted the analyses.

Conflicts of interest

The authors declare no conflict of interest.