# Genomic surveillance of enterovirus associated with aseptic meningitis cases in southern Spain, 2015–2018

Fabiana Gámbaro, Ana Belén Pérez, Eduardo Agüera, Matthieu Prot, Luis Martínez-Martínez, María Cabrerizo, Etienne Simon-Loriere, Maria Dolores Fernandez-Garcia

► To cite this version:

Fabiana Gámbaro, Ana Belén Pérez, Eduardo Agüera, Matthieu Prot, Luis Martínez-Martínez, et al. Genomic surveillance of enterovirus associated with aseptic meningitis cases in southern Spain, 2015–2018. Scientific Reports, 2021, 11 (1), pp.21523. 10.1038/s41598-021-01053-4. pasteur-03720850

HAL Id: pasteur-03720850

https://pasteur.hal.science/pasteur-03720850

Submitted on 12 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Genomic surveillance of enterovirus associated with aseptic meningitis cases in southern Spain, 2015–2018

Fabiana Gámbaro<fnref n="1" />,⁹, Ana Belén Pérez²,³,⁹, Eduardo Agüera²,³, Matthieu Prot<fnref n="1" />, Luis Martínez-Martínez²,³,⁴, María Cabrerizo⁵,⁶,⁷, Etienne Simon-Loriere<fnref n="1" />,¹⁰& Maria Dolores Fernandez-Garcia³,⁸,¹⁰

New circulating Enterovirus (EV) strains often emerge through recombination. Upsurges of recombinant non-polio enteroviruses (NPEVs) associated with neurologic manifestations such as EVA71 or Echovirus 30 (E30) are a growing public health concern in Europe. Only a few complete genomes of EVs circulating in Spain are available in public databases, making it difficult to address the emergence of recombinant EVs, understand their evolutionary relatedness and the possible implication in human disease. We have used metagenomic (untargeted) NGS to generate full-length EV genomes from CSF samples of EV-positive aseptic meningitis cases in Southern Spain between 2015 and 2018. Our analyses reveal the co-circulation of multiple Enterovirus B (EV-B) types (E6, E11, E13 and E30), including a novel E13 recombinant form. We observed a genetic turnover where emergent lineages (C1 for E6 and I [tentatively proposed in this study] for E30) replaced previous lineages circulating in Spain, some concomitant with outbreaks in other parts of Europe. Metagenomic sequencing provides an effective approach for the analysis of EV genomes directly from PCR-positive CSF samples. The detection of a novel, disease-associated, recombinant form emphasizes the importance of genomic surveillance to monitor spread and evolution of EVs.

Enteroviruses (EVs) belong to the family *Picornaviridae*, genus *Enterovirus*. Their positive single-strand RNA genome is about 7,500 nucleotides (nt), and is composed of a large open reading frame (ORF) flanked by 5' and 3' untranslated regions (UTRs)<fnref n="1" />. The 5' part of the ORF encodes the structural proteins that form the capsid (among them the VP1 which is the most external), while the 3' part of the ORF encodes the non-structural proteins. There are more than 100 types of EVs infecting humans which are classified into four species A to D². EVs are associated with a wide spectrum of clinical symptoms ranging from nonspecific febrile illness or upper respiratory illness to severe neurological conditions, including meningitis, encephalitis or acute flaccid paralysis<fnref n="1" />.

Genetic recombination is a major process in EVs evolution³. Since new recombinant EVs may present different properties (variations in virulence or transmissibility), it is essential to capture complete genomes to detect the occurrence of these recombination events for surveillance and public health purposes. Application of metagenomic next-generation sequencing (mNGS) approaches to characterize the whole-genome sequence of EVs has become a valuable tool for detecting multiple EVs that may be present during co-infection, especially caused by two EV types of same species. Currently, Sanger sequencing of the VP1 capsid protein gene is the gold standard for EV genotyping. But as mNGS becomes increasingly affordable, accessible and cost-effective, sequencing the whole genome of EVs will likely allow improved typing in the near future⁴. Therefore, it is crucial that microbiology reference laboratories progress towards implementation of NGS methods to fully characterize and respond to EVs and publicly share these sequence data⁴. In most countries in Europe, non-polio EV (NPEV) infections are not notifiable and surveillance is mainly passive⁵. Because there is no specific EV surveillance system, surveillance captures primarily EVs detected from hospitalized patients with neurological infections⁵. Identification and full

<footnote n="1">¹Institut Pasteur, Paris, France. ²Hospital Universitario Reina Sofía, Córdoba, Spain. ³Instituto Maimónides de Investigación Biomédica de Córdoba (IMIBIC), Córdoba, Spain. ⁴Universidad de Córdoba, Córdoba, Spain. ⁵National Centre for Microbiology, Instituto de Salud Carlos III, Madrid, Spain. ⁶CIBER de epidemiología y Salud Pública (CIBERESP), Madrid, Spain. ⁷Red de Investigación Translacional en Infectología Pediátrica (RITIP), IdiPaz, Madrid, Spain. ⁸Present address: National Centre for Microbiology, Instituto de Salud Carlos III, Madrid, Spain. ⁹These authors contributed equally: Fabiana Gámbaro and Ana Belén Pérez. ¹⁰These authors jointly supervised this work: Etienne Simon-Loriere and María Dolores Fernández-García. & email: etienne.simon-loriere@pasteur.fr; mdfernandez@isciii.es</footnote>
<table><caption>**Table 1.** Clinical features and laboratory findings for the cases of enteroviral meningitis, Córdoba, Spain, 2015–2018. EV enterovirus, CSF cerebrospinal fluid, GI gastrointestinal symptoms, Gluco glucose, NA not available, Prot proteins, WBC White Blood Cells.</caption><thead><tr><th rowspan="2">EV type</th><th rowspan="2">Sample N°</th><th rowspan="2">Age (y)</th><th rowspan="2">Sex</th><th rowspan="2">Clinical features</th><th rowspan="2">Date sample collection</th><th rowspan="2">Days of hospitalization</th><th colspan="4">CSF analysis</th></tr><tr><th>WBC mm³</th><th>Lymph %</th><th>Gluco (mg/dl)</th><th>Prot (mg/dl)</th></tr></thead><tbody><tr><td rowspan="4">E6</td><td>138</td><td>19</td><td>F</td><td>Fever, headache, vomiting, neck stiffness, GI</td><td>May-2015</td><td>2</td><td>300</td><td>48</td><td>94</td><td>41</td></tr><tr><td>255</td><td>31</td><td>M</td><td>Fever, headache, nausea, GI</td><td>Apr-2015</td><td>4</td><td>233</td><td>42</td><td>59</td><td>69</td></tr><tr><td>268</td><td>27</td><td>F</td><td>Fever, headache, nausea, vomiting</td><td>Apr-2015</td><td>6</td><td>118</td><td>63</td><td>55</td><td>46</td></tr><tr><td>365</td><td>23</td><td>F</td><td>Fever, headache</td><td>May-2015</td><td>2</td><td>203</td><td>61</td><td>63</td><td>41</td></tr><tr><td rowspan="3">E11</td><td>484</td><td>29</td><td>F</td><td>Fever, headache, nausea</td><td>Nov-2015</td><td>2</td><td>86</td><td>90</td><td>53</td><td>47</td></tr><tr><td>1059</td><td>29</td><td>M</td><td>Fever, headache, photophobia, GI</td><td>Dic-2018</td><td>3</td><td>383</td><td>72</td><td>NA</td><td>68</td></tr><tr><td>1106</td><td>36</td><td>F</td><td>Fever, headache</td><td>Dic-2018</td><td>2</td><td>312</td><td>90</td><td>55</td><td>63</td></tr><tr><td>E13</td><td>53</td><td>35</td><td>F</td><td>Fever, headache, nausea, GI</td><td>May-2016</td><td>10</td><td>812</td><td>72</td><td>50</td><td>113</td></tr><tr><td rowspan="4">E30</td><td>265</td><td>30</td><td>F</td><td>Fever, headache, neck stiffness</td><td>Feb-2016</td><td>1</td><td>101</td><td>90</td><td>53</td><td>45</td></tr><tr><td>519</td><td>31</td><td>M</td><td>Fever, headache, vomiting</td><td>Mar-2018</td><td>2</td><td>422</td><td>91</td><td>NA</td><td>63</td></tr><tr><td>520</td><td>25</td><td>M</td><td>Headache, vomiting, photophobia, GI</td><td>May-2018</td><td>4</td><td>310</td><td>82</td><td>NA</td><td>41</td></tr><tr><td>675</td><td>17</td><td>M</td><td>Fever, headache, vomiting, photophobia</td><td>Mar-2017</td><td>3</td><td>156</td><td>63</td><td>60</td><td>58</td></tr></tbody></table>

characterization of the EV types involved in these patients are essential to (1) identify EV types associated with neurological infections and estimate better their associated disease burden, (2) monitor the emergence of new EV strains or new recombinant forms, and (3) gain a better understanding of circulating NPEVs⁴.

Despite the common occurrence of EV-associated meningitis cases in Spain⁶⁻⁸, there are limited studies defining their molecular epidemiology and addressing the emergence of recombinant EVs. These studies have sequenced either multiple parts of the genome (the complete VP1 and 3Dpol proteins)⁹⁻¹⁵ or the whole genome¹⁶,¹⁷. However, very few full genomes of EVs circulating in Spain apart from EVA71 or EVD68 are available in public databases (Echovirus 30 [n = 5], Coxsackie B2 [n = 1] and Coxsackie A6 [n = 3]) making it difficult to understand the evolutionary relatedness of the different EV types and the possible implication in human disease.

Here, we aimed to characterize the genetic diversity of EVs circulating in Southern Spain during 2015–2018. We performed mNGS on 12 CSF samples of laboratory-diagnosed enteroviral meningitis. Our viral genomic analysis revealed the co-circulation of E6, E11, E13 and E30 types including a novel E13 recombinant form.

## Results

**Clinical, epidemiological and laboratory features of patients with enteroviral meningitis.** From 2015 to 2018, twelve adult patients were laboratory-diagnosed with enteroviral meningitis at the University Hospital Reina Sofía (Córdoba, Spain). All samples were identified as EVs positive but no other typing data was available for these strains. In this work, we investigated the twelve CSF samples collected from these patients by using ribosomal RNA depletion before mNGS on total extracted RNAs. Epidemiological and clinical information about each patient is listed in Table 1. The median age of the patients was 27.6 years (range 17–36 years) with a female-to-male ratio of 1:0.7. All patients were living in the Cordoba Province, located in southern Spain. Cases showed a seasonal pattern, with E6, E13 and E30 infections occurring during spring (March–May) months, while E11 cases occurred during winter (November and December). The most common symptoms in patients were headache (100%) and fever (90.9%). Seven patients (58.3%) reported contact with potential EV-infected children at home that presented fever or gastrointestinal symptoms. All patients had pleocytosis, defined as CSF WBC counts of > 5/mm³, with a median WBC count of 286/mm³ and predominance of lymphocytes (from 42 to 91% of WBC) (Table 1).

**Complete genome analysis.** After quality trimming, we assembled the reads using MetaSPAdes v3.12, and used Diamond to query the contigs against a non-redundant protein database (NCBI). For each sample, we obtained a large contig corresponding to a near-complete enterovirus genome, or fragments that could be assembled into a genomic scaffold. We did not detect co-infection with other pathogens, including those that have been found in the CSF of patients with meningoencephalitis such as Dengue virus, Chikungunya virus or West Nile virus¹⁸. Results of the metagenomic analysis can be found in Table 2. De novo assembly resulted in the reconstruction of complete or near-complete EVs genome for all 12 samples. All viruses presented high sequence coverage throughout the genome as shown in Supplementary Fig. S1. Overall, 4 different EV types belonging to species-B were identified. Four samples were found to contain E6, three samples E11, four samples
<table><caption>**Table 2.** RNA-sequencing results using mNGS from CSF samples.</caption><thead><tr><th></th><th>Sample N°</th><th>Number of mapped viral reads</th><th>Mapped viral reads (%)</th><th>Genome covered</th><th>Average coverage depth</th></tr></thead><tbody><tr><td rowspan="4">E6</td><td>138</td><td>127,810</td><td>1.2</td><td>99.30%</td><td>1273</td></tr><tr><td>255</td><td>70,030</td><td>0.6</td><td>99.88%</td><td>697</td></tr><tr><td>268</td><td>56,357</td><td>0.59</td><td>99.86%</td><td>560</td></tr><tr><td>365</td><td>6524</td><td>0.14</td><td>99.31%</td><td>64</td></tr><tr><td rowspan="3">E11</td><td>484</td><td>1046</td><td>0.01</td><td>95.76%</td><td>10</td></tr><tr><td>1059</td><td>240,901</td><td>2.76</td><td>99.40%</td><td>2383</td></tr><tr><td>1106</td><td>366,776</td><td>3.13</td><td>99.70%</td><td>3624</td></tr><tr><td rowspan="2">E13</td><td>53</td><td>159,567</td><td>1.61</td><td>100%</td><td>1593</td></tr><tr><td>265</td><td>68,949</td><td>0.35</td><td>98.91%</td><td>684</td></tr><tr><td rowspan="3">E30</td><td>519</td><td>54,571</td><td>0.33</td><td>99.30%</td><td>537</td></tr><tr><td>520</td><td>628,623</td><td>4.39</td><td>99.60%</td><td>6214</td></tr><tr><td>675</td><td>91,887</td><td>0.76</td><td>99.50%</td><td>909</td></tr></tbody></table>

E30 and one sample E13. Genomes ranged from 7213 to 7451 nt in length, encoding polyproteins of 2188, 2191, 2194 and 2195 aa for E13, E6, E30 and E11, respectively.

**Phylogenetic analysis.** In order to characterize the study strains and evaluate their evolutionary relationship to previously characterized homotypic strains detected in Spain and globally in the last decades, phylogenetic analyses were performed for the four encountered genotypes first using the VP1 region (Fig. 1).

E6 lineages were assigned according to Smura et al. and Cabrerizo et al. into lineages A-C<sup>10,19</sup>. Phylogenetic analysis showed that all E6 strains detected in 2015 clustered within sublineage C1 (bootstrap value of 97%). Previously reported E6 strains from Spain (detected between 2004 and 2010) fell into the sublineage C9, clustering into three distinct groups (C9a, C9f, and C9h) (Fig. 1A). All four E6 study strains showed together high nucleotide similarity in the complete VP1 region (mean 97.8% (96.2–99.2%) nt and 99.2% (98.6–100%) aa identity) and branched together within the diversity of strains detected in France and Poland in 2011.

E30 lineages were defined according to Bailly et al. into lineages A-H<sup>20</sup>. The E30 study strains clustered in two different lineages. Sequences LCR265 (2016) and LCR675 (2017) clustered within lineage F, a group that contained most of E30 strains previously identified in Spain (from 1996 to 2018). Analysis with complete VP1 sequences (Supplementary Fig. S2) showed that while LCR675 strain forms a subcluster (bootstrap value of 100%) with other E30 strains from Spain isolated in 2016 in Catalonia (displaying identities of at least 96.4% and 91.1% at nucleotide and peptide levels, respectively), LCR265 cluster with viruses detected in France, Luxembourg and the Netherlands in 2016–2017 suggesting the circulation of two different sublineages. In contrast, strains LCR519 and LCR520 from 2018 grouped into a distinguishable lineage, tentatively called here “lineage I” (corresponding to clade G6 in the classification of Benschop et al.<sup>21</sup> which contains sequences from Netherlands, France, Ireland and Spain, sampled in the same time period (2017–2018) (Fig. 1B). Similar pattern of clustering was obtained when the complete VP1 region was used for phylogenetic inference (Supplementary Fig. S2). A notable observation is that “lineage I” is more closely related to lineage B, which was not detected beyond the 1970s<sup>20</sup> (mean *p*-distance of 0.17), than to more recent lineages like E and F which contain worldwide E30 strains including all E30 strains previously reported in Spain (mean *p*-distance between lineage I and lineages E and F of 0.25 and 0.26, respectively) (Supplementary Table S1).

All E11 strains obtained in this study belong to subgenotype D5 according to the classification introduced by Li et al.<sup>22</sup> which included stains detected in Europe and Asia that have circulated for over 20 years (Fig. 1C). E11 study strains were all highly similar (nt identity ≥ 95.3%; aa identity ≥ 97.6%) and were closely related to strains obtained in France in 2014 from meningitis cases in children (mean 96.5% (95.7–97.8%) nt and 98.6% (98.2–99.3%) aa similarities in the complete VP1 region).

Finally, Fig. 1D shows that E13 study strain LCR53 from 2016 cluster closely (bootstrap value of 100%) with two strains collected from children with meningitis in France in 2015 (nt identity ≥ 98.8%; aa identity ≥ 99.6%).

**Recombination analysis.** Considering the prevalence of recombination in the evolutionary history of enteroviruses, we next looked for evidence of recombination in the study strains. Recombination analyses were carried out using a combination of six methods implemented in RDP5.5<sup>23</sup>. We detected multiple events of recombination in the evolutionary history of the EV types studied. This was also evidenced when comparing the topologies of phylogenies built using the 3 coding regions (P1, P2 and P3) of the sequences from this study and representative EV-B reference genomes available in GenBank (Fig. 2). The phylogenetic analysis using the P1 capsid coding region showed that the E6, E11, E13 and E30 sequences obtained from this study cluster together with their respective reference genome which coincides with what we observed using the partial VP1 region (Fig. 2A). However, this was no longer observed when constructing the phylogeny based on the P2 and P3 non-capsid genomic regions (Fig. 2B,C). These incongruent tree topologies between the capsid and non-capsid regions indicate that genetic exchanges through intertypic recombination with other EV-B types might have occurred. We also used similarity plot and bootscanning analyses comparing the study strains and closely related types to visualize the mosaic genomes reported here. The sequences of closely related types included in
<figure><img src="image_3.png" /><figcaption>(A)</figcaption></figure>

■ Strains from Spain (2006-10, n=64) and France (2001-08, n=8)

■ Strains from Spain (2004-05, n=2) and France (2000-02, n=4)

■ Strains from Spain (2009, n=8) and Australia (2009, n=4)

○ Poland 2011

○ France 2011

<figure><img src="image_4.png" /><figcaption>(B)</figcaption></figure>

■ LCR675 and strains from Spain (1997-98, 2000, 2006 and 2016-18), Netherlands 2017 (n=1)

◇ Ireland 2017

■ Strains from Spain (2016-17, n=10)

■ Spain 2016-18

○ France 2016-17

¥ Austria 2016

▲ Netherlands 2017-18

<figure><img src="image_5.png" /><figcaption>(C)</figcaption></figure>

■ Strains from Spain (2002-08, n=28), Finland (2005-07, n=6), Netherlands (1997-2008, n=20), Russia (2004-09, n=11), Vietnam (2005, n=1), France (2001-08, n=8) and China (2016-19, n=4)

○ France 2012-14

△ Japan 2013-15

* Russia 2012

■ USA 2012

<figure><img src="image_6.png" /><figcaption>(D)</figcaption></figure>

■ Strains from China (2000-12, n=9)

■ Strains from India (2006-09, n=8)

■ Strains from Spain (2006, n=5) and France (2005-06, n=5)

■ Strains from Spain 2000 (n=64)

○ France 2011-15

**Figure 1.** Maximum-likelihood trees based on a partial VP1 inferred by IQ-TREE v.2.0.6 for the four genotypes study here: E6 (A), E30 (B), E11 (C) and E13 (D). Sequences reported in our study are indicated by black circles. Numbers on nodes indicate the bootstrap support of the node (> 80) and the scale bars represents the expected number of nucleotide substitutions per site. Non collapsed trees can be found in Supplementary Fig. 4A–D. Trees are annotated with the classification proposed by Smura et al. and Cabrerizo et al.<sup>10,19</sup> for E6, the classification proposed by Bailly et al.<sup>20</sup> for E30, and the classification proposed by Li et al.<sup>22</sup> for E11. Only bootstrap values above 80% are indicated in branch nodes. Scale bars indicate nucleotide substitutions per site.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Phylogenetic analysis based on the P1, P2 and P3 coding sequences of the 12 study strains and other fully sequenced EV-B global strains. Maximum-likelihood phylogenetic trees were inferred using IQ-TREE2 on an alignment of ~2565, 1734 and 2268 nt sequences corresponding to the P1 (A), P2 (B), and P3 (C) coding regions respectively. The numbers at the nodes indicate bootstrap support values > 80 for that node. Scale bar represents nucleotide substitutions per site. GenBank accession numbers for published sequences are shown in the tree.</figcaption></figure>

the analysis were those with highest score when using the genomic regions identified in the RDP5.5 analysis as of recombinant origin as queries for BLASTn (http://www.ncbi.nlm.nih.gov/) (Supplementary Table S2).

Genome sequences of E11 study strains were compared with all E6 study strains, the E11 prototype strain and sequences closely related to different parts of LCR484 which was identified as mosaic by all 6 methods within RDP5.5. While most of LCR484 genome presented high identity with the other E11 strains reported here (LCR1106, LCR1059), the 3' end of the genome (from the 3C region) presented high identity with non-E11 sequences. Genomes with similar recombinant genomic organization have been detected previously (Supplementary Fig. S3).

The LCR53 (E13) genome also presented strong signal of recombination, both intra and intertypic. In particular, we noted signal of recombination with other E13 strains in the 5' (P1 and P2 regions) of the genome. Most importantly, we detected intertypic recombination signal in the 3'end of the genome (from the 3C region) (Fig. 3). No similar recombinant genomes have been reported.

Among the E30 strains reported here, those from lineage I (LCR 519 and LCR520) and those from lineage F (LCR675 and LCR265) were detected as of recombinant origin, with signal in the 3'end of the genome (Supplementary Fig. S3). However, like for E11, genomes with similar mosaicism have been reported in the Netherland during the same year (Supplementary Fig. S3).

## Discussion

Application of mNGS methods to characterize the whole-genome sequence of viruses has become a valuable tool for tracking viral evolution, to elucidate the mechanisms shaping their diversity and to better understand the biology of viruses associated with human disease²⁴. In the present study, we used a previously described²⁵ protocol for mNGS of clinical samples on twelve CSF samples confirmed as EV-positive and collected from patients with aseptic meningitis during 2015–2018 in Southern Spain. It is important to note that the clinical samples were previously frozen and thawed at least once and stored at –45 °C. Our approach allowed to obtain complete or near-genomes despite the length of storage (samples were collected from 2015–2018 and processed in 2019), freeze–thaw cycles and the storage temperature above –480 °C. This suggests that this method can be effectively implemented in public health reference laboratories to detect and reconstruct complete viral genomes from CSF samples even if there has been a progressive degradation of RNAs. Moreover, unlike PCR amplicon-based or capture-based enrichment methods, this metagenomic approach does not require previous knowledge of the pathogen sequence, making this approach suitable to characterize novel viruses or recombinant forms. This is important as upsurges of emerging recombinant NPEVs associated with neurologic manifestations such
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Plot of similarity (A) and bootscanning analysis (B) of E13/Spain_LCR53/2016 study strain with prototype and closely related strains. The enterovirus genetic organization is shown in the top panel. Analyses were conducted by using SimPlot 3.5.1 (Kimura distance model, window size 400 bp moving in 40 nucleotides steps).</figcaption></figure>

as EVA71 or E30 are a growing public health concern in Europe²⁶,²⁷. A real-time detection and rapid sharing of whole-genome sequence information on these emerging NPEVs from public health reference laboratories may help detect outbreaks of severe neurologic disease and implement early intervention strategies ensuring rapid control of disease spread.

The types identified in this study were E6, E11, E13 and E30. Echoviruses belong to the EV-B species and constitute the largest group of the EV genus, with 28 types. A meta-analysis of the worldwide distribution of EVs showed that E30, E6 and E13 were the most commonly detected human EV types in meningitis cases²⁸. In Europe, between 2015 and 2017, E30, E6 and E11 were among the ten most frequently reported types representing 12%, 12% and 4%, respectively, of all typed EV-positive samples⁵. In 2018, an E30 upsurge was observed in Denmark, Germany, The Netherlands, Norway and Sweden associated with meningitis or meningoencephalitis²⁶. A recent study identified that this upsurge was caused by co-circulation of E30 viruses from two different lineages: E and I (corresponding to clades G1 and G6, respectively) that replaced E30 viruses from lineage F (corresponding to clade G2) which predominated in 2016 and 2017²¹. We similarly note that the E30 strains detected in Spain in 2017 are from lineage F, and that the subsequent strains (from 2018) correspond to the emerging lineage I. The close genetic relatedness of the E30 sequences reported in lineage I during 2017–2018 in the Netherlands, France, Ireland and Spain is further indicative of rapid widespread transmission of this lineage²⁹.

In Spain, from 1988 to 2008, E30 has been the predominant EV type (33.7% of the total typed EVs)⁸,³⁰. However, in 2009, an upsurge of E6 (60%) replaced E30 as the predominant type. In more recent years (2016–2019), E30 and E6 were among the ten most prevalent types of all typed EVs in Spain⁶,³¹. All E6 study strains analyzed in this study were from 2015 and segregated into the same sublineage C1, together with previously described
strains from France and Poland isolated in 2011, indicating that related C1 strains were circulating during those years over a broad geographic area in Europe. Interestingly, E6 strains circulating in Spain (n = 86) in the previous decade (2000–2010), all clustered into sublineage C9 suggesting again, as was observed for E30, a genetic turnover where C9 viruses might have subsequently be replaced by C1 viruses.

In this study, phylogenetic analysis showed that E11 study strains from 2015 and from 2018 belong to the same genomic cluster (subgenotype D5) which is distributed worldwide. Study strains have close relationship with each other (> 95% identity) suggesting sustained local circulation in the region during the study period. The fact that strains detected in Spain from 2002 to 2008 also cluster in D5 supports the idea that this subgenotype has persisted for almost two decades in the country and is in agreement with previous studies describing relative genetic stability of E11 strains, associated with long-term endemic patterns¹³,³². Of note, although E11 study strains grouped within the previously described D5 cluster, they formed a separate subcluster suggesting genetic evolution of this endemic E11 viral population from a common ancestor.

Consistently with Europe where E13 is among the less reported types (<1%), E13 is also rarely reported in Spain⁵,⁶,³³. It appears to circulate in a cyclical fashion with epidemic years followed by years with low detection rates³³. In Spain, E13 reporting increased in 2000 and 2019 (45% and 8% of the total typed EVs, respectively) with almost no detection in between⁶,³³. In this study, we detected a sporadic case of aseptic meningitis associated with an E13 strain from 2016 closely related to E13 strains from France also collected in 2016 from children with meningitis. This close genetic relatedness between E13 sequences reported in 2016 in both France and Spain suggests rapid widespread transmission.

Recombination is known to constitute a widespread mechanism of EV evolution, through the continual exchange of fragments between co-circulating viruses which may impact viral replication and pathogenicity. Full genome sequencing is crucial to detect these recombination events. Our analysis revealed that the E6, E11, E13 and E30 circulating types present complex mosaic genomes constituted of fragments of different origins, implicating intertypic recombination events.

Our analysis highlights the continued circulation of recombinant E11 viruses carrying non-structural sequences similar to other (non-E11) types. Evidence of evolution of E6 and E11 through intertypic recombination has been suggested before³⁴, and the detection of strains with similar mosaic organization in multiple parts of the world (USA, Japan, Italy) indicates that this lineage has been circulating worldwide in recent years, in parallel to other E11 lineages.

Within E30, the lineage F strains reported here (LCR265, 2016 and LCR675, 2017) present a mosaic pattern previously detected in genomes from Spain (Catalonia) in 2016, as well as genomes noted in the Netherlands in 2017. In a similar way, the lineage I study strains (LCR519 and LCR520, 2018) also present mosaic patterns previously detected in genomes from the Netherlands and Ireland in 2017 and from Austria in 2016. These results suggest that these recombinant lineages have been circulating for some time and started to diversify. Multiple reports highlight the involvement of genetic recombination in the evolutionary process of E30 with other EV-B strains in the non-structural region¹³,²⁹,³⁵.

The E13 study strain (LCR53) presents a mosaic organization with no reported homologs, suggesting that the recombination event leading to the circulation of this virus might be recent, or that it has been circulating at low noise and was not captured up to now. The detected recombination breakpoint lies at the beginning of the P3 region, providing further support to the idea that the structural and non-structural genome regions of EV-B types evolve independantly⁶. To better monitor recombination events and more accurately identify potential recombination partners or location of emergence (including potential geographical hotspots), stronger genomic surveillance programs for EVs are needed.

In Europe, EVs are more commonly detected in late summer and autumn⁵. Our study revealed that most types (E6, E13 and E30) were detected during the spring season (March-May) which is consistent with previous reports describing the incidence of EV-B in Spain peaking during spring⁸,³⁰,³⁷. Rather unexpectedly, all three cases associated to E11 occurred in winter. Different seasonality in circulating EV types and the fact that seasonality in temperate zones such as Spain is less marked than in tropical climates could explain these differences.

One limitation of this study is the limited number of samples here analyzed within a 3-year period (only severe cases requiring hospitalization). For this reason, we are likely underestimating the genetic diversity of EVs in Southern Spain. This emphasizes the need to improve the genomic surveillance of EVs in Spain by screening larger cohorts of patients and obtaining detailed epidemiological and clinical data. This could be complemented with increased nationwide surveillance and further exploration of environmental samples using mNGS to monitor emergent and/or novel recombinant types and to better understand their circulation patterns and potential recombination events among co-circulating EVs.

In conclusion, this study documents for the first time the whole-genome sequences of E6, E11 and E13, as well as E30 circulating in Spain. The detection of a novel, disease-associated, and recombinant form highlights the importance of genomic surveillance characterizing full-length EV genomes and public sharing of whole-genome sequence data.

## Materials and methods

**Ethics statement.** Ethical approval for the experimental protocol was given by the Research Ethics Committee from Cordoba (Comité de Ética de Investigación de Córdoba, reference 201999903552445). The procedures were carried out in accordance with approved guidelines, regulations and the principles of the Declaration of Helsinki. Written informed consent was obtained from all patients or from a parent and/or legal guardian.

**Study background.** From January 2015 to December 2018, eleven adult patients were diagnosed in the Department of Neurology in the University Hospital Reina Sofia (Córdoba, Spain) with a laboratory-confirmed
enteroviral meningitis. The Department of Neurology is responsible for the treatment of patients from 14 years old and is the reference unit for patients with suspected CNS infections in the province. The University Hospital Reina Sofia is a 1280-bed tertiary referral hospital for the southwestern Spanish province of Córdoba (population 461.078). A case was defined as meningitis-suspected based upon the presence of these symptoms: fever, headache, vomiting, neck stiffness, nausea and sometimes accompanied by photophobia, abdominal pain or diarrhea. CSF samples were obtained at patient admission according to standard procedures and analyzed at the Department of Microbiology. Samples were tested for the presence of EVs using a qualitative multiplexed PCR (FilmArray® Meningitis/Encephalitis Panel, BioFire Diagnostics) or RT-PCR (Xpert EV, Cepheid). In all CSF samples, bacterial, herpes simplex virus types 1 and 2, and varicella-zoster virus infections of the CNS were excluded by culture or PCR. Bacterial screening was performed by culture on chocolate agar, blood agar and thioglycollate broth, while herpes simplex virus and varicella zoster virus screening was performed by a qualitative multiplexed PCR (FilmArray® Meningitis/Encephalitis Panel, BioFire Diagnostics) or RT-PCR (RealCycler® herpesvirus type 1 + herpesvirus type 2 + varicella-zoster virus, Progenie Molecular). Sample remnants were stored at -45 °C. Clinical, epidemiological and laboratory data were collected retrospectively from medical records.

**Library preparation and metagenomic sequencing.** RNA was extracted from 140 µl of CSF using the QIAamp Viral RNA Mini Kit (Qiagen) followed by Turbo DNAse treatment (Ambion) and purification with Agencourt RNAClean XP beads. We next depleted host rRNA using custom probes and RNAs e H treatment as previously described²⁵. Depleted samples were purified using AMPure RNA clean beads (Beckman Coulter Genomics) and eluted in 10 ul of RNase-free water for cDNA synthesis. RNA was converted to double-stranded cDNA in two steps. First, RNA was retro transcribed using random primers and SuperScript IV (Invitrogen). Second-strand cDNA was generated using E. coli DNA ligase, RNAs e H and DNA polymerase (New England Biolabs) and purified using Agencourt AMPure XP beads (Beckman Coulter). Libraries were then prepared using the Nextera XT DNA Library Prep Kit (Illumina) and sequenced on an Illumina NextSeq500 (2×75 cycles).

**Genome assembly and analysis of sequence data.** Raw reads were trimmed using Trimmomatic v0.36 to remove low-quality reads and Illumina adaptors. Reads were de novo assembled using the metaspades option from SPAdes v3.12 and the contigs obtained were used as blast queries on Virus Pathogen Resource (ViPR)³⁸. We used clc-assembly-cell v5.1.0 to perform mapping and extract consensus. A minimum of 3X read depth coverage was used, and N were added in case of lower coverage. Samtools v1.3 was used to sort the aligned BAM files and to generate alignment statistics. All alignments and consensus sequences were manually inspected using Geneious Prime 2020.2 (https://www.geneious.com/).

**Recombination analysis.** Recombination analyses were carried out on the study strain using a background of randomly sampled EV genomes retrieved from GenBank on May 2021, using a combination of six methods implemented in RDP5.5²³ (RDP, GENECONV, MaxChi, Bootscan, SisScan and 3SEQ) and we considered recombination signals detected by more than three methods for breakpoint identification. Except for specifying that sequences are linear, all settings were kept to their defaults. In addition, similarity plot and bootscanning analysis were performed by using the SimPlot program, version 3.5.1, with a 400-nt window moving in 40-nt steps and using a Kimura 2-parameter method with a transitions-transversions ratio of 2 with 1000 resampling.

**Phylogenetic analysis.** We performed phylogenetic analysis of the 12 genomes from our dataset, together with published genomes available on NCBI GenBank retrieved on May 2021. Because the majority of the reported sequences from strains circulating in Spain consisted of partial VP1 sequences, we first performed phylogenetic analysis based on the alignment of partial 3'-half VP1 region (E6, E30, and E13) or in the 5'-half VP1 region (E11). Nt sequence alignment was performed by using ClustalW multiple alignment program³⁹ within the BioEdit sequence Alignment Editor package, version 7.0.9.0. Maximum-likelihood trees were constructed using IQ-TREE2 (v.2.0.6)⁴⁰ with 1000 bootstrap replicates⁴¹. We constructed the trees using the best-fit model as determined by ModelFinder⁴² implemented in IQ-TREE2. Genetic distances and sequence divergence were calculated using the Molecular Evolutionary Genetics Analysis (MEGA-X) software package (http://megasoftwa.re.net/).

**Data availability**

Raw sequencing reads were deposited on the European Nucleotide Archive (https://www.ebi.ac.uk/ena/brows er/home) with study number PRJEB45068. The GenBank (https://www.ncbi.nlm.nih.gov/genbank/) accession numbers of the assembled virus genomes are MZ389224–MZ389234 and MZ436966. Both ENA and GenBank are part of the International Nucleotide Sequence Database Collaboration (INSDC https://www.insdc.org/).

Received: 16 July 2021; Accepted: 20 October 2021

Published online: 02 November 2021

**References**

1. Palacios, G. & Oberste, M. S. Enteroviruses as agents of emerging infectious diseases. *J. Neurovirol.* **11**, 424–433 (2005).

2. The Online (10th) Report of the International Committee on Taxonomy of Viruses (2017). Available at: https://talk.ictvonline.org/ictv-reports/ictv_online_report/positive-sense-rna-viruses/picornavirales/w/picornaviridae/681/genus-enterovirus.

3. Oberste, M. S., Maher, K. & Pallansch, M. A. Evidence for frequent recombination within species human enterovirus B based on complete genomic sequences of all thirty-seven serotypes. *J. Virol.* **78**, 855–867 (2004).
4. Harvala, H. et al. Recommendations for enterovirus diagnostics and characterisation within and beyond Europe. *J. Clin. Virol.* https://doi.org/10.1016/j.jcv.2018.01.008 (2018).

5. Bubba, L. et al. Circulation of non-polio enteroviruses in 24 EU and EEA countries between 2015 and 2017: A retrospective surveillance study. *Lancet. Infect. Dis.* https://doi.org/10.1016/S1473-3099(19)30566-3 (2020).

6. Annual Epidemiological reports of Acute Flaccid Paralysis Surveillance and Enterovirus Surveillance, Spain, 2016–2019. Available at: https://www.isciii.es/QueHacemos/Servicios/VigilanciaSaludPublicaRENAVE/EnfermedadesTransmisibles/Paginas/Resultados_Vigilancia_Polio.aspx.

7. Cabrerizo, M., Echevarria, J. E., González, I., De Miguel, T. & Trallero, G. Molecular epidemiological study of HEV-B enteroviruses involved in the increase in meningitis cases occurred in Spain during 2006. *J. Med. Virol.* https://doi.org/10.1002/jmv.21197 (2008).

8. Trallero, G. et al. Enteroviruses in Spain over the decade 1998–2007: Virological and epidemiological studies. *J. Clin. Virol.* https://doi.org/10.1016/j.jcv.2009.11.013 (2010).

9. Andrés, C. et al. Recombinant CV-A6 strains related to hand-foot-mouth disease and herpangina at primary care centers (Barcelona, Spain). *Fut. Microbiol.* https://doi.org/10.2217/fmb-2018-0336 (2019).

10. Cabrerizo, M., Trallero, G. & Simmonds, P. Recombination and evolutionary dynamics of human echovirus 6. *J. Med. Virol.* https://doi.org/10.1002/jmv.23741 (2014).

11. Calderón, K. I. et al. Molecular epidemiology of coxsackievirus B3 infection in Spain, 2004–2014. *Adv. Virol.* https://doi.org/10.1007/s00705-016-2783-1 (2016).

12. Leitch, E. C. M. et al. Transmission networks and population turnover of echovirus 30. *J. Virol.* https://doi.org/10.1128/jvi.02109-08 (2009).

13. McWilliam Leitch, E. C. et al. Evolutionary dynamics and temporal/geographical correlates of recombination in the human enterovirus echovirus types 9, 11, and 30. *J. Virol.* https://doi.org/10.1128/jvi.00783-10 (2010).

14. McWilliam Leitch, E. C. et al. The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71. *J. Virol.* 86, 2676–2685 (2012).

15. González-Sanz, R. et al. Molecular epidemiology of an enterovirus A71 outbreak associated with severe neurological disease, Spain, 2016. *Eurosurveillance* https://doi.org/10.2807/1560-7917.ES.2019.24.7.1800089 (2019).

16. Leon, K. E. et al. Genomic and serologic characterization of enterovirus A71 brainstem encephalitis. *Neurol. Neuroimmunol. Neuroinflam.* https://doi.org/10.1212/NXI.0000000000000703 (2020).

17. Puenpa, J. et al. Molecular epidemiology and the evolution of human coxsackievirus A6. *J. Gen. Virol.* https://doi.org/10.1099/jgv.0.000619 (2016).

18. de Lima, S. T. S. et al. Fatal outcome of chikungunya virus infection in Brazil. *Clin. Infect. Dis.* https://doi.org/10.1093/cid/ciaa1038 (2020).

19. Smura, T. et al. Molecular evolution and epidemiology of echovirus 6 in Finland. *Infect. Genet. Evol.* 16, 234–247 (2013).

20. Bailly, J. L. et al. Phylogeography of circulating populations of human echovirus 30 over 50 years: Nucleotide polymorphism and signature of purifying selection in the VP1 capsid protein gene. *Infect. Genet. Evol.* https://doi.org/10.1016/j.meegid.2008.04.009 (2009).

21. Benschop, K. et al. Molecular epidemiology and evolutionary trajectory of emerging echovirus 30. *Eur. Emerg. Infect. Dis.* 27(6), 1616–1626. https://doi.org/10.3201/eid2706.203096 (2021).

22. Li, J. et al. Multiple genotypes of Echovirus 11 circulated in mainland China between 1994 and 2017. *Sci. Rep.* https://doi.org/10.1038/s41598-019-46870-w (2019).

23. Martin, D. P. et al. RDP5: A computer program for analyzing recombination in, and removing signals of recombination from, nucleotide sequence datasets. *Virus Evol.* https://doi.org/10.1093/ve/veea087 (2021).

24. Houldcroft, C. J., Beale, M. A. & Breuer, J. Clinical and biological insights from viral genome sequencing. *Nat. Rev. Microbiol.* https://doi.org/10.1038/nrmicro.2016.182 (2017).

25. Matranga, C. B. et al. Enhanced methods for unbiased deep sequencing of Lassa and Ebola RNA viruses from clinical and biological samples. *Genome Biol.* 15, 519 (2014).

26. Broberg, E. K. et al. Upsurge in echovirus 30 detections in five EU/EEA countries, April to September, 2018. *Eurosurveillance* https://doi.org/10.2807/1560-7917.ES.2018.23.44.1800537 (2018).

27. Ngangas, S. T. et al. Multirecombinant enterovirus A71 subgenogroup C1 isolates associated with neurologic disease, France, 2016–2017. *Emerg. Infect. Dis.* https://doi.org/10.3201/eid2506.181460 (2019).

28. Suresh, S., Rawlinson, W. D., Andrews, P. I. & Stelzer-Braid, S. Global epidemiology of nonpolio enteroviruses causing severe neurological complications: A systematic review and meta-analysis. *Rev. Med. Virol.* https://doi.org/10.1002/rmv.2082 (2020).

29. Mirand, A., Henquell, C., Archimbaud, C., Piegue-Lafeilleu, H. & Bailly, J. L. Emergence of recent echovirus 30 lineages is marked by serial genetic recombination events. *J. Gen. Virol.* https://doi.org/10.1099/vir.0.82146-0 (2007).

30. Cabrerizo, M., Echevarria, J. E., Gonzalez, I., de Miguel, T. & Trallero, G. Molecular epidemiological study of HEV-B enteroviruses involved in the increase in meningitis cases occurred in Spain during 2006. *J. Med. Virol.* 80, 1018–1024 (2008).

31. Cabrerizo, M. et al. Molecular epidemiology of enterovirus and parechovirus infections according to patient age over a 4-year period in Spain. *J. Med. Virol.* https://doi.org/10.1002/jmv.24658 (2017).

32. Oberste, M. S., Nix, W. A., Kilpatrick, D. R., Flemister, M. R. & Pallansch, M. A. Molecular epidemiology and type-specific detection of echovirus 11 isolates from the Americas, Europe, Africa, Australia, southern Asia and the Middle East. *Virus Res.* https://doi.org/10.1016/S0168-1702(02)00291-5 (2003).

33. Trallero, G. et al. First epidemic of aseptic meningitis due to echovirus type 13 among Spanish children. *Epidemiol. Infect.* https://doi.org/10.1017/S0950268802008191 (2003).

34. Su, T. et al. Molecular characterization of a new human echovirus 11 isolate associated with severe hand, foot and mouth disease in Yunnan, China, in 2010. *Adv. Virol.* https://doi.org/10.1007/s00705-015-2496-x (2015).

35. Lukashev, A. N., Ivanova, O. E., Eremeeva, T. P. & Gmyl, L. V. Analysis of echovirus 30 isolates from Russia and new independent states revealing frequent recombination and reemergence of ancient lineages. *J. Clin. Microbiol.* https://doi.org/10.1128/JCM.02386-06 (2008).

36. Lukashev, A. N. et al. Recombination in circulating Human enterovirus B: Independent evolution of structural and non-structural genome regions. *J. Gen. Virol.* https://doi.org/10.1099/vir.0.81264-0 (2005).

37. Rodà, D. et al. Clinical characteristics and molecular epidemiology of Enterovirus infection in infants <3 months in a referral paediatric hospital of Barcelona. *Eur. J. Pediatr.* https://doi.org/10.1007/s00431-015-2571-z (2015).

38. Pickett, B. E. et al. ViPR: An open bioinformatics database and analysis resource for virology research. *Nucleic Acids Res.* 40, 2 (2012).

39. Thompson, J. D., Higgins, D. G. & Gibson, T. J. CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. *Nucleic Acids Res.* https://doi.org/10.1093/nar/22.22.4673 (1994).

40. Minh, B. Q. et al. IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. *Mol. Biol. Evol.* https://doi.org/10.1093/molbev/msaa015 (2020).

41. Hoang, D. T., Chernomor, O., Von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: Improving the ultrafast bootstrap approximation. *Mol. Biol. Evol.* https://doi.org/10.1093/molbev/msx281 (2018).
42. Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., Von Haeseler, A. & Jermiin, L. S. ModelFinder: Fast model selection for accurate phylogenetic estimates. *Nat. Methods* https://doi.org/10.1038/nmeth.4285 (2017).

## Acknowledgements

This work was funded by the Junta de Andalucía, Consejería de Salud y Familias (Project number PI-0216-2019) and by the Acción Estratégica en Salud Intramural (Project number PI20CIII/00005). MD Fernandez-Garcia received a Miguel Servet Research Contract (CP18/00067) from the Strategic Action in Health 2018 and funded by National Institute of Health Carlos III (ISCIII). ESL acknowledges funding from the INCEPTION programme (Investissements d'Avenir grant ANR-16-CONV-0005).

## Author contributions

Conceptualization, M.D.F.G.; Coordinated and supervised the study, E.S.L., M.D.F.G.; Funding acquisition, E.S.L., M.D.F.G.; Collection of clinical samples and clinical information, E.A.; Performed experiments, F.G., A.B.P., M.P., M.D.F.G.; Analyzed and interpreted the data, E.S.L., F.G., M.D.F.G.; Writing—Original draft preparation, E.S.L., M.D.F.G., F.G.; Writing—Review and editing, L.M.M., M.C., A.B.P., E.A.; All authors have read and agreed to the published version of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-01053-4.

**Correspondence** and requests for materials should be addressed to E.S.-L. or M.D.F.-G.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021, corrected publication 2021