Potato virus A (PVA) isolates were collected from major potato-growing areas in China. Each isolate was maintained a plant of Nicotiana debneyi in the lab. The presence of PVA was confirmed by DAS-ELISA (Agdia, Elkhart, United States). Total RNA was extracted from each N. debneyi sample using Trizol (Invitrogen, Carlsbad, CA, United States) and reverse-transcribed following the manufacturer’s protocol (Promega, Madison, WI, United States). The complete genome of PVA was obtained by amplifying 10 overlapping fragments (nucleotides 1–409, 200–1428, 1259–2530, 2384–3630, 3476–4730, 4572–5824, 5725–7040, 6889–8160, 7963–9345, and 8986–9567) using 10 pairs of degenerate primers (Supplementary Table 1), which were designed from highly conserved regions of published PVA genomes (accession numbers Nos. KF977085, MT521081, MT435487, MT435489, KM365068, MT502380, and MT502370).

Polymerase chain reaction (PCR) amplifications of cDNA were performed in a total volume of 50.00 μL, containing 2.00 μL of cDNA templates, 25.00 μL of Premix Taq (TaKaRa), 1.00 μL of forward primer (10.00 μmol/L), 1.00 μL of reverse primer (10.00 μmol/L), and 21.00 μL of ddH₂O. The PCR program comprised 5 min at 94°C; 35 cycles of 94°C for 30 s, 48°C–57°C for 30 s (Supplementary Table 1), and 72°C for 1 min; followed by a final extension of 10 min at 72°C. Samples were amplified using a DNA Engine Peltier Thermal Cycler (Bio-Rad Laboratories, Hercules, CA, United States). The PCR products were separated on 1.5% agarose gels in Tris-acetate-EDTA (TAE) buffer and visualized by UV illumination.

PCR products were purified and ligated to pESI-T vector, which was provided in the Hieff Clone ^® Zero TOPO-TA Cloning Kit (Yeasen, China), and propagated in cells of Escherichia coli strain TOP10. The cloned DNA fragments of recombinant plasmids were sequenced in both directions by Sangon Biological Co. Ltd (Shanghai, China). At least 3–5 independent cDNA clones for each segment were sequenced to assemble consensus sequences.

Eleven complete or nearly complete genome sequences of PVA isolates, including nine from China, one from Peru, and one from the Netherlands, were obtained in this study and deposited in GenBank under accession numbers MW592838–MW592842 and MW616801–MW616806 (see Supplementary Table 1 for the list of primers used for sequencing). In addition to the novel sequences, 55 complete genome sequences of PVA isolates were downloaded from GenBank (Supplementary Table 2). The sequences had been collected from 14 countries and had known host species and geographic locations. To increase post-analysis interpretability, the isolates were grouped according to their geographic origins. The combined sequence data included China (n = 10), Europe (n = 15), South America (n = 33), Oceania (n = 6), and North America (n = 2) and were used for the subsequent analyses. Sequences were aligned using the MEGA X (Kumar et al., 2018) and the polyprotein ORF of each sequence was extracted from the alignment. Codon-based sequence alignment was then performed using the MAFFT algorithm (Katoh and Standley, 2013) implemented in PhyloSuite v1.2.2 (Zhang et al., 2020). The program was run using the FFT-NS-I iterative refinement method with the following parameters: mafft –thread 8 –threadtb 5 – threadit 0 –reorder –leavegappyregion. Ambiguously aligned regions were trimmed using the program Gblock 0.91b (Talavera and Castresana, 2007) implemented in PhyloSuite, with the “codon” mode, half gaps allowed, and all other parameters at default settings. The resulting sequence alignment had a length of 9180 nucleotides and used for subsequent population genetics analysis.

To assess how the diversity varied across geographical and host populations, haplotype diversity (H_d) and nucleotide diversity (P_i) were calculated using DnaSP v5.0 (Librado and Rozas, 2009). Analysis of molecular variance (AMOVA) was also carried out using Arlequin v3.5 (Excoffier and Lischer, 2010). The significance of φ-statistics was tested by 1023 random permutations of sequences among the population.

Pairwise among-populations fixation indices (F_ST) were calculated using Arlequin v3.5 (Excoffier and Lischer, 2010), and the significance was obtained with 1000 permutations. A sliding-window analysis was used as an additional approach for evaluating genetic population differentiation. This analysis was performed using the PopGenome package (Ver. 2.7.5; Pfeifer et al., 2014) in R software (ver. 3.5.1), with a window size of 100 nt and a step size of 30 nt. In addition, discriminant analysis of principal components (DAPC) was used to infer clusters of genetically related individuals. This new multivariate method pioneered by Jombart et al. (2010) was designed to investigate the genetic structure of biological populations without assuming panmixia. In this study, we only performed the DAPC analysis based on pre-defined geographic groups using the adegenet package (Ver. 2.0.1; Jombart, 2008) in R software (ver. 3.5.1), and therefore, the populations of North America (two potato isolates) and Oceania (three tamarillo isolates) were excluded from the analysis due to an inadequate sample size (n ≤ 3).

A recent study found evidence for intragenic recombination within the PVA genome (Fuentes et al., 2021). To investigate the role of recombination, we used two different methods to investigate the occurrence of recombination events in CP sequences. A phylogenetic network was first reconstructed using the neighbor-net method implemented in SplitsTree v4.14.8 (Huson, 1998) with default settings. The pairwise homoplasty index (PHI) test implemented in SplitsTree was also carried out to test signals of recombination (p < 0.05, significant evidence of recombination). In addition, to confirm the occurrence of recombination in our dataset, localization of recombination breakpoints and identification of likely parental sequences were achieved with the RDP v4.101 package, which incorporates the algorithms RDP, Geneconv, Bootscan, Maxchi, Chimaera, Siscan, and 3Seq (Martin et al., 2015). Recombination events supported by at least four different algorithms of analysis and with p values < 1.0 × 10^–5, viral isolates were identified as recombinants. Because recombinants may result in misleading results in selection analysis, as reported by Anisimova et al. (2003) and Sironi et al. (2015), the recombinants were excluded from subsequent selection analysis.

To measure the selection pressure in the complete PVA genome, we calculated the ratio (ω) of non-synonymous (dN) to synonymous (dS) substitutions, as it was done in most adaptive evolution studies, using the CodeML program in the PAML package (Yang, 2007) implemented in EasyCodeML v1.4 (Gao et al., 2019). After all recombinants were removed, our selection analysis was based on 58 polyprotein coding-region sequences. The positive selection models (M2a, M8) and their respective null models (M1a, M7) implemented in the site models were used to conduct the adaptive evolution analysis. Likelihood ratio tests (LRTs) were performed twice to compare the difference in the log-likelihoods between the nested codon-based models against an x² distribution with the degree of freedom equal to the differences in the number of parameters between the models (Yang, 1998). When the LRTs yielded significant results, the Bayes empirical Bayes (BEB) method was used to identify the codons that were the most likely to be under positive selection (Yang et al., 2005).

To assess the geographical and host effects on the PVA population, three statistics, the association index (AI), parsimony score (PS), and maximum monophyletic clade size (MC), were calculated from the posterior tree samples using BaTS v2.0 (Parker et al., 2008). For this analysis, phylogenetic uncertainty was used to investigate phylogeny-trait correlations, with 1000 random permutations of tip locations to estimate a null distribution for each statistic. The results that generated a low AI index and PS and high MC scores with p < 0.05 suggested a strong phylogeny–trait association.

A data set consisting of 66 complete sequences was included in the analysis. After trimming the ambiguously regions from the alignment, we found that all mutations in the PVA genome are substations, with two sequences as exceptions, which had one (accession number: AJ131403) and two codon deletions (accession number: MT502380), respectively. The 66 PVA isolates in this study were composed of 66 haplotypes with an overall haplotype diversity of 1.00 and nucleotide diversity of 0.077 (Table 1). When the viral isolates were categorized according to geographic origin, the highest nucleotide diversity (0.103 ± 0.019) was found in the Oceania population and the lowest (0.012 ± 0.002) was discovered in the Chinese population. When the isolates were grouped according to host origin, higher genetic variation was observed in viral isolates collected from tamarillo compared with those isolated from potato (Table 1), although only four isolates isolated from tamarillo were analyzed in this study. Stepwise diversity analysis indicated that the fragment spanning nucleotide 1-189 and 2791-2950 (Supplementary Figure 1), which corresponds to the coding region for P1 and PIPO, respectively, is the most variable and conserved region on the genome of PVA. Sixty haplotypes were identified in the 66 nucleotide sequences of P1 with an overall haplotype diversity of 0.994 and nucleotide diversity of 0.104 (Supplementary Table 3). In PIPO, 28 haplotypes were identified in the 66 sequences, with an overall haplotype diversity of 0.876 and nucleotide diversity of 0.025 (Supplementary Table 3).

The genetic differentiation between all populations of geographic regions was significant, with the F_ST values ranging from 0.22 to 0.46 (Table 2), indicating significant genetic differentiation between geographic groups of PVA. Similarly, the genetic differentiation between viral isolates with different host species had an F_ST value of 0.45. The results of sliding-window analysis of the pairwise F_ST values of population differentiation among geographic regions and host species are illustrated in Figure 1. The pairwise F_ST values estimated based on the geographic groupings were similar to those based on the host species groupings. This was in agreement with the pairwise F_ST analysis above (Table 2).

An AMOVA analysis also revealed a significant level of genetic differentiation between the PVA sequences originating from either different geographic origins or host species. As shown in Table 3, the variation among geographic regions accounted for 28.07% of the total variation (Φ_ST = 0.281, p < 0.001), while the variation within regions accounted for 71.93%. When performing the AMOVA only on viral isolates from potato and tamarillo, similar results were obtained; significant variation among groups made up 45.81% of the total variation (Φ_ST = 0.458, p < 0.01, Table 3), which accounted for nearly 50% of the total genetic variance of PVA. Taken together, it seems that the effect of host species on the genetic variance of PVA is greater than that of geography.

The results of the DAPC analysis showed similar patterns of population differentiation as those revealed by the pairwise F_ST analysis. Discriminant analysis of principal components scatter plots indicated that the population of China were relatively distinct from the other populations along the first discriminant function axis, while the population of Europe exhibited more subtle structure along the second discriminant function axis (Figure 2). Discriminant analysis of principal components scatterplots also showed that the analyzed PVA isolates were divided into three genetic clusters (Figure 2), corresponding to the geographic regions. All the three genetic clusters were clearly differentiated. Cluster 1 contained all PVA isolates from China. Cluster 2 contained 15 individuals from Europe, whereas Cluster 3 comprised 33 individuals from South America.

Our phylogenetic network analysis showed that PVA isolates were clustered into three lineages (Figure 3), corresponding to the groups W (World), A (Andean), and T (Tamarillo) from the phylogenetic analysis by Fuentes et al. (2021). Three isolates from tamarillo were placed into the T lineage and 14 isolates were placed into the A lineage. The W lineage contained isolates with a considerable diversity of sampling locations. Within the W lineage, all Chinese isolates formed a highly divergent sub-lineage (Figure 3). We also found several conflicting phylogenetic signals that may have been due to recombination (Figure 3), which was supported by the PHI test with statistically significant evidence of recombination (p < 0.001). Using the RDP package, 8 PVA isolates were identified as recombinants by at least 4 algorithms (Table 4). In one recombinant (Apu046, accession number MT502353), the breakpoints were detected in the P3 cistron, whereas in all other 7 recombinants (accession numbers GU144321, MT435486, MT435487, MT435489, MT435495, MT502353, MT502377 and MT521083), the breakpoints were detected in the CP cistron (Table 4). The sequences of these eight recombinants were excluded from the selection analysis presented below.

Fifty-eight non-recombinant sequences were included in the selection and phylogeny-trait association analyses. The ratio of mean dN/dS (less than 1) of the polyprotein coding region showed that the majority of polymorphic sites (98.85%) were under purifying selection (Figure 4), suggesting that most of mutations in the genome were deleterious and consequently being weeded out by natural selection. However, the LRT indicated that the positive selection models (M2a and M8) were significantly better than the control models (M1a and M7), providing evidence for the presence of codons under positive selection. Further analysis from BEB scores indicated a strong positive selection pressure on nine codons, including the P1 (codon sites 34, 46, and 146), NIa (codon site 2268), and NIb cistrons (codon sites 2557, 2560, 2561, 2563, and 2591), with high posterior probability ≥0.95 (Table 5).

With the exception of viral isolates from Europe (MC = 2.00, p > 0.05), significant signal for geographic clustering was found when PVA isolates were grouped by their sampling regions based on tests of phylogeny-trait association analysis (MC: p < 0.05, Table 6), indicting a great spatial structure of the pathogen. However, we accepted the null hypothesis of no association between host species and phylogenetic relationships when the PVA isolates were grouped by their host origins (MC: p > 0.05, Table 6). Taken together, the BaTS results suggested that geographic effects contributed to the diversification of the virus, which could be explained by geography-driven adaptation.