Total RNA was extracted from alfalfa leaves systemically infected by AMV. The extracted RNA was used as a template for RT-PCR amplification in a 20-μL reaction containing total RNA (50–5000 ng), 1 μL gDNA Remover, 1 μL primer AMV-R ( Supplementary Table S4 ), 10 μL of 2× Reaction Mix, 1 μL Enzyme Mix with 5 mM dNTP (each), and RNase-free water to 20 µL. The reaction was incubated at 25°C for 10 min, 42°C for 45 min, and 85°C for 5 s, then the PCR commenced with the system of 2× EasyTaq^® PCR Super Mix (+ dye) 25 µL, cDNA 2 µL, pimer cp-F 1 µL, primer cp-R 1 µL, and 21 µL of ddH₂O to total 50 µL (TransGen, China). PCR cycling conditions comprised an initial denaturation step at 94°C for 3 min, 35 cycles of denaturation at 94°C for 1 min, annealing at 57°C for 1 min, and extension at 72°C for 3 min, followed by a final extension at 72°C for 5 min. The expected RT-PCR products of the cp gene were 666 bp, using primer pairs cp F/cp R ( Supplementary Table S5 ), and together covered the entire length of the viral genome. The PCR product segments were electrophoresed in 1.0% agarose gels and purified by the BioTeq PCR quick Gel Extraction Kit (BioTeq, USA). The purified fragments were cloned into the pMD18-T vector (Takara, Dalian, China) and used to transform Escherichia coli JM110. Insert sequences were determined for at least three clones for each fragment using either the ABI (ABI BigDye 3.1, Applied Biosystems) or Beckman (GeXP with Genome Lab DTCS sequencing kit) system. Sequence data were assembled using DNASIS version 3.5 (Hitachi, Tokyo, Japan), Laser gene (DNASTAR, Madison WI, USA), or BIOEDIT version 5.0.9 (Hall, 1999).

The published cp gene sequences of AMV isolates from Iran, Spain, and other countries were collected from the GenBank database ( Supplementary Table S1 ).

All of the cp genes of AMV populations from China, Iran, and Spain, respectively, were aligned with MUSCLE (Edgar, 2004) as implemented in MEGA version 5.0 (Tamura et al., 2011). The best model of nucleotide substitution was determined by MODELTEST version 3.7 (Posada and Crandall, 2001). SplitsTree version 4 was used for the split-decomposition network analysis (Huson and Bryant, 2006). Maximum credibility clade (MCC) phylogenetic reconstructions were conducted using BEAST version 1.5.4 (Drummond and Rambaut, 2007). Markov Chain Monte Carlo (MCMC) simulations were run for 10⁷ generations to ensure convergence of all parameters. Branches with a posterior support probability of 0.50 were collapsed.

Population genetics analyses and diversity measures for all populations (Tajima’s D, Ewans-Watterson, Chakraborty, and Fu’s F) were calculated by ARLEQUIN (http://cmpg.unibe.ch/software/ arlequin3) (Excoffier et al., 2005). The best alignment of the sequences of every population was inputted into the ARLEQUIN version 3.0 and the corresponding analysis was performed.

Selective pressures operating at each codon were evaluated based on the difference between synonymous (d_S ) and nonsynonymous (d_N ) substitution rates for the cp gene calculated by MEGA X (Tamura et al., 2011). As a first approach, values of d_N −d_S <0, = 0, or >0 indicate purifying selection, neutral evolution, and positive selection, respectively.

Data were divided into different groups among bioclimatic regions ( Supplementary Figure S1 ). To determine the extent of geographic structure in AMV populations, BaTS version 1.0b2 (Parker et al., 2008) was used to compute the parsimony score (PS) (Slatkin and Maddison, 1989), the association index (AI) (Wang et al., 2001), and the maximum monophyletic clade size (MC) (Parker et al., 2008), and to assess the significance of these three statistics. The first 10% of sampled trees were discarded as burn-in and 104 randomizations were performed to estimate the null distributions of the three statistics.

The data sets from the three countries were analyzed respectively for the recombination using the methods of RDP, BOOTSCAN, CHIMERA, GENECONV, LARD, MAXCHI, SISCAN and 3SEQ implemented in RDP4 in the default configuration (Martin et al., 2010) as well as the GARD and PHI test in SplitsTree (Huson and Bryant, 2006). Only those recombination events predicted by at least five of these methods and with a p value < 0.05 were regarded as valid.

The parsimony networks revealed the existence of two major clades, one more divergent (Clade III) and one less so (Clades I and II) ( Supplementary Figure S2 ), which was consistent with the clustering identified by the MCC tree obtained with BEAST ( Supplementary Figure S3 ). In the MCC tree, Clade A (yellow) contained isolates from three climatic regions (plateau mountain climate, temperate continental climate, temperate monsoon climate), while Clade B (blue) comprised a single isolate from Inner Mongolia, and Clade C (purple) contained isolates from all of the climatic regions. Clade I in Supplementary Figure S2 corresponded with Clade A in Supplementary Figure S3 , and Clades II and III were consistent with Clades B and C, respectively. The software package RDP 4, as well as GARD and PHI tests in SplitsTree showed no evidence of recombination, therefore migration among populations would explain the reticulated pattern.

To explore the possible existence of a molecular clock, TempEst (Rambaut et al., 2016) was employed, and the algorithm found a significant correlation between genetic divergence and time ( Supplementary Figure S4 ). The relevant parameter estimates from the Bayesian analysis were time to most recent common ancestor (TMRCA)=1951.2 ± 2.5 years ago; an effective number of infections of 114.0±6.6; an exponential growth rate of 0.023 ± 0.001 per year; and a rate of molecular evolution of 9.1 ± 0.4×10⁻⁴ substitutions/site/year ( Table 1 ). The plots of the effective number of infections versus time are shown in Figure 1 .

The BaTS algorithm was run to explore whether the MCC clustering was significantly associated to (i) geographic origin, (ii) host species, and (iii) bioclimatic zone. Firstly, using geographic origin as a trait, a highly significant association was found (AI = 8.8954, P <0.0001; PS = 69.3392, P <0.0001) ( Table 2 ). In the single collection sites of Cangzhou (MC = 1.996, P = 0.02), Dali (MC = 1.9891, P = 0.036), Qiqihaer (MC = 1.6233, P = 0.015), Shihezi (MC = 2, P = 0.013), and Hangzhuo (MC = 2, P = 0.005) ( Table 2 ), the isolates were predominantly driving the differentiation. Based on the provinces (AI = 5.99, P <0.0001; PS = 52.3072, P <0.0001) ( Table 3 ), significant differences were found among Anhui, Hebei, Yunnan, Guizhou, Zhejiang, Heilongjiang, Tibet, Qinghai, Ningxia, Shanxi, Shannxi, Inner Mongolia, and Xinjiang. Secondly, using host species (M. sativa, G. pentaphyllum, N. tabacum, N. glutinosa, C. japonica, J. procumbens, V. persica, and T. pretense) as a discrete trait, no significant associations were found ( Table 4 ). Viruses isolated from members of the family Solanaceae clustered together (MC = 2, P = 0.005), while viruses isolated from other families were well mixed in the MCC tree ( Table 5 ). Finally, significant associations were found with bioclimatic zones (AI = 3.9158, P <0.0001; PS = 37.053, P <0.0001) ( Table 6 ). Furthermore, within the subtropical monsoon zone, a significant difference was found among locations (MC = 4.0014, P <0.0035) ( Table 6 ).

Figure 2A and Supplementary Table S5 show that there are quite a lot of sites under purifying selection.

As a complementary approach to the association analyses conducted with BaTS, population genetics AMOVA was also performed. The data were divided into different groups based on different provinces and bioclimatic regions of China ( Supplementary Figure S1 ); the bioclimatic regions of China included four groups—plateau mountain, temperate continental, temperate monsoon, and subtropical monsoon. A population genetics approach was used to further evaluate the geographic population structure suggested by BaTS. BaTS exhibited significant differences among provinces (AI = 5.99, P <0.0001; PS = 52.3072, P <0.0001) ( Table 3 ), and also among localities within province (AI = 8.8954, P <0.0001; PS = 69.3392, P <0.0001) ( Table 2 ), which is similar to the AMOVA results with P = 0.0289 ( Supplementary Table S6 ). Furthermore, a significant difference was observed among different samples taken from the same locality ( Supplementary Tables S6 , S7 ). For instance, in Hangzhou in Zhejiang province, significant variance was detected among the samples (P <0.005) ( Table 2 ). None of the neutrality tests found a significant effect of selection. For the different hosts of AMV, there was a difference only on M. sativa; there was no evidence of differences with the other hosts ( Table 4 ). The samples from four different bioclimatic zones were also analyzed for population structure. The P values indicate that the virus clearly shows differences between bioclimatic regions and among localities within each of the four bioclimatic regions ( Table 5 and Supplementary Table S7 ).

AMV isolates from Iran and Spain were selected for population genetics analyses since they were the most thoroughly sampled isolates. For the Iranian population, the network suggests a single clade with very few divergent genotypes ( Supplementary Figure S5 ), and the pattern is distinct from that of the Chinese network. Recombination among the populations from Iran, Spain, and China was not detected by PHI. The parsimony network MCC tree for isolates from Spain ( Supplementary Figure S6 ) was more similar to the Chinese one, with three well-differentiated groups and the suggestion of recombination between and within the groups by the PHI test. Isolates from Zaragoza province were all over the network. Andalucian isolates were mostly homogenously grouped and the group included an isolate from Zaragoza, suggesting that an introduction from Andalucia may be the origin of the local epidemic in Zaragoza ( Supplementary Figure S6 ).

TempEst found a significant correlation between genetic divergence and time for the Iranian population, allowing for the estimation of molecular evolutionary rates ( Supplementary Figure S7 ). The TMRCA of the Iranian population was 1996.3±1.4 years, with an effective number of infections of 7.4±4.6×10⁶, an exponential growth rate of 0.67±0.04 per year, and a rate of molecular evolution of 8.05±0.50×10⁻⁴ substitutions/site/year ( Table 7 ). These data suggest that the epidemic is older in Iran compared with that in China, and that it expanded at a 10-fold faster rate, with many more infections. For the Spanish population, TempEst found a significant time-stamp in genetic divergence ( Supplementary Figure S9 ), allowing for the estimation of rates of evolution. The TMRCA was 1891.0±5.5 years, with an effective number of infections of 1016.9±73.1, an exponential growth rate of 0.043±0.002 per year, and a rate of molecular evolution of 2.79±0.08×10⁻⁴ substitutions/site/year ( Table 8 ). Therefore, the Spanish epidemic appears to be the most recent ongoing epidemic, followed by China and then Iran. The effective number of infections in Spain is higher compared with that in China, but the epidemic growth rate is slightly slower than in China. Likewise, the rate of molecular evolution of the Spanish epidemic is also the lowest among all three epidemics.

BaTS found evidence of population structure at the geographic level in both Spain and Iran (Iran: AI = 7.192, P = 0.04, PS = 47.730, P = 0.04; Spain: MC = 1.9711, P = 0.014) ( Table 9 ). No evidence of association of isolates by host species or families was found in either case. This again suggests a very recent introduction of the virus and a rapid spread among susceptible hosts without time to adapt to each possible host. For the Iranian population, no differences were detected by AMOVA among provinces or among and within localities ( Supplementary Table S9 ). For the Spanish population, BaTS revealed a significant association between MCC tree clusters and geographic origin (AI = 2.366, P = 0.005), which was mostly driven by isolates from Zaragoza province (MC = 11.7938, P = 0.004) forming their own clusters ( Table 10 ). This maybe indicate the idea that Zaragoza province is the geographic origin of the epidemic. BaTS did not find an association between the MCC tree clusters and the host species or families.

Selection analysis of the Iranian population shows several cases of negative selection and a distinct case of positive selection in Figure 2B and Supplementary Table S10 . For the Spanish population, the selection is predominantly negative ( Figure 2C and Supplementary Table S11 ). Furthermore, using AMOVA, a significant difference was found within Spanish provinces ( P = 0.0026) ( Supplementary Table S12 ) and potentially a few sites with weakly positive selection, but not the same as that suggested for China and Iran.

In the Iranian population, the very rapid spread of the epidemic has not allowed the establishment of significant population differentiation among geographic locations ( Table 1 ); essentially it is a single pandemic population. In the Spanish population, significant differences only exist between samples taken from the same geographic location and not among locations from the same province (likewise the case in China) ( Tables 2 , 3 ). Moreover, based on the host at species and family levels, the Spanish population exhibits significant differences among the hosts ( Supplementary Tables S13 , S14 ), while the Chinese population only shows a significant difference at the host family level ( Tables 4 , 5 ); for the Iranian population, there was no differences among the hosts at species or family levels ( Supplementary Tables S15 , S16 ).