FULL PAPER Porcine Circovirus 3 (PCV-3) www.advancedscience.com A Shift in Porcine Circovirus 3 (PCV-3) History Paradigm: Phylodynamic Analyses Reveal an Ancient Origin and Prolonged Undetected Circulation in the Worldwide Swine Population Giovanni Franzo,* Wanting He, Florencia Correa-Fiz, Gairu Li, Matteo Legnardi, Shuo Su, and Joaquim Segalés 1. Introduction The identification of a new circovirus (Porcine circovirus 3, PCV-3) has raised a remarkable concern because of some analogies with Porcine circovirus 2 The Circovirus genus claimed vet- (PCV-2). Preliminary results suggest an extremely recent PCV-3 emergence erinarian attention shortly after its and high mutation rate. Retrospective studies prove its circulation at least discovery, since its members were con- sidered responsible of relevant diseases since the early 1990s, revealing that PCV-3 could have been infecting pigs in birds.[1–3] However, it was only in for an even longer period. Therefore, a new evaluation, based on an updated the middle 1990s, when Porcine circo- collection of PCV-3 sequences spanning more than 20 years, is performed virus 2 (PCV-2) and its related clinically using a phylodynamic approach. The obtained results overrule the previous and economically relevant syndromes PCV-3 history concept, indicating an ancient origin. These evidences are (thereafter named Porcine circovirus dis- eases; PCVD) were recognized, that this associated with an evolutionary rate far lower (10−5–10−6 substitution/site/ genus became the focus of an intensive year) than the PCV-2 one. Accordingly, the action of selective pressures on research activity.[4,5] PCV-3 open reading frames (ORFs) seems to be remarkably lower compared Since then, a remarkable collection of to those acting on PCV-2, suggesting either a reduced PCV-3 plasticity or a data and knowledge has been gathered less efficient host-induced natural selection. A complex and not-directional on PCV-2 biology, pathogenesis, epide- viral flow network is evidenced through phylogeographic analysis, indicating miology, control, and evolution.[6,7] One of the most astonishing findings was a long lasting circulation rather than a recent emergence followed by the evidence that, in spite of the sudden spreading. Being recent emergence has been ruled out, efforts should be PCVD emergence, PCV-2 has shared a devoted to understand whether its recent discovery is simply due to improved long path with domestic and wild swine detection capabilities or to the breaking of a previous equilibrium. populations.[8,9] Not only this but other swine infection scenarios have pointed Dr. G. Franzo, M. Legnardi Dr. F. Correa-Fiz Department of Animal Medicine IRTA Production and Health (MAPS) Centre de Recerca en Sanitat Animal (CReSA, IRTA-UAB) University of Padua Campus de la Universitat Autònoma de Barcelona Viale, dell’Università 16, 35020 Legnaro (PD), Italy Bellaterra 08913, Spain E-mail: giovanni.franzo@unipd.it Prof. J. Segalés Dr. W. He, Dr. G. Li, Prof. S. Su UAB MOE International Joint Collaborative Research Laboratory Centre de Recerca en Sanitat Animal (CReSA, IRTA-UAB) for Animal Health & Food Safety Campus de la Universitat Autònoma de Barcelona Jiangsu Engineering Laboratory of Animal Immunology Bellaterra 08913, Spain Institute of Immunology Prof. J. Segalés College of Veterinary Medicine Departament de Sanitat i Anatomia Animals Nanjing Agricultural University Facultat de Veterinària Nanjing 210000, China Universitat Autònoma de Barcelona The ORCID identification number(s) for the author(s) of this article Bellaterra 08913, Spain can be found under https://doi.org/10.1002/advs.201901004. © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and repro- duction in any medium, provided the original work is properly cited. DOI: 10.1002/advs.201901004 Adv. Sci. 2019, 6, 1901004 1901004 (1 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com out the role of the modern farming system in the rising of new long, spanning the region between nucleotides 1347 and 1803 multifactorial diseases.[4,10] based on the reference genome KT869077.1) in the dataset3. The history seems to repeat itself for a new, recently discov- All databases comprised sequences originating from 14 coun- ered porcine circovirus: Porcine circovirus 3 (PCV-3).[11] This virus tries over a time period between 1996 and 2018, although with is featured by a circular ssDNA genome of ≈2000 bases con- different number per country (Datasets S1–S3, Supporting taining three open reading frames (ORFs) identified so far,[11] Information). although only ORF1 and ORF2 have been characterized. ORF1, located on the positive strand, apparently codes for a single repli- case protein of 296–297 aa.[11,12] ORF2 is located on the negative 2.2. PCV-3 Origin, Evolutionary Rate, and Population Dynamics DNA viral strand and encodes the caspsid protein.[12] Despite the common genomic organization, PCV-3 is distantly related to The analysis performed on the dataset1 provided a time to most other known circoviruses, although a certain relation with bat recent common ancestor (tMRCA) of 1372.23 years before pre- and avian circoviruses has been suggested based on phyloge- sent (ybp), although with a high uncertainness [95 high poste- netics, codon bias and genome composition analysis.[13,14] rior density [95HPD]: 666.30–3780.91] (Figure 1a and Figure S1, PCV-3 was first identified in the USA in 2015 using a Supporting Information). The lower tMRCA estimate was pro- metagenomics approach in tissues from animals displaying vided by dataset3 (156.9 ybp [95HPD: 75.72–318.28]) (Figure 1b porcine dermatitis and nephropathy syndrome (PDNS) and and Figure S2, Supporting Information). Dataset2 provided an reproductive disorders.[11] Thereafter, it has been identified all intermediate tMRCA estimation (≈300 ybp) when the results of over the world, including Asia,[15–18] Europe[19–22] and South the ten independent runs were averaged. However, a substan- America,[11,23] in presence of several clinical syndromes like tial overlap in the confidence intervals compared to dataset1 PDNS,[11] reproductive disorders,[24,25] respiratory disease,[26,27] was observed [95HPD: 138.79–1333.4]. The different runs, and myocarditis.[28] Moreover, PCV-3 genome has been detected performed on 10 randomly generated sequence datasets pro- through in situ hybridization and immunohistochemistry in vided consistent results (Figure 1c and Figure S3, Supporting different tissue lesions,[11,28,29] supporting the potential etiolog- Information). ical role of PCV-3. However, its identification with comparable Molecular clock estimates for the complete genome prevalence both in healthy pigs[30,31] and wild boar[32–34] ques- dataset reported 2.35 × 10−5 [95HPD: 8.79 × 10−6–4.71 × 10−5], tions the pathogenic role of this virus or at least suggests the 4.71 × 10−5 [95HPD: 1.75 × 10−5–9.38 × 10−5], and 5.13 × 10−5 need of other concomitant factors to trigger overt disease. [95HPD: 1.32 × 10−5–1.92 × 10−5] substitution/site/year for the Nevertheless, based on the “PCV-2 experience,” a pressing ORF1, ORF2, and intergenic regions, respectively (Figure 1a). concern has been directed toward the study of the evolution Accordingly, the dataset2 estimates showed an evolutionary and origin of PCV-3: is PCV-3 a newly originated viral species rate in the range 10−4–10−5 substitution/site/year, consistent or is it an ancient one that has been circulating for a long time, among different runs (Figure 1c). Finally, the evolutionary only recently emerging as a potential threat for swine industry? rate estimated using the dataset3 was 2.88 × 10−4 [95HPD: The first attempts to answer this question, based on a lim- 1.85 × 10−4–4.41 × 10−4] (Figure 1b). ited number of samples collected over a short time period, The relative genetic diversity (Ne × t) was featured by a broad pointed out a recent PCV-3 origin (approximately in the new 95HPD, independently from the considered dataset or run. millennium) coupled with a noteworthy evolutionary rate.[35] Nevertheless, a trend toward a rise in the viral population size However, retrospective studies performed in Sweden[20] in 1993 was observed, peaking approximately in the 1980s (Figure 2). and Spain[36] and China[15] in 1996 demonstrated that PCV-3 has been circulating during several decades in domestic pigs. Remarkably, PCV-3 has been detected in the oldest samples 2.3. Phylogeographic Analysis so far tested in these studies, showing a marked limit in our knowledge due to scarce data availability and suggesting that Reconstruction of the viral spread over time demonstrated a this virus could have infected pigs for even a longer period.[12] relevant uncertainness (i.e., posterior probability lower than Because of the relevance in terms of PCV-3 epidemiology 0.9) in the ancestral country estimation (Figure 3). Therefore, understanding and control strategies development, the pre- accounting for this uncertainness, relationship among strains sent study attempts to re-evaluate PCV-3 history, population were estimated in terms of well-supported (i.e., Bayesian Factor dynamics and spreading patterns based on a wider sequence (BF) > 10) contacts among countries. dataset, spanning the broader collection time window currently Except for a certain interdataset variability, the overall available. The herein reported results throw a new light on picture supports the presence of three interconnected nuclei PCV-3 history and its evolutionary pathways. of viral spread, corresponding to Asia, Europe, and America, with Asia acting as an intermediary between Europe and America (Figure 4 and Figures S4 and S5, Supporting Infor- 2. Results mation). Within those main areas, several local migration routes were proven significant, being China the most likely 2.1. Dataset responsible for viral spread to other Asian countries (e.g., South Korea, Thailand, Japan, and Russia). A more complex A total of 187 complete genome sequences were included and less directional network could be inferred for European in the dataset1, 208 in the dataset2 and 421 (427 nucleotide countries. Viral spread appears also to occur between USA Adv. Sci. 2019, 6, 1901004 1901004 (2 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Figure 1. a) PCV-3 tMRCA and evolutionary rate based on dataset1. Left figure: density plot of the MRCA posterior probability. Right figure: density plot of the mean evolutionary rate posterior probability. Evolutionary rates of different genomic regions have been color coded. The 95HPD intervals are reported for both figures. b) PCV-3 tMRCA and evolutionary rate based on dataset3. Left figure: density plot of the MRCA posterior probability. Right figure: density plot of the mean evolutionary rate posterior probability. The 95HPD intervals are reported for both figures. c) PCV-3 tMRCA and evolutionary rate based on dataset2. Upper figure: box plot (left) and density plot (right) of the MRCA posterior probability. Lower figure: box plot (left) and density plot (right) of the mean evolutionary rate posterior probability. Results have been estimated performing ten independent runs based on randomly sampled sequences. The 95HPD intervals are reported for both figures. and Mexico. However, the contact with the South American 2.4. Selective Pressure Analysis country included in the study (i.e., Brazil) was more likely mediated by viral strain exchange with China rather than Selective pressure analysis demonstrated a clear dominance of other American countries. sites under negative pressure both in the Rep and Cap proteins Adv. Sci. 2019, 6, 1901004 1901004 (3 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Figure 2. PCV-3 genotype population dynamics reconstructed based on dataset2. Upper figure: mean relative genetic diversity (N e x t) of the worldwide PCV-3 population overtime. The results of the ten independent runs have been color coded. Lower figure: median and upper and lower 95HPD values are reported for each run. (Figure 5). Limited evidences were detected of pervasive 3. Discussion diversifying selection in both proteins: positions 19, 44, 45, and 122 of the Rep protein were proven statistically significant with The identification of PCV-3, a new porcine circovirus resem- the Fast Unconstrained Bayesian AppRoximation (FUBAR) bling the significant pig pathogen PCV-2 from several per- method only, while position 5 of the Cap was detected by fixed spectives, has raised a great interest toward the epidemiology effects likelihood (FEL) and FUBAR, and sites 56 and 137 by and biology of this virus. One of the most relevant questions the latter method only. to solve is about its origin and evolution. The number of rec- On the other hand, episodic diversifying selection was ognized swine viruses has remarkably increased in the last detected in positions 5, 44, and 97 in the Rep, and in positions years[10,12] raising a dichotomous question to this phenomenon: 5, 56, and 102 in the Cap proteins. have these viruses been recently introduced in a new host A comparison of the selective pressures acting on the two species?, or have they circulated for a long period, undetected genes showed a statistically significant different selective regimen in the domestic pig population and emerged only recently as a (i.e., different selection strength and proportion of sites under major threat because of concomitant factors or improved detec- selection) acting on both genes (p = 0.003), being the overall dN/ tion and research technologies? While the second explanation dS and proportion of sites higher in the ORF2 alignment. appears more likely for several pathogens, including PCV-2,[4,9] Homology modeling allowed an approximate reconstruction their natural prior-to-emergence history remains largely of the protein tertiary structure. Two sites (i.e., 5 and 44) under unknown. episodic diversifying selection were located on the Rep protein Predictably, similar questions arose on the PCV-3 origin and surface, as well as positions 19 and 45 (detected by FUBAR) some studies have meritoriously attempted to investigate this (Figure 6a). The remaining sites under diversifying selection issue.[35,37] Preliminary results suggested a recent PCV-3 emer- were buried in the protein predicted structure. In the Cap gence, which was located approximately at the beginning of the protein, all the amino-acids under episodic diversifying selection new millennium. However, all the mentioned studies included were located on the capsid surface (Figure 6b), while sites only sequences of PCV-3 strains collected over a short time detected by FUBAR were located inside the protein structure. frame; i.e., after 2015. Therefore, a poor precision on tMRCA Adv. Sci. 2019, 6, 1901004 1901004 (4 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Figure 3. Ancestral location scatter plot. Scatter plot representing the posterior probability of each ancestral location (color coded) prediction over time. The results of ten independent BEAST run are reported. estimation could be forecasted, especially considering the low period, imposing a relatively fast evolutionary rate. Accord- resolution in the collection date annotation (i.e., collection year ingly, PCV-3 was reported to be the fasted evolving circovirus,[35] resolution). In fact, tMRCA underestimation can also severely which appears suspicious considering the limited genetic bias the substitution rate estimation. A recent tMRCA implies variability reported so far. The detection of PCV-3 in retrospec- the current genetic heterogeneity originated over a limited time tive samples collected during the 1990s in several countries from Europe and Asia confirmed its underestimated origin and claimed further analysis based on an updated genetic information availability.[15,20,36] Remarkably, although in Fu et al., estimation of the “PCV-3 only clade” still showed a recent origin, they were able to antici- pate PCV-3 emergence in the middle of the previous century by including the PorkNW2/USA/2009 strain in the coalescent analysis.[37] This strain is genetically similar to PCV-3 but shows a lower genome size and, although it could be reasonably classified as a defective PCV-3 or a replicative intermediate,[38] its classification remains controversial, mining the reliability of the results. Substantially comparable results were obtained by Saraiva et al.;[38] however, the PCV-3 strains included in this report were collected after 2015 and only Asia and America were represented. Based on these premises, the present study aimed to recon- struct the evolution history and population dynamics of PCV-3 based on a high quality and updated sequence dataset spanning a sampling time longer than 20 years (1996–2018). Despite our efforts, the currently available data are still limited and poten- Figure 4. Network reporting the well-supported migration routes (BF > 10). PCV-3 spreading path among different countries estimated tially biased by the different diagnostic and sequencing activi- using ten independent BEAST runs (color coded) based on dataset2. The ties performed in different countries. Therefore, particular arrows’ size is proportional to the BF value. care was dedicated to evaluate result consistency by analyzing Adv. Sci. 2019, 6, 1901004 1901004 (5 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Figure 5. Line graph reporting dN–dS values estimated for each codon position of ORF1 (upper figure) and ORF2 gene (lower figure) using different methods. Sites detected to be under statistically significant positive selection by MEME are reported as red circles. different genome regions and randomly generated, down-sam- the reconstruction of more ancient events. Moreover, the use pled, datasets with different features (i.e., alignment length, of the protein coding region allowed to implement a model number of sequences, presence of complete coding regions, that, although it cannot be considered an actual codon model etc.).[39] Independently from the considered dataset, PCV-3 (accounting for differential synonymous and nonsynonyms origin was always backdated before 1900 at least. However, with substitution rates), depicted the heterogeneous substitution the dataset3 being the only exception, a far more ancient origin rates among different alignment regions and codon positions was supported, in the order of several centuries (dataset2) or in a more effective way. Several studies have pointed out the millennia (dataset1). The large overlapping in 95HPD between underestimation of the origin of rapidly evolving viruses and the two datasets and analysis runs, further support the consist- the occurrence of a “time-dependent rate phenomenon,” where ency of the results. The trend toward an increase in tMRCA viral evolutionary rates appear to vary over time, continuously with increasing sequence length supports the usefulness of decreasing along with the timescale of rate measurement.[41–43] adding informative sites to improve parameter estimation Among the possible causes of these phenomena, substitution accuracy.[40] Particularly, the full genome dataset displays some saturation, poor modeling of natural selection and inability to features that could have additionally contributed to tMRCA deal with the vast majority of substitutions occurring multiple back-estimate. The inclusion and independent modeling of times at a limited subset of sites (i.e., high rate heterogeneity more conserved regions, like the ORF1, could have allowed among sites), have been advocated.[41–44] Therefore, because Figure 6. a) Tertiary structure of the Rep protein estimated through homology model. When located on the protein surfaces, sites under pervasive and episodic diversifying selection have been highlighted in orange and red, respectively. b) Tertiary structure of the Cap protein estimated through homology model. When located on the protein surfaces, sites under pervasive and episodic diversifying selection have been highlighted in red. Adv. Sci. 2019, 6, 1901004 1901004 (6 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com of the higher number of informative sites and more realistic Despite the long lasting PCV-3 circulation in the swine popu- model, the complete genome-based estimations can likely be lation, an increase in the viral relative genetic diversity (Ne x t; considered a more reliable tMRCA approximation. However, i.e., a proxy of population size dynamics) was observed in a potential underestimation of PCV-3 tMRCA can still not be the last decades using a skyline plot, similarly to what previ- excluded based on its high diversity compared to other known ously described for PCV-2.[9] Common epidemiological causes circoviruses, stressing the need for further improvements in could thus be hypothesized, like the alteration in the long our mathematical modeling capabilities. lasting equilibrium due to pig raising conditions in the con- Apart from these considerations, which are far beyond the test of expanding intensive farming.[4] The estimated increase scope of the study, the achieved results consistently demon- in viral population size, largely anticipate the identification of strated that PCV-3 origin should have occurred centuries ago. PCV-3 in animals showing clinical signs. Therefore, the answer This scenario is further supported by the worldwide distribu- to whether the detected rise mirrors an increased pathogenic tion pattern of the virus, featured by strain collected in dif- role of PCV-3 or is simply due to a wider viral circulation (of ferent countries widely interspersed in the phylogenetic tree no clinical significance) in a bigger and more connected animal (Figures S1–S3, Supporting Information), similarly to what has population world, remains elusive. Although some evidences previously been reported by other authors.[12,38,45,46] appear to support a certain association between PCV-3 and The root and ancestral node location posterior probability clinical disease, contradictory reports have been published up was often low (as indicated in Figure 3), revealing a largely to date and more extensive studies should be performed. expected uncertainty considering the large time frame between The ancient origin of PCV-3 implies also a lower evolutionary the estimated tMRCA and the oldest available sequences. Long rate compared to PCV-2 and previous PCV-3 estimation, both branches and the lack of historical data hinder the inference reporting a substitution rate in the order of magnitude of 10−3 of the spatial history of older viral lineages with confidence. substitution/site/year.[9,35,50] On the contrary, the present study Additionally, the long branch length is likely to conceal addi- results showed a far lower rate (i.e., ≈10−5 substitution/site/ tional spatial movements between multiple locations.[47] Conse- year), which is consistent with the limited genetic variability quently, migration patterns were evaluated in terms of ‘contact’ so far observed and with the high similarity between recent among countries, minimizing the risk of over-interpretation of sequences and those obtained from early-mid 1990s.[12] This their timing and directionality. Altogether, three major nuclei of evidence could suggest a lower intensity of diversifying selec- local transmission (i.e., North America, Europe and Asia), con- tive pressures shaping the evolution of this virus. nected by long distance transmission events were consistently Accordingly, compared with previous studies evaluating identified. Such uncontrolled viral circulation can be easily the forces acting on PCV-2, the number of sites under posi- explained by the recent PCV-3 identification and by its frequent tive selection in the Cap was remarkably lower in PCV-3.[9,51] detection in healthy animals. The strain distribution along the Although further confirmation will be needed, this scenario tree and the best fitting of a symmetric migration model over could be due to a lower plasticity of PCV-3 or to a less intense the asymmetric one, poses in favor of a long lasting viral circu- host-induced natural selection, which could be tentatively lation rather than a recent emergence followed by progressive considered as an evidence of a lower PCV-3 virulence and/ introduction in different countries. Nevertheless, based on the or prolonged virus-host co-evolution, leading to a decreased network structure and following a parsimony criterion, China immune response stimulation. Unfortunately, no reliable seems to have played a pivotal role in the spreading of PCV-3 experimental in vitro and/or in vivo model is currently available both within and between continents, which could seem sur- to investigate the PCV-3 immunopathogenesis and interaction prising being China a minor exporter of live swine. Despite our within the host, which are likely the most relevant determi- attempt to limit the effect of uneven sequence availability from nants of its evolution. Overall, the field of PCV-3 immunology different countries by down-sampling the original dataset and is still in its infancy and further confirmation must be provided creating randomly generated ones, a certain bias in spreading to shed light on this hypothesis plausibility. Nevertheless, the pattern inference due to the more intense sequencing activity phylodynamic approach implemented in the present study, in China cannot be totally excluded. However, fully compa- based on a global viral sampling spanning more than 20 years, rable long and short range spreading pattern has already been can provide a useful and consistent depiction of the overall pat- described for other livestock pathogens, like PCV-2 and Infec- terns and determinants of PCV-3 evolution, avoiding assump- tious bronchitis virus (IBV),[9,48,49] supporting the plausibility of a tions and constraints of viral biology induced by experimental comparable scenario for PCV-3. Therefore, the actual presence conditions. of some preferential livestock-virus “highways” can be sug- In fact, the ORF1 gene demonstrated a tendency toward a gested and further efforts should be dedicated to highlight the lower substitution rate compared to the ORF2, which was underlying causes, investigating the factors affecting viral dis- reflected by a statistically significant difference in diversifying persal and introducing effective control strategies, if necessary. selection acting on the two coding regions. These results, cou- While our hypothesis followed a parsimony criterion, other pled with the location on the protein surface of sites under patterns could represent a more accurate depiction of PCV-3 diversifying selection, suggest at least a limited action of the spreading. However, the current lack of accurate data in pig host immune response in shaping PCV-3 evolution. However, flow from most of the countries considered in the present study it must be stressed that the homology based estimation of Rep (especially in the time period when initial PCV-3 dispersal likely and Cap protein conformation should be evaluated with caution occurred) prevents additional investigations of an association because of the absence of closely related experimentally derived between swine or swine products trades and viral spread. tertiary structures. Adv. Sci. 2019, 6, 1901004 1901004 (7 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Overall, the present study provides an updated represen- partition was used for dataset2 and dataset3.The substitution model was tation of PCV-3 origin, population dynamics and evolution, selected based on the BIC scores calculated using Jmodeltest version pointing out a quite ancient viral origin and a low evolutionary 2.1.7. Molecular clock, population dynamics model, and discrete trait substitution model (i.e., symmetric vs asymmetric migration rate) were rate compared to other circoviruses of clinical relevance. These selected by evaluation of the BF (i.e., the ratio of the compared model results could contribute to the evaluation of the actual PCV-3 marginal likelihoods, estimated using a Path sampling and Stepping relevance for swine industry and, possibly, to the planning of stone approach) as suggested by Baele et al.[57] Relaxed lognormal effective control strategies. It must be emphasized that the molecular clock,[58] skyline population model,[59] and symmetric present work just scratched the surface of PCV-3 history and migration rate were selected.[56] biology and future and constant re-evaluation of the present The tree and the model parameters were estimated over a 100 million generation Markov Chain Monte Carlo (MCMC) chain, sampling them results will be mandatory to update and improve the knowledge every 10 000 generations. Run results were accepted only if the mixing and of this emergent virus behavior. convergence, visually inspected using Tracer,[60] were adequate and the estimated sample size (ESS) was higher than 200, after discharging the first 20% generations as burn-in. Parameter estimation were summarized as median and 95HPD. The maximum clade credibility tree was estimated 4. Experimental Section using the treeannotator tool of the BEAST 1.8.2[55] package. BF of well- Dataset Preparation: All currently available PCV-3 sequences were supported migration rates was calculated using SpreaD3.[61] downloaded from Genbank (accessed on 2018 November 29) and Selective Pressure Analysis: The presence of protein sites under annotated with collection year and country when available. Sequences diversifying selection was evaluated using several dN–dS based lacking of these data were removed from the dataset. Despite the fact methods. To this purpose, all complete ORF1 and ORF2 sequences that the obtained sequence collection was the broader currently available, were downloaded from Genbank and aligned at codon level using the molecular epidemiology information was still limited and sparse. TranslatorX.[53] The achieved alignments were tested for pervasive To deal with unbalanced sequence availability, which could bias the diversifying selection using single-likelihood ancestor counting (SLAC), results, different sequence datasets were prepared to alternatively benefit FEL, and FUBAR[62,63] methods implemented in HyPhy.[64] Significance from the higher sequence length or number (i.e., representativeness): level was set to p < 0.05 for SLAC and FEL and to a posterior probability dataset1 included all available complete genome; dataset2 comprised higher than 0.9 for FUBAR. Sites were considered under pervasive long PCV-3 sequences (1000bp), including the Spanish ones diversifying selection when detected by at least two methods. obtained during a retrospective study conducted from the mid-1990s Since adaptive evolution often occurs in episodic bursts,[65] i.e., onward;[36]dataset3 was based on a region where the higher number of affecting a small subset of branches at individual sites, episodic PCV-3 sequence had a full coverage. diversifying selection was also investigated using MEME.[66] The MEME Each dataset was designed using the following approach: significance level was set p < 0.05. Dataset1: Complete genome sequences were divided in ORF1, ORF2, The tertiary structure of PCV-3 Rep and Cap proteins were estimated and intergenic regions (both intergenic regions were merged in a single by homology modeling using Phyre 2.0[67] to predict sites location under partition). Coding regions were aligned at amino-acid level and then diversifying selection. back translated to nucleotide sequence using the MAFFT algorithm[52] The action of selective pressures was also compared among different implemented in TranslatorX.[53] Poorly aligned sequences as well as those genes (ORF1 and ORF2) using the dNdSDistributionComparison. showing premature stop codons or frame-shift mutation were excluded bf implemented in HyPhy.[64] Particularly, the presence of a difference from the analysis. Recombination analysis was performed using RDP4[54] between the two ORFs in selective pressure strength, proportion of sites on complete genome alignment and each ORF independently. under selective pressure and selective regime (i.e. a combination of the Dataset2: This dataset was primarily designed to benefit from all the two factors) was evaluated. 1990s samples, including the partial ones describe by Klaumann et al.[36] Additionally, dataset2 was used to evaluate the impact of sampling bias in analysis results. To this purpose, ten independent datasets were generated by randomly sampling a maximum of ten sequences for Supporting Information each country-year pair, as described by Franzo et al.[48,49] The obtained Supporting Information is available from the Wiley Online Library or sequences were aligned using MAFFT version 7.271[52] and scanned for from the author. recombination events using RDP4.[54] Dataset3: All the available, partial and complete, PCV-3 sequences were aligned using MAFFT[52] and the alignment region with the highest sequence coverage was selected for further analysis. All sequences spanning the selected region were extracted, realigned, and scanned for Acknowledgements recombination events using RDP4.[54] G.F., W.H., S.S., and J.S. contributed equally to this work. W.T.H., G.R.L., Population Parameters Estimation/Phylodynamic Analysis: On each of and S.S. were financially supported by the National Key Research and the above-mentioned datasets, a tip dated serial coalescent analysis was Development Program of China (2017YFD0500101), the Natural Science performed using the Bayesian approach implemented in BEAST 1.8.2[55] Foundation of Jiangsu Province (BK20170721), and the China Association to primarily estimate tMRCA, evolutionary rate, and population dynamics for Science and Technology Youth Talent Lift Project (2017–2019). over time. Additionally, a discrete state phylogeography was performed J.S. and F.C.-F. thank the financial support of the E-RTA2017-00007-00-00 as described by Lemey et al.[56] Additionally, the implementation of the project from the Instituto Nacional de Investigación y Tecnologia Bayesian stochastic search variable selection (BSSVS) allowed a BF Agraria y Alimentaria (Spanish Government). G.F. thanks the financial test that identified the most parsimonious description of the spreading support of the Department of Animal Medicine, Production and Health, process. A BF > 10 was considered as suggestive of a significant University of Padua, (grant BIRD187958/18). migration pattern between country pairs. For the dataset1, the three partitions (ORF1, ORF2, and intergenic regions) were allowed independent substitution and clock models, while a single tree topology for the three regions was constrained. Additionally, Conflict of Interest ORF1 and ORF2 regions were further partitioned allowing independent evolution models for each codon position. On the other hand, a single The authors declare no conflict of interest. Adv. Sci. 2019, 6, 1901004 1901004 (8 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com Keywords [27] H. Shen, X. Liu, P. Zhang, L. Wang, Y. Liu, L. Zhang, P. Liang, C. Song, Transboundary Emerging Dis. 2018, 65, 264. epidemiology, evolution, porcine circovirus 3 (PCV-3), phylodynamic, [28] T. G. Phan, F. Giannitti, S. Rossow, D. Marthaler, T. Knutson, L. Li, phylogeography X. Deng, T. Resende, F. Vannucci, E. Delwart, Virol. J. 2016, 13, 184. [29] B. Arruda, P. Piñeyro, R. Derscheid, B. Hause, E. Byers, K. Dion, Received: April 29, 2019 D. Long, C. Sievers, J. Tangen, T. Williams, K. Schwartz, Emerging Revised: July 12, 2019 Microbes Infect. 2019, 8, 684. Published online: September 30, 2019 [30] S. Zheng, X. Wu, L. Zhang, C. Xin, Y. Liu, J. Shi, Z. Peng, S. Xu, F. Fu, J. Yu, W. Sun, S. Xu, J. Li, J. Wang, Transboundary Emerging Dis. 2017, 64, 1337. [31] G. Franzo, M. Legnardi, C. M. Tucciarone, M. Drigo, F. Klaumann, [1] I. Tischer, R. Rasch, G. Tochtermann, Zentralbl. Bakteriol., M. Sohrmann, J. Segalés, Vet. Rec. 2018, 182, 83. Parasitenkd., Infektionskrankh. Hyg., Abt. 1: Orig., Reihe A 1974, 226, [32] G. Franzo, C. M. Tucciarone, M. Drigo, M. Cecchinato, M. Martini, 153. A. Mondin, M. L. Menandro, Transboundary Emerging Dis. 2018, 65, [2] I. Tischer, H. Gelderblom, W. Vettermann, M. A. Koch, Nature 1982, 957. 295, 64. [33] F. Klaumann, A. Dias-Alves, O. Cabezón, G. Mentaberre, [3] D. Todd, Vet. Microbiol. 2004, 98, 169. R. Castillo-Contreras, M. López-Béjar, E. Casas-Díaz, M. Sibila, [4] J. Segalés, T. Kekarainen, M. Cortey, Vet. Microbiol. 2013, 165, 13. F. Correa-Fiz, J. Segalés, Transboundary Emerging Dis. 2019, 66, 91. [5] X. J. Meng, Virus Res. 2012, 164, 1. [34] G. Franzo, L. Grassi, C. M. Tucciarone, M. Drigo, M. Martini, [6] F. Madec, N. Rose, B. Grasland, R. Cariolet, A. Jestin, D. Pasotto, A. Mondin, M. L. Menandro, Transboundary Emerging Transanboundary Emerging Dis. 2008, 55, 273. Dis. 2019, 66, 1548. [7] M. A. Ssemadaali, M. Ilha, S. Ramamoorthy, Res. Vet. Sci. 2015, 103, [35] G. Li, W. He, H. Zhu, Y. Bi, R. Wang, G. Xing, C. Zhang, J. Zhou, 179. K.-Y. Yuen, G. F. Gao, S. Su, Adv. Sci. 2018, 5, 1800275. [8] G. Franzo, M. Cortey, A. M. M. G. de Castro, U. Piovezan, [36] F. Klaumann, G. Franzo, M. SohrmannF. Correa-Fiz, M. Drigo, M. P. J. Szabo, M. Drigo, J. Segalés, L. J. Richtzenhain, Vet. J. I. Núñez, M. Sibila, J. Segalés, Transboundary Emerging Dis. 2018, Microbiol. 2015, 178, 158. 65, 1290. [9] G. Franzo, M. Cortey, J. Segalés, J. Hughes, M. Drigo, Mol. [37] X. Fu, B. Fang, J. Ma, Y. Liu, D. Bu, P. Zhou, H. Wang, K. Jia, Phylogenet. Evol. 2016, 100, 269. G. Zhang, Transboundary Emerging Dis. 2018, 65, e296. [10] G. Fournié, L. Kearsley-Fleet, J. Otte, D. U. Pfeiffer, Vet. Res. 2015, [38] G. L. Saraiva, P. M. P. Vidigal, J. L. R. Fietto, G. C. Bressan, A. Silva 46, 114. Júnior, M. R. de Almeida, Virus Genes 2018, 54, 376. [11] R. Palinski, P. Piñeyro, P. Shang, F. Yuan, R. Guo, Y. Fang, E. Byers, [39] M. D. Hall, M. E. J. Woolhouse, A. Rambaut, Virus Evol. 2016, 2, B. M. Hause, J. Virol. 2017, 91, e01879. vew003. [12] F. Klaumann, F. Correa-Fiz, G. Franzo, M. Sibila, J. I. Núñez, [40] J. Felsenstein, Mol. Biol. Evol. 2006, 23, 691. J. Segalés, Front. Vet. Sci. 2018, 5, 315. [41] J. O. Wertheim, D. K. W. Chu, J. S. M. Peiris, S. L. Kosakovsky Pond, [13] G. Franzo, J. Segales, C. M. Tucciarone, M. Cecchinato, M. Drigo, L. L. M. Poon, J. Virol. 2013, 87, 7039. PLoS One 2018, 13, e0199950. [42] E. C. Holmes, J. Virol. 2003, 77, 3893. [14] G. Li, H. Wang, S. Wang, G. Xing, C. Zhang, W. Zhang, J. Liu, [43] P. Aiewsakun, A. Katzourakis, J. Virol. 2016, 90, 7184. J. Zhang, S. Su, J. Zhou, Virulence 2018, 9, 1301. [44] S. Y. W. Ho, R. Lanfear, L. Bromham, M. J. Phillips, J. Soubrier, [15] J. Sun, L. Wei, Z. Lu, S. Mi, F. Bao, H. Guo, C. Tu, Y. Zhu, W. Gong, A. G. Rodrigo, A. Cooper, Mol. Ecol. 2011, 20, 3087. Transboundary Emerging Dis. 2018, 65, 607. [45] G. Franzo, M. Legnardi, C. K. Hjulsager, F. Klaumann, L. E. Larsen, [16] X. Ku, F. Chen, P. Li, Y. Wang, X. Yu, S. Fan, P. Qian, M. Wu, Q. He, J. Segales, M. Drigo, Transboundary Emerging Dis. 2018, 65, 602. Transboundary Emerging Dis. 2017, 64, 703. [46] R. Fux, C. Söckler, E. K. Link, C. Renken, R. Krejci, G. Sutter, [17] S. Qi, M. Su, D. Guo, C. Li, S. Wei, L. Feng, D. Sun, Transboundary M. Ritzmann, M. Eddicks, Virol. J. 2018, 15, 25. Emerging Dis. 2019, 66, 1004. [47] M. I. Nelson, C. Viboud, A. L. Vincent, M. R. Culhane, S. E. Detmer, [18] W. Wang, W. Sun, L. Cao, M. Zheng, Y. Zhu, W. Li, C. Liu, X. Zhuang, D. E. Wentworth, A. Rambaut, M. A. Suchard, E. C. Holmes, J. Xing, H. Lu, T. Luo, N. Jin, BMC Vet. Res. 2019, 15, 60. P. Lemey, Nat. Commun. 2015, 6, 6696. [19] G. Franzo, M. Legnardi, C. K. Hjulsager, F. Klaumann, L. E. Larsen, [48] G. Franzo, P. Massi, T. M. Tucciarone, I. Barbieri, G. Tosi, J. Segales, M. Drigo, Transboundary Emerging Dis. 2018, 65, 602. L. Fiorentini, M. Ciccozzi, A. Lavazza, M. Cecchinato, A. Moreno, [20] X. Ye, M. Berg, C. Fossum, P. Wallgren, A. L. Blomström, Virus PLoS One 2017, 12, e0184401. Genes 2018, 54, 466. [49] G. Franzo, M. Cecchinato, G. Tosi, L. Fiorentini, F. Faccin, [21] T. Stadejek, A. Woźniak, D. Miłek, K. Biernacka, Transboundary C. M. Tucciarone, T. Trogu, I. Barbieri, P. Massi, A. Moreno, PLoS Emerging Dis. 2017, 64, 1350. One 2018, 13, e0203513. [22] P. J. Collins, J. McKillen, G. Allan, Vet. Rec. 2017, 181, 599. [50] C. Firth, M. A. Charleston, S. Duffy, B. Shapiro, E. C. Holmes, [23] C. Tochetto, D. A. Lima, A. P. M. Varela, M. R. Loiko, W. P. Paim, J. Virol. 2009, 83, 12813. C. M. Scheffer, J. I. Herpich, C. Cerva, C. Schmitd, S. P. Cibulski, [51] G. Franzo, C. M. Tucciarone, M. Cecchinato, M. Drigo, Sci. Rep. A. C. Santos, F. Q. Mayer, P. M. Roehe, Transboundary Emerging Dis. 2016, 6, 39458. 2018, 65, 5. [52] K. Katoh, D. Standley, Mol. Biol. Evol. 2013, 30, 772. [24] Y. Zou, N. Zhang, J. Zhang, S. Zhang, Y. Jiang, D. Wang, Q. Tan, [53] F. Abascal, R. Zardoya, M. J. Telford, Nucleic Acids Res. 2010, Y. Yang, N. Wang, Arch. Virol. 2018, 163, 2841. 38, W7. [25] S. Faccini, I. Barbieri, A. Gilioli, G. Sala, L. R. Gibelli, A. Moreno, [54] D. P. Martin, B. Murrell, M. Golden, A. Khoosal, B. Muhire, Virus C. Sacchi, C. Rosignoli, G. Franzini, A. Nigrelli, Transboundary Evol. 2015, 1, vev003. Emerging Dis. 2017, 64, 1661. [55] A. J. Drummond, M. A. Suchard, D. Xie, A. Rambaut, Mol. Biol. [26] S. L. Zhai, X. Zhou, H. Zhang, B. M. Hause, T. Lin, R. Liu, Evol. 2012, 29, 1969. Q.-L. Chen, W.-K. Wei, D. H. Lv, X. H. Wen, F. Li, D. Wang, Virol. J. [56] P. Lemey, A. Rambaut, A. J. Drummond, M. A. Suchard, PLoS 2017, 14, 222. Comput. Biol. 2009, 5, e1000520. Adv. Sci. 2019, 6, 1901004 1901004 (9 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.advancedsciencenews.com www.advancedscience.com [57] G. Baele, P. Lemey, T. Bedford, A. Rambaut, M. A. Suchard, [63] B. Murrell, S. Moola, A. Mabona, T. Weighill, D. Sheward, A. V. Alekseyenko, Mol. Biol. Evol. 2012, 29, 2157. S. L. Kosakovsky Pond, K. Scheffler, Mol. Biol. Evol. 2013, 30, [58] A. J. Drummond, S. Y. W. W. Ho, M. J. Phillips, A. Rambaut, PLoS 1196. Biol. 2006, 4, e88. [64] S. L. K. Pond, S. D. W. Frost, S. V. Muse, Bioinformatics 2005, 21, [59] A. J. Drummond, A. Rambaut, B. Shapiro, O. G. Pybus, Mol. Biol. 676. Evol. 2005, 22, 1185. [65] S. L. K. Pond, B. Murrell, M. Fourment, S. D. W. Frost, W. Delport, [60] A. Rambaut, A. J. Drummond, D. Xie, G. Baele, M. A. Suchard, Syst. K. Scheffler, Mol. Biol. Evol. 2011, 28, 3033. Biol. 2018, 67, 901. [66] B. Murrell, J. O. Wertheim, S. Moola, T. Weighill, K. Scheffler, [61] F. Bielejec, G. Baele, B. Vrancken, M. A. Suchard, A. Rambaut, S. L. K. Pond, PLoS Genet. 2012, 8, e1002764. P. Lemey, Mol. Biol. Evol. 2016, 33, 2167. [67] L. A. Kelley, S. Mezulis, C. M. Yates, M. N. Wass, M. J. E. Sternberg, [62] S. L. K. Pond, S. D. W. Frost, Mol. Biol. Evol. 2005, 22, 1208. Nat. Protoc. 2015, 10, 845. Adv. Sci. 2019, 6, 1901004 1901004 (10 of 10) © 2019 The Authors. Published by WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim