The complete genomic sequences and coding sequence (CDS) of 1006 EV-A strains were obtained from the Virus Pathogen Resource (ViPR) database (Pickett et al., 2012).

Phylogenetic analysis requires that the evolution of a serious of sequences may be properly delineated by a single phylogenetic tree, whereas recombination events in these sequences may severely compromise the accuracy of phylogenetic trees. Hence, non-recombinant datasets produced by recombination analysis are desirable prior to phylogenetic analysis. We have conducted the recombination analysis by Recombination Detection Program (version 4.95) using RDP methods (Martin et al., 2015). Phylogenetic analysis was then carried out on non-recombinant sequences in the MEGA software (version 7.0) using the maximum likelihood model (bootstrap value=1000) (Kumar et al., 2016).

The codon compositions in 3rd position (A3%, U3%, C3% and G3%) have been calculated with Codon W 1.4.2 (Peden and John, 2000). The frequencies of A, U, C and G (%), GC/AU and GC1, GC2, GC12, GC3 content (frequency of occurrence of GC-dinucleotides at each location) have also been calculated. Five codons were excluded from the codon usage bias analyses, including stop codons (UAA, UAG and UGA) and the only codon encoding for Met(AUG) and Trp(UGG).

The frequencies of dinucleotides are known to be influenced by corresponding individual nucleotide abundance. To avoid this bias, the odds ratio index (ρXY) is used (Karlin et al., 1997). X and Y stand for each nucleotide (A, U, G or C). F_X is the frequency of X while F_Y is the frequency of Y. F_XY is the frequency of dinucleotide XY. The calculation of ρXY is as follows:

ρXY=1 indicates that the dinucleotide is present at the corresponding single nucleotide frequencies. Values of ρXY>1 indicate an over-representation of dinucleotide, and ρXY<1 indicates an under-representation.

The relative synonymous codon usage (RSCU) values for all the EV-A CDS were calculated to determine the characteristics of synonymous codon usage. The RSCU values represent the ratio of the observed frequency of one codon to its expected frequency in the synonymous codon family, considering that all codons for an particular amino acid are equally used. Thus, the RSCU values may eliminate the influences from the amino acid compositions or the coding sequence sizes. When the RSCU value is 1.0, there is no bias for codon usage (Sharp and Li, 1986). Synonymous codons with RSCU values >1.0 display a positive codon usage bias while those with RSCU values <1.0 display a negative codon usage bias. Synonymous codons with RSCU values above 1.6 and below 0.6 were considered to be over-represented or under-represented codons, respectively (Wong et al., 2010). The RSCU of individual codons has been calculated by CodonW 1.4.2.

The effective number of codons (ENC) value describes the degree of unbalanced use of synonymous codons (Wright, 1990). ENC values range between 20 and 61. The lower ENC value indicates a greater extent of codon bias within a gene. Typically, an ENC value below 35 denotes a strong codon bias while an ENC value above 50 denotes a weak codon bias. ENC values were calculated through CodonW 1.4.2.

The ENC-GC3 may be used to explore the pattern for synonymous codon usage. Genes where the codon selection is limited only by the G&C mutation will be on or just below the the expected ENC curve (Comeron and Aguadé, 1998).

Correspondence analysis (COA) is a multivariate analysis that provides a geometrical representation of a contingency table, which is used to analyze the major trends of codon usage patterns among genes (Grantham et al., 1980). RSCU is used in COA to minimize the effect of amino acid composition on the analysis of codon usage. Each coding sequence is represented as a vector with 59 dimensions which correspond to the RSCU value of each codon except for Met, Trp and the stop codons. COA reduces the high-dimensional codon-frequency data to a limited number of variables, referred to as principal axes. The axes preserve much of the information about the variability of codon usage among the genes. COA plot selects the first two principal axes as X and Y axes and shows the viral strains as scatters in the two-dimensional diagram. COA analysis is conducted using CodonW 1.4.2.

Correlation analyses were conducted to identify the relationships between the nucleotide composition and the first two COA principal axes using Spearman’s rank correlation test, using the R package “Corrgrams” (Friendly, 2002).

Parity rule 2 (PR2) analysis, exploring the bias between A3/(A3 + U3) and G3/(G3 + C3) in the amino acid family with four-codon, was used to demonstrate the effects of mutation pressure and natural selection on the codon usage. The dots at the center of the plot (A=U, G=C) indicate equal effects of mutation and selection (Sueoka, 1995; Sueoka, 1999).

In the neutrality plot analysis, the regression coefficient (CG12 against GC3) is considered to be the mutation-selection equilibrium coefficient. This analysis was performed to determine the effects of mutation pressure and natural selection on the codon usage patterns of the EV-A coding sequences (Sueoka, 1988). The slope of the regression lines (ϵ) may indicate the evolutionary rates of mutation pressure, whereas 1-ϵ represents selective pressure. If ϵ is close to 1, it means a strong correlation between GC12 and GC3, indicating that mutation pressure is dominant. While ϵ is near 0, it indicates a low mutation pressure.

The codon adaptation index (CAI) is a common measure used to quantify the similarity of the synonymous codon usage between samples and reference. It may also be used to predict the gene expression level based on its codon sequence. The CAI values range from 0 through 1. The higher the CAI value, the higher the degree of gene expression and the higher the adaptation to specific hosts (Sharp and Li, 1987). We also borrowed two other indices reported in recent studies, the relative codon deoptimization index (RCDI) and the similarity index (SiD), to confirm the effects of EV-A codon usage patterns on hosts or vectors. When the RCDI value is closer to 1, the codon usage pattern between the virus and its hosts is more similar. When RCDI is greater than 1, the higher RCDI value indicates greater deoptimization of the virus codon usage patterns relative to hosts (Mueller et al., 2006). The SiD value ranges from 0 to 1. When the SiD value is closer to 1, the codon usage patterns between the virus and the hosts differ significantly, and the effect of host codon usage on the virus is stronger (Zhou et al., 2013).

The CAI and RCDI analysis of the EV-A coding sequences was carried out using the CAIcal server (Puigbò et al., 2008). Here, synonymous codon usage patterns of selected potential hosts were used as references and obtained from the Codon Usage database (Nakamura et al., 2000).

EV-A genomic sequences were obtained from the National Center for Biotechnology Information, including EV71 (KU936132), CVA6 (KX064292), CVA10 (AY421767), CVA16 (MG957117), EVA76 (AY697458.1) and EVA89 (KT277550.1). The viral capsid encoding segment was obtained through in vitro synthesis, and was cloned into the expression vector by a Seamless Cloning Kit (D7010M, Beyotime, Shanghai, China) according to the manufacturer’s instruction. The expression vector was derived from the previously reported “pcDNA6.0A-FY-capsid-GFP” (Chen et al., 2012). The vector is based on pcDNA6.0 backbone and contained a CMV promoter. A fused green fluorescent protein (EGFP) was introduced between the promoter and the viral capsid, and a EVA71 2A protease recognition site (-AITTL-) was inserted between EGFP and VP4. Thus, EGFP was fused to the capsid protein and co-expressed. The amount of green fluorescence could be employed to signal the production of capsid protein.

To avoid the interference by the recombination events in phylogenetic analysis, recombination analysis was first performed on 1006 EV-A genomes to obtain non-recombinant datasets. 878 sequences showing a recombinant signal and 3 containing several stop codons were excluded from the analysis. The P1 region of 125 EV-A genomes and three EVs as outgroup (EVD68, accession numbers: AY426531, CVB3, accession numbers: M33854, PV1, accession numbers: V01150) were subjected to a phylogenetic analysis. A maximum likehood (ML) tree was built, and verified with 1000 bootstrap replicates. EV-A clustered into three phylogenetic clades: clade 1 (n=36), clade 2 (n=80) and clade 3 (n=9). The clade 1 includes CVA7, CVA14, CVA16, EVA71 and EVA120. Clade 2 includes CVA2, CVA3, CVA4, CVA5, CVA6, CVA8, CVA10 and CVA12. Clade 3 comprises EVA76, EVA89, EVA90, EVA91 and EVA92. The descriptions of the 125 strains are given in Supplementary Table S1 and the phylogenetic tree is shown in Supplementary Figure S1 .

To understand the common features and differences of the EV-A genomes, we determined the contents of nucleotides and dinucleotides ( Supplementary Table S2 ). The composition for nucleotide A (28.44 ± 0.35%) is highest, followed by U (24.49 ± 0.32%), C (24.23 ± 0.32%), and G (22.84 ± 0.38%). The proportions of different nucleotides at the 3rd codon position (A3, U3, G3 and C3) show that U3 (29.04 ± 1.04%) and C3 (27.04 ± 0.86%) are higher than A3 (24.65 ± 0.72%) and G3 (19.26 ± 0.90%). Compositions for GC and AU are 47.07 ± 0.61% and 52.93 ± 0.61% respectively, while compositions for GC3 and AU3 are 46.31 ± 1.60% and 53.69 ± 1.59%. These data show that EV-A coding sequences are rich in AU. The Kruskal-Wallis test, one way ANOVA and multiple comparison by LSD test was performed on three EV-A clades for each nucleotide index ( Supplementary Table S3 ). Three of the 15 nucleotide indices show significant differences across the three EV-A clades (p<0.01). Other comparisons show significant differences between clade 1 and clade 2 in 7 out of 15 indices (p<0.05), whereas clade 3 shows significant differences in 15 indices compared to clade 1 and clade 2 (p<0.01).

Analysis of dinucleotide shows that the relative frequency of CpG is significantly lower than other dinucleotides ( Figure 1A ), which explains in part why G and G3 are less frequent. The composition of UpA is also significantly lower than that of other dinucleotides. CpG levels for clade 3 are significantly lower than other clades (p <0.05).

We found that clade 3 has a particular nucleotide composition that may affect codon bias. To estimate the extent of the codon usage bias in EV-A coding sequences, the ENC values were calculated ( Figure 1B ). The average ENC value of all EV-A strains is 55.97 ± 1.43, while those of clade 1, 2 and 3 are 56.44 ± 0.41, 56.30 ± 0.45, 51.16 ± 0.96. For the CDS of each gene, the ENC values of clade 3 were significantly lower than the other two clades (p<0.01), while the ENC values of clade 1 and clade 2 are similar. We also noticed that ENC values of P2P3 regions were significantly lower than P1 (p<0.01). Overall, ENC values indicated that the codon usage bias in the EV-A genomes varied in different clades.

RSCU was calculated to determine the synonymous codons usage in the EV-A coding sequences ( Table 1 ) and the preferred codon with the highest RSCU value for each animo acid was bolded. We noted that the commonly/uncommonly preferred codons ratios among EV-A clade 1:clade 2, clade 1:clade 3 and clade 2:clade 3 were 14:4, 12:6 and 12:6, respectively. By comparing the codon usage patterns of Homo sapiens and EV-A, the ratio of coincident/antagonistic preferred codons is 7/11 in clade 1, 9/9 in clade 2, while only 4/14 in clade 3 among the 18 preferred codons ( Table 1 ). This clade-specific RSCU pattern highlights the distinct evolutionary dynamics of the EV-A clades.

Of the 18 preferred codons, eight are G/C-ended (four C-ended; four G-ended) and the remaining ten are A/U-ended (five A-ended; five U-ended) in the EV-A coding sequences. The findings show that EV-As prefer A/U at the end of codons rather than G/C. We have summarized the RSCU values of the preferred codons by each amino acid in Figure 2 . It also shows that EV-As prefer A/U at the end of codons rather than G/C. We noted that EV-A clade 3 has a stronger AU3 bias than other clades, and the AU3 bias in the P2P3 region of EV-As is stronger than the P1 region.

We built COA plots of CDS and each gene segment to show trends in codon usage variations among different EV-A strains ( Figure 3 ). The first (f’1) and second (f’2) principal axes represent respectively 28.17% and 16.15% of the total variation of the EV-A coding sequences. The EV-A strains were grouped into three distinct clusters on the plot which were shown to be largely correlated with their phylogenetic relationships. Overall, clade 1 and clade 2 formed cluster I, except for CVA5 and CVA6, which formed cluster II; and cluster III comprised all clade 3 strains. These data further support that clade 3 has a unique codon usage pattern. At the level of the individual genes, the clade-specific clustering was also observed ( Supplementary Figure S2 ). For genes encoding 3C, 3D and VP3, clade 3 strains were clustered independently, indicating the effect of divergence events on individual genes. For the genes encoding 2A, 2BC, 3AB, VP1, VP4 and VP2, overlapping patterns in codon usage were observed. Based on these data and together with ENC and RSCU, we confirm that clade 3 has a unique codon usage pattern and likely has an independent evolutionary origin. Conversely, the overlap between clade 1 and clade 2 indicates that they have similar codon usage patterns and a closer relationship. Together, clade 3 stands apart from other clades, which may be affected by unique evolutionary forces.

Mutational pressure and translational selection were reported as the primary factors influencing the codon usage bias. Here, we examined whether they have affected the EV-A codon usage and which factor is dominant through correlation analysis, PR2 plot and ENC-GC3 plot analysis. In the correlation analysis, we compared the correlation between the nucleotide compositions and the principal axes of the COA plot, and noted a significant positive or negative correlation between nucleotide composition, such as A/A3 or A/G. We also observed that most nucleotide composition constraints were significantly correlated with f’1 and f’2 ( Figure 4 ). Together, these results supports that mutation pressure from the nucleotide composition may have affected the EV-A codon biases.

In PR2 plot, the relationship between A/U and G/C composition on the 3rd codon position in four degenerated codon families (Ala, Arg, Gly, Leu, Pro, Ser, Thr, and Val) was analyzed. We observed that U3 was used more frequently than A3, while G3 is approximately equal to C3 in the entire coding sequences ( Figure 5 ). G3 and U3 were more frequent in P1 regions, while C3 and U3 were more frequent in P2-P3 regions. The A3/U3 and G3/C3 bias denotes the influence of both selective pressure and mutation effect.

To understand whether pressure dominates in the EV-A codon usage bias, ENC–GC3 plots were constructed. We observed that all the EV-A strains clustered together under the expected ENC curve ( Figure 6A ). This suggests that natural selection may affect the genomic evolution of EV-A strains. The EV-A clades deviated from the expected curve differently and the clade 3 strains were lower below the curve, suggesting that clade 3 strains are subjected to higher selective pressure than other clades. The clade-specific difference can also be seen in P1 and P2P3 coding regions ( Figures 6B, C ). However, the effects of mutation pressure and natural selection on individual proteins varied ( Supplementary Figure S3 ). All 2BC and 3D gene data points fell below the expected curve, while the 3AB, 3C, VP3 and VP4-VP2 coding sequences of some strains were aggregated on the curve, suggesting a dominant influence due to mutation pressure. It is noteworthy that clade 3 strains are the lower points under the curve in 3AB, 3D and VP3 plots, indicating a greater effect by natural selection in these coding regions.

Mutation pressure and natural selection help to shape the codon usage patterns. To study the magnitude of each force, we used the neutrality plot analysis. We observed a significant positive correlation between GC12 and GC3 values across the EV-A coding sequences (slope =0.113, r =0.695, p=0.007). The linear regression slope suggested that the relative neutrality (mutation pressure) was 11.31%, and the relative constraint on GC3 (natural selection) was 89.69%, indicating that the EV-A codon usage patterns were primarily shaped by natural selection. Regression slopes in clade 1 and clade 2 are 0.1461 and 0.0417, while the correlation is not significant in clade 3, consistently indicating high levels of influence on codon patterns by natural selection ( Figure 7A ). The results of the neutrality plot analysis in P1 and P2P3 regions were similar to the whole coding sequences ( Figures 7B, C ). We analyzed individual proteins in three clades, and observed significant correlations between the GC12 and GC3 values of the VP1, 2A, 3D sequences in clade 1 and the VP3, 2BC, 3D coding sequences in clade 2 strains. All slope values for each protein are below 0.5 (from 0.007 to 0.4076), which supports the stronger effects of selection pressure than mutation pressure ( Supplementary Figure S4 ). Notably, the VP1 coding sequences in clade 1 showed an exceptional slope value= 0.8696, indicating a weak selection force in this region ( Figure 7D ). Overall, our data suggest that the influence of natural selection prevails over the coding sequences of EV-A strains.

The CAI analysis reflects the fitness of the viral gene for the host cell. The CAI values were calculated for each coding sequence against 38 model organisms based on the CUTG database. The six highest average CAI values obtained for the entire EV-A coding sequence were assigned to Xenopus laevis (0.82 ± 0.006), Ciona intestinalis (0.78 ± 0.009), Gallus gallus (0.77 ± 0.005), Danio rerio (0.76 ± 0.005), Homo sapiens (0.74 ± 0.005), and Mus musculus (0.74 ± 0.005) ( Figure 8 and Supplementary Table S4A ). For all six matched species, the CAI analyses of the P1 and P2P3 regions are consistent with the entire coding sequence. We also observed that the CAI values of clade 3 is significantly higher than clade 1 and clade 2 in Xenopus laevis and Ciona intestinalis (p<0.01), indicating that EV-A strains show clade-specific codon adaptation patterns.

In addition, we also calculate the RCDI and SiD values to show the degree of unsuitability and potential impact of model organism codon usage patterns on EV-A. ( Supplementary Figure S5 , and Supplementary Table S4 ). The six lowest average RCDI values of the entire EV-A coding sequence include Xenopus laevis (1.06), Ciona interstinals (1.09), Gallus gallus (1.11), Danio rerio (1.11), Mus musculus (1.12), Homo sapiens (1.12). Among the 38 model animals, the six species with the lowest mean SiD values also matched the CAI and RCDI analysis, indicating that they have similar codon bias patterns and may be most appropriate hosts for EV-A. The data of the P1 and P2P3 regions are consistent with the entire coding sequence.

As above, we concluded that there are significant differences in codon bias patterns between conventional and unconventional EV-A strains. To explore whether the codon bias can affect the protein translation of EV-A, we first tested if there are differences in protein translation between viruses representative of three clades. The expression of GFP-tagged viral capsid protein in HEK293T cells was examined by fluorescence microscopy ( Figures 9A, B ). It was found that EV71 and CVA16 of clade 1, CVA6 and CVA10 of clade 2 had similar levels of expression, significantly higher than EVA76 and EVA89 of clade 3. This result indicates that the expression levels of conventional EV-A capsids are substantially higher than those of unconventional EV-A.

To confirm if the lower expression levels of unconventional clades are caused by codon bias, we have optimized the codon of the capsid region of EVA89 and EVA76 according to Homo sapiens. The optimized fragments EVA76-Cap-opti and EVA89-Cap-opti were expressed in HEK293T and RD cell lines. Using immunofluorescence microscopy and Western blot, we found that the optimized capsid protein of EVA89 and EVA76 expressed at significantly higher levels ( Figure 9C ), showing that the codon bias of clade 3 does contribute to the lower protein expression.