Dear Editor,

The genus Enterovirus (EV), belonging to the family Picornaviridae, order Picornavirales, comprises 15 species: Enterovirus A–L and rhinovirus A–C (Zell et al. 2017). Enterovirus A (EV-A) has 25 serotypes, including several pathogens, such as EV-A71, coxsackievirus A16 (CV-A16), and CV-A6 that infect humans, as well as four serotypes (EV-A122, 123, 124, and 125) isolated from simians (Zell et al. 2017). EV genomes share a similar structure with a length of 7.5 kb, and two open reading frames (ORFs) flanked by 5′- and 3′-untranslated regions (UTRs) in some EVs (Oberste et al. 2013b; Zell et al. 2017; Guo et al. 2019; Lulla et al. 2019).

Although EV-A71, CV-A16, and CV-A6 have been well investigated, the genomic data, epidemiological characteristics of several other enterovirus serotypes remain to be identified (Puenpa et al. 2018; Han et al. 2020a). Other than humans, novel enterovirus serotypes have been identified in non-human primates (NHPs), rodents, dromedary camels, pigs, and bovines (Oberste et al. 2007; Zell et al. 2017; Li et al. 2020). EV-A122 [previously known as Simian enterovirus 19 (SV19)], detected in Macaca fascicularis, showed distant phylogenetic characteristics compared with other EV-A serotypes (Oberste et al. 2007). Surprisingly, certain enterovirus serotypes that circulate in humans have also been detected in NHPs, suggesting the possibility of cross-species transmission (Sadeuh-Mba et al. 2014; Mombo et al. 2015). Enteroviruses identified in NHPs are significant among the emerging enteroviruses that infect humans (Oberste et al. 2013b). Therefore, investigating pathogenic NHP enteroviruses is essential.

In this study, 60 specimens were collected from 30 monkeys, including throat swabs (n = 30, 50%) and rectal swabs (n = 30, 50%). The specimens were collected from several species of monkeys, including Rhesus macaque or Macaca mulatta (n = 16, 53.3%), Macaca thibetana (n = 5, 16.7%), Macaca nemestrina (n = 5, 16.7%), and Macaca assamensis (n = 4, 13.3%), located in two counties of Yunnan Province, China (Supplementary Table S1). We constructed two libraries following the Illumina TruSeq DNA Preparation Protocol, which were pooled and subjected to 150 bp paired-end sequencing using the Novaseq 6000 platform (Illumina, San Diego, CA, USA). We obtained 172, 681, 140 and 176, 787, 228 raw reads from these two libraries (Supplementary Table S1).

The order Caudovirales, which includes viral species infecting bacteria, showed the highest proportion between the two pools, followed by the known eukaryotic viral assignment, Riboviria (Fig. 1A), with 56, 547 and 34, 969 assembled contigs (Trinity version 2.8.4) longer than 300 bp, respectively. Based on annotation against the nucleotide database, we identified 8475 and 1457 contigs in the viral kingdom. Manual analysis resulted in a total of 4 contigs, a nearly complete genome, showing 77% nucleotide identity with the M19s (P2) strain (GenBank accession number AF326754). The clean reads were mapped to the M19s (P2) genome and the sequencing depth was calculated, resulting in two consensus sequences from two pools (Fig. 1B). Based on the full-length genome, both consensus sequences clustered with the RJG7 strain (GenBank accession number EF667344) of the EV-A92 prototype, and also clustered with EV-A species along with the EV-A122 prototype (GenBank accession number AF326754) based on theP1 coding region. (Supplementary Fig. S1A-1B).

Figure 1.  A The taxonomy of viral reads identified in each pool. The y-axis represents the proportion of different viruses. B The sequencing depth and coverage of two pools along with the reference genome. C Similarity plot of novel EV-A122 variants and EV-A122 prototypes. The x-axis shows the nucleotide position, whereas the y-axis represents the genomic similarity between novel EV-A122 variants and the query strain. D Bayesian phylogenetic tree based on the partial VP1 coding region of EV-A122 strains. Different colours indicate the information of countries. The EV-A123 prototype (Genbank no. AF326761) was used as an outgroup. The number at each node represents the posterior probability, and scale bars represent the substitutions per site.

Real-time RT-PCR assays were used to distinguish the positive specimens in the two libraries, with specific primers and probes targeting the assembled genomic sequences (Supplementary Table S2). Finally, we identified nine positive samples with different cycle threshold values in the two pools (Table 1). Among the positive specimens, two throat swabs showed a positive reaction for EV-A122, even though simian enteroviruses are typically detected in stool specimens (Oberste et al. 2007, 2013a). We implemented the ''primer-walking'' strategy to obtain the fulllength genomes (Oberste et al. 2007; Han et al. 2018, 2020b), seven of nine were successfully obtained and were deposited in the GenBank database under accession numbers MT649085–MT649091. The metagenomic data were submitted to NCBI Sequence Read Archive (SRA) under accession number SRP270853.

All EV-A122 variant genomes were 7316–7341 nucleotides long, with a long poly(A) tail. We did not identify the second small ORF (ORF2) in the EV-A122 strains, despite the typical genomic characteristic of several EVs harbouring two ORFs (Supplementary Fig. S2A). Several conserved sites were found in ORF2 of EV-A (Supplementary Fig. S2B). Although ORF2 was not detected in these EV-A122 strains, the conserved domains might play an important role in receptor binding, which determines viral intestinal infection (Guo et al. 2019), the highly conserved WIGHPV domain is required for EV-A71 particle release (Guo et al. 2019; Lulla et al. 2019).

The seven variants showed approximately 79% and 80% nucleotide similarity with the EV-A122 prototype in terms of the VP1 and P1 coding regions, respectively, higher than that of other serotypes (Supplementary Table S3). The seven variants showed the highest nucleotide similarity with the EV-A122 prototype in the P1 coding region and partial 5′-UTR, whereas the similarity between the 3′-UTR and P2 and P3 coding regions was low (Fig. 1C).

The maximum-likelihood phylogenetic tree was constructed using IQ-TREE software, with 1000 bootstrap and SH-like approximate likelihood ratio test (SH-aLRT) replicates. Bayesian inference method (Mrbayes v.3.2.6) was also used to infer phylogenetic relationships. The phylogenetic tree for ORFs indicated that the seven variants clustered with the EV-A92 prototype, and indirectly aggregated with the EV-A122 and EV-A123 prototypes (Supplementary Fig. S2D). The phylogenetic tree of entire P2 and P3 coding regions of the seven variants with the EV-A prototypes showed that these two genogroups clustered with the EV-A92 prototype, implying potential recombination events (Supplementary Fig. S3). Due to limited publicly available genomic sequences, we constructed the evolutionary history of EV-A122 strains according to partial VP1 sequences and identified three lineages; the strains isolated from China in this study and Bangladesh indirectly clustered together, indicating possible circulation of EV-A122 at these locations (Fig. 1D), whereas the strains from the USA clustered together and formed a single lineage (designated as lineage 1). The three lineages showed 22.9%–28.4% nucleotide divergence. These results revealed two different evolutionary directions with two genogroups of isolates in this study. Interestingly, the strains identified from these two genogroups were collected from two different counties of Yunnan Province. Thus, geographic obstruction may play important roles in the evolutionary trajectories of the simian enteroviruses.

Based on the bootscanning plots of the prototypes and variants (Supplementary Fig. S2 and S3), the P1 coding region corresponded to the EV-A122 prototype, whereas the P2 and P3 coding regions showed large-scale recombination activity with EV-A92-like strains. We screened the possible EV-A122-like and EV-A92-like sequences in GenBank with P2 and P3 coding regions. The P1 genomic region of strain cg4006 (EV-A122-like strain), which was identified in the faeces of Rhesus macaque in 2014, showed 77.6%–79% nucleotide identity with that of strains identified in this study (Supplementary Fig. S2–S6). However, the P2 and P3 coding regions of strain CHN/BJ/2018-1A (GenBank accession number MN427525), an EV-A92-like strain collected from faecal specimens of an African green monkey in 2018 in China, were identified as a potential recombinant donors from these limited simian enteroviruses genomes deposited in GenBank, which showed 86.7%–88.5% nucleotide identity with the EV-A122 strains of this study, higher than that of other strains. We thus constructed the phylogenetic trees of the P1 and P2-P3 genomic modules to confirm these recombination events (Supplementary Fig. S7). The EV-A122 variants in this study directly clustered with strain cg4006 in the P1 coding region, whereas they clustered with strain CHN/BJ/2018- 1A in the P2-P3 genomic module. The potential recombination donor strain CHN/BJ/2018-1A belonged to the serotype of EV-A92 in terms of the VP1 coding region. Actually, the significant recombinants identified through this method could be not the real recombinant donors in natural condition, because they might occur more complex genetic exchanges among different serotypes. At least, the strain CHN/BJ/2018-1A harbored the highest phylogenetic identity compared with any other known strains deposited in GenBank. These results demonstrated that Genbank data are sparse regarding simian EVs and that all these simian serotypes have evolved by exchanging genetic materials through inter-serotype recombination. Till present, a largescale recombination event in EV-A122 strains has not been reported previously (Oberste et al. 2007; Harvala et al. 2011; Van Nguyen et al. 2014).

In this study, we did not identify the EV-A92 strains in the specimens of several NHP species in two distant counties of Yunnan Province. The recombination events of EV-A122 strains might have already occurred prior to spreading, emphasising the essentiality of persistent surveillance. Full recombination possibly involves several intermediates of different serotypes, which is yet to be reported. With the increase of serotypes and novel species of enteroviruses detected from several host species, these animals significantly influence the spread of enterovirus disease and cross-species transmission (Oberste et al. 2007; Van Nguyen et al. 2014; Mombo et al. 2015). Enteroviruses, especially serotypes of NHPs, warrant further investigations. A more impeccable and extensive surveillance network is essential for the unabridged enterovirus library, and establishment of necessary surveillance systems for enterovirus and enterovirus diseases is required.