1. Background {#sec0005}
=============

Human enterovirus D68 (EV-D68) was originally isolated in 1962, and for many years it was considered a rare cause of respiratory disease because it was identified in isolated cases and in small outbreaks \[[@bib0005]\]. Although outbreaks of EV-D68 infection occurred in the Philippines \[[@bib0010]\], Italy \[[@bib0015], [@bib0020]\], the Netherlands \[[@bib0025]\] and Japan \[[@bib0030]\] between 2005 and 2012, it was only in 2014 that the virus gained epidemiological and clinical relevance. In that year, a large-scale outbreak of severe respiratory infections due to this agent, in some cases associated with central nervous system diseases, occurred in the USA and Canada \[[@bib0035]\], with subsequent cases in other countries \[[@bib0040], [@bib0045], [@bib0050], [@bib0055], [@bib0060], [@bib0065], [@bib0070], [@bib0075], [@bib0080]\]. Since that period, several outbreaks in different continents have been reported, and it was confirmed that the virus can cause upper and lower respiratory tract diseases and be associated with the development of neurological problems, including acute flaccid paralysis (AFP), encephalitis and aseptic meningitis \[[@bib0085], [@bib0090], [@bib0095]\].

Phylogenetic analysis of EV-D68 strains detected worldwide has led to the identification of four distinct viral clades that are designated groups A--D according to the characteristics of the structural protein VP1, the most variable region of the EV genome \[[@bib0100]\]. Moreover, within clades, subgroups were identified according to amino acid substitutions in structural and nonstructural proteins. Several studies have shown that these different clades and subgroups can circulate or co-circulate during different periods and be associated with both mild and severe diseases \[[@bib0105]\]. However, little is known about the molecular evolution of this virus during outbreaks, and it is not precisely defined whether specific genetic variances are responsible for the most severe cases. Studies regarding the molecular phylogeny, diversity, and epidemiology of EV-D68 are essential to understanding the evolution of this virus and its adaptation to the human host. Moreover, they can contribute to the development of effective preventive and therapeutic measures.

2. Objectives {#sec0010}
=============

The main aim of this research was to study the genetic characteristics of EV-D68 strains that were circulating in Italy during the period March-October 2016.

3. Study design {#sec0015}
===============

3.1. Patients and sample {#sec0020}
------------------------

A retrospective study of the circulation of different types and subtypes of EV-D68 was performed. Nasopharyngeal swabs collected from March 2016 through October 2016 in children admitted to the Emergency Room with respiratory diseases and stored in a freezer at −80 °C were tested. The sample collection had occurred in the Pediatric Department of the Fondazione IRCCS Policlinico San Matteo, University of Pavia (Italy) and in the Pediatric Highly Intensive Care Unit of the Fondazione IRCCS Ca' Granda Policlinico, University of Milan (Italy). The Ethics Committees of both Institutions approved the study. Written informed consent was obtained from a parent of a legal guardian of each enrolled child. Children aged ≥8 years were asked for their assent. Clinical data were retrieved from archived clinical charts.

3.2. EV-D68 detection {#sec0025}
---------------------

Respiratory samples were extracted using a QIAsymphony^®^ instrument with a QIAsymphony^®^ DSP Virus/Pathogen Midi Kit (Complex 400 protocol) according to the manufacturer's instructions (QIAGEN, Qiagen, Hilden, Germany). Panels of laboratory-developed real-time RT-PCR or real-time PCR \[[@bib0110], [@bib0115]\] were used to detect and quantify the following viruses: influenza virus A and B (including subtype determination), human Rhinovirus (hRV), human parainfluenza virus (hPIV) 3, respiratory syncytial virus (RSV) types A and B, human coronavirus (hCoV)-OC43, -229E, -NL63, and -HKU1, human metapneumovirus (hMPV) and human adenoviruses (hAdVs). Extracted samples were also amplified using EV-D68-specific real-time RT-PCR as previously described \[[@bib0120]\].

3.3. VP1 sequencing and phylogenetic analysis {#sec0030}
---------------------------------------------

Complete VP1 sequences were amplified for all EV-D68-positive samples according to the method of Piralla et al. \[[@bib0020]\]. Purified PCR products were sequenced using a BigDye Terminator Cycle-Sequencing Kit (Applied Biosystems, Foster City, USA) in an ABI Prism 3130xl Genetic Analyzer (Applied Biosystems, Foster City, USA).

The sequences were assembled using Sequencher software, version 4.6 (Gene Codes Corporation, Ann Arbor, USA). Sequence alignment and phylogenetic analysis were performed using the MEGA software, version 7 with the TN92 + G nucleotide substitution model, which was selected using the hierarchical likelihood ratio test \[[@bib0125]\]. Statistical support of specific clades in the inferred tree was evaluated by using the bootstrap analysis (1000 replicates).The VP1 secondary structure was investigated using the Jpred4 server \[[@bib0130]\].

The nucleotide dataset of the VP1 gene was obtained by including sequences from complete genome (n = 84), partial VP1 (n = 20) from the Netherlands outbreak 2016 \[[@bib0135]\] as well as from EV-D68 strain sequences that were circulating in Italy during the period 2008--2014 (n = 30). The EV-D68 sequences originated in this study were deposited in the GenBank database (accession numbers: MF073335--MF073354).

3.4. Selective pressure analysis {#sec0035}
--------------------------------

Tests for positive selection were conducted using single-likelihood ancestor counting (SLAC), fixed-effects likelihood (FEL), random-effects likelihood (REL), internal branch fixed-effects likelihood (IFEL), the mixed-effects model of evolution (MEME), and fast unconstrained Bayesian approximation (FUBAR) on the Datamonkey server \[[@bib0140]\], with the dN/dS ratios being calculated using the SLAC and FEL codon-based maximum likelihood approaches. SLAC counts the number of non-synonymous changes per non-synonymous site (dN) and tests whether it is significantly different from the number of synonymous changes per synonymous site (dS). FEL estimates the ratios of non-synonymous to synonymous changes for each site in an alignment \[[@bib0145]\]. The IFEL method is similar to FEL but tests for site-by-site selection only along internal branches of the phylogeny. To avoid an excessive false-positive rate, sites with SLAC, FEL, IFEL and MEME p-values of \<0.1 and a FUBAR posterior probability of \>0.90 were accepted as candidates for selection.

4. Results {#sec0040}
==========

4.1. Clinical and epidemiological data {#sec0045}
--------------------------------------

Samples obtained from 390 children, 270 from Pavia and 120 from Milan, were tested. Twenty-two patients (13 males; 59.1%), mean age 47 months (range, 1--174 months) were found to be infected by EV-D68 ([Table 1](#tbl0005){ref-type="table"} ). The majority of cases (14/22; 63.6%) occurred in June (n = 9) and July (n = 5). The remaining cases were observed in April (n = 1), May (n = 2), September (n = 2), and October (n = 3). Pneumonia was diagnosed in 12 (54.5%) children, whereas viral wheezing, rhinopharyngitis, asthma and bronchiolitis, were diagnosed in four children, three children, two children, and one child, respectively. The previous clinical history was negative in most of the cases (16/22; 72.7%). However, two (9.1%) children previously suffered from leukemia and had been recently transplanted with hematopoietic stem cells; 2 (9.1%) children had a history of recurrent episodes of viral wheezing; and 2 (9.1%) children had allergic asthma. The disease severity was judged to be mild-to-moderate in all the patients. Lower than 92% oxygen saturation (SaO~2~) values at the first visit in the emergency room were found only in 5 (22.7%) children. Hospitalization occurred in 19/22 (86.4%) children, with a mean duration of 4 days (range, 4--8 days). None of the patients had to be admitted to the intensive care unit. In 16/22 (72.7%) cases, EV-D68 was the only infectious agent detected in the respiratory samples. Coinfection with human adenovirus (hAdV) and with human rhinovirus (hRV) was detected in 4 (18.2%) and 2 (9.1%) of the children, respectively. In all these cases, coinfecting agents were detected in very low concentrations (almost all below 1.0 × 10^3^ copies/mL). The median EV-D68 load observed was 1.3 × 10^5^ RNA copies/mL of respiratory sample (range, 9.0 × 10^2^--5.1 × 10^6^ RNA copies/mL of respiratory sample). No association was found between disease severity (measured by the SaO~2~ values) and EV-D68 load. The 368 samples negative for EV-D68 resulted positive for hRV (60.2%), RSV (26.6%), bocavirus (10.1%), hMPV (8.2%), and hAdV (7.9%). In 48 (13%) cases, particularly when hRV was detected, coinfection between two viruses was evidenced. Ten children with RSV infection, two with hMPV and one with hAdV were hospitalized for severe disease and three of those with RSV required admission to the intensive care unit.Table 1Epidemiologic and clinical data of EV-D68 positive patients diagnosed in Italy during the period March 1- October 30, 2016.Table 1Patient \#Sampling dateSex/Age (months)DiagnosisSaO~2~ (%)Viral coinfection (VL in copies/mL)Underlying diseaseDays of hospitalizationEv-d68 VL1Apr 18F/8Bronchiolitis93NoneNone58.7 × 10^5^2May 20M/174Rhinopharyngitis99NoneALL; previous HSCTNone4.0 × 10^3^3May 23M/50Rhinopharyngitis99NoneALL, previous HSTCNone9.0 × 10^2^4Jun 1M/33Pneumonia96HRV (7.0 × 10^5^)None45.1 × 10^6^5Jun 3F/60Asthma88NoneSensitivity to *Ambrosia* and *Dermatophagoides pteronyssinus*57.2 × 10^4^6Jun 5F/1Rhinopharyngitis99NoneNone44.8 × 10^5^7Jun 6M/66Pneumonia96NoneNone43.1 × 10^4^8Jun 20F/41Pneumonia94hAdV (NA)None41.9 × 10^6^9Jun 27M/22Pneumonia94NoneNone32.7 × 10^6^10Jun 28F/59Viral wheezing96hAdV (90)None41.8 × 10^5^11Jun 29M/26Pneumonia79NoneNone82.5 × 10^4^12Jun 30M/20Pneumonia92NoneNone44.5 × 10^3^13Jul 1F/74Viral wheezing88hAdV (135)Recurrent wheezing54.6 × 10^4^14Jul 11F/36Viral wheezing84hAdV (180)Recurrent wheezing52.2 × 10^3^15Jul 13F/53Pneumonia94NoneNone38.1 × 10^4^16Jul 18M/22Pneumonia95NoneNone33.7 × 10^5^17Jul 22F/35Pneumonia95NoneNone31.5 × 10^4^18Sep 20M/68Pneumonia96NoneNone3NA19Sep 23M/35Viral wheezing98NoneNoneNone2.4 × 10^6^20Oct 3M/10Pneumonia94NoneNone39.2 × 10^5^21Oct 10M/106Asthma85HRV(NA)Sensitivity to *Dermatophagoides pteronyssinus*64.4 × 10^5^22Oct 30M/35Pneumonia95NoneNone3NA[^1]

4.2. Phylogenetic and sequence analyses {#sec0050}
---------------------------------------

Complete VP1 sequencing was obtained for all the EV-D68-positive patients. Phylogenetic analysis of the VP1 region showed that all the strains identified in this study belonged to clade B3, while EV-D68 strains circulating in Italy in the period 2008--2014, highlighted in red, clustered in the others clades and subclades ([Fig. 1](#fig0005){ref-type="fig"} ). Within B3 subclade, the Italian EV-D68 strains were most closely related to strains detected in US and the Netherlands in 2016 and Southern China and Japan in 2014--2015.Fig. 1Maximum likelihood phylogenetic tree of EV-D68 strains. The evolutionary history was inferred by using the Maximum Likelihood method based on the Tamura 3-parameter model. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories +G, parameter = 0,5031). Significant bootstrap values (\>70) are indicated on the branch. EV-D68 strain sequences originally found in this study are indicated by black circle (n = 20); sequences from the Netherlands (n = 20) and USA 2016 outbreaks are indicated by blue and red circles, respectively. The other EV-D68 Italian sequences circulated in 2008--2014 period are colored in red (n = 30). The patient numbers presented in [Table 1](#tbl0005){ref-type="table"} are reported in brackets after the strain name. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 1

The average nucleotide identity between Italian EV-D68 subclade B3 strains was 98.9% (range, 96.9--100%), with a maximum of 21 nucleotides and 3 amino acid changes Eight β-sheets (βB to βI) were predicted in the VP1 sequences by the Jpred4 server. Among B3 subclade were also clustered EV-D68 strains identified in Japan and China between 2014 and 2015. In addition, the average nucleotide identity between EV-D68 subclade B3 strains and B1 and B2 subclade strains was 92.3% and 94.7%, respectively. Inside the motifs corresponding to the BC and DE loops, no changes were observed compared with other EV-D68 strains belonging to clade B1-3 ([Fig. 2](#fig0010){ref-type="fig"} ). The EV-D68 strains belonging to the new circulated B3 clade are characterized by 2 unique substitutions (A6G, S218T; [Fig. 2](#fig0010){ref-type="fig"}). Codon 218 is included in the GH loop that is also involved in the sialic acid receptor binding site \[[@bib0150]\]. In addition, four strains harbored the I258 V mutation (highlighted with an asterisk in [Fig. 1](#fig0005){ref-type="fig"}) that was only present in the EV-D68 ATCC strains (EVD68/Homo sapiens/XXX/ATCC VR-1197; KT725431). One strain (ITA/MI-ESC128/2016) was characterized by the K256I change, which had never been previously observed. Of note, one strain (ITA/PV-38993/2016) was characterized by the E168K change. This mutation was also present in the Fermon reference strain (AY4263531) but was more present in a few EV-D68 sequences belonging to all the three EV-D68 clades.Fig. 2Amino acid changes in VP1 sequences according to different clades. Amino acid residues characterizing clade B3 EV-D68 strains are indicated with black square. Amino acids included in VP1 motifs corresponding to protein loops are highlighted with a dashed line. The gaps are indicated by a dash (−) and the conserved amino acid residues by a dot (.). Numbering of amino acidic changes is according to strain Enterovirus D68 strain USA/ID/2014-19068 (KT995555) belong to B1 clade.Fig. 2

To assess the signatures of evolution in the VP1 gene region, the overall dN/dS ratio was estimated as 0.08. Several codons (more than 150) were identified to be under negative selection by at least three of the methods used (SLAC, FEL, REL, FUBAR or MEME). Similarly, no positively selected codons were identified by the IFEL model used to determine the selection pressure acting on the codons along the internal branches of the tree.

5. Discussion {#sec0055}
=============

This retrospective study did not test all the nasopharyngeal samples collected in children with respiratory symptoms during the period March-October 2016 in Pavia and Milan. Consequently, no definitive data regarding the epidemiology and clinical relevance of this infectious agent in children could be drawn. However, the collected data enables the integration of available information regarding the circulation of and genetic differences among different EV-D68 clades and can help elucidate the molecular evolution of this virus. The four EV-D68 clades have exhibited a differing predominance in recent years. The Italian EV-D68 strains identified in 2008 belonged to clade A and clade C. Clade A was again identified in 2010, but in 2012, clade B became predominant \[[@bib0020]\]. During the 2014 EV-D68 epidemic in Europe, most of the 205 sequenced viruses collected in 17 countries in the second semester of the year were found to belong to B1 and B2 clades, whereas A clades, although co-circulating, were detected in a smaller number of patients \[[@bib0155]\]. Since April 2016, EV-D68 infection outbreaks have been reported in several European countries, particularly in Northern Europe \[[@bib0155], [@bib0160], [@bib0165]\]. Phylogenetic analyses showed that these viruses belonged to the B3 lineage \[[@bib0160]\]. The same results were obtained with this study, suggesting that all the 2016 European outbreaks originated from a common source and that the B3 clade has emerged and become preeminent, eclipsing clades A and clades B1 and B2. The worldwide importance of clade B3 appears to be confirmed by the evidence that outbreaks of infections due to EV-D68 infection belonging to clade B3 have been reported in 2015 in China \[[@bib0170]\], and almost simultaneously with the European epidemics, in the USA \[[@bib0175]\].

In most of the cases, EV-D68 was detected as a single infectious agent. Moreover, when a coinfection was demonstrated, the co-infecting virus was found at a very low concentration, suggesting a bystander role. This result highlights the pathogenetic role of EV-D68, although a strict relationship between viral load and infection site or disease severity is not evidenced. The relatively low number of tested children with EV-D68 might explain this limitation.

In this study, EV-D68 infections were diagnosed mainly in June and July and in pre-school children. These findings agree with those reported by other authors, who found a higher incidence of EV-D68 diseases in warm months and in younger subjects \[[@bib0180], [@bib0185], [@bib0190], [@bib0195]\]. However, in this study, differently from other viral infections, particularly those due to RSV, the EV-D68-related respiratory infections were all of mild or moderate severity, and no patients showed signs and symptoms of nervous system involvement. Even when SaO~2~ levels were \<92%, the respiratory signs and symptoms rapidly subsided, and no patients were admitted to the intensive care unit. This result is in contrast to the findings reported by Dyrdak et al. who, among the 74 patients with EV-D68 clade B3 infection they identified, observed ten cases with severe respiratory or neurological symptoms and one death \[[@bib0165]\]. Knoester et al. also reported a high incidence of severe disease among children infected by EV-D68 \[[@bib0160]\]. These authors had to admit to the intensive care unit 13/17 (76.4%) of the children seen with respiratory disease. Moreover, in one case, AFP developed. By contrast, a child with EV-D68 B3 clade respiratory infections associated with acute flaccid paralysis has been recently described in Italy, in the same region and in the same period of the cases reported in this study \[[@bib0095]\].

In the past, it has been proposed that the development of a severe clinical picture following EV-D68 infection is the consequence of mutations occurring in the virus. In a study in which EV-D68 B1 strains associated with AFP cases diagnosed in the USA in 2014 were examined, 6 polymorphisms were identified and potentially associated with neuropathogenesis \[[@bib0200]\]. Other viral variants were identified in EV-D68 strains associated with severe disease in Taiwan \[[@bib0205]\]; moreover, in this case, it was proposed that the variants were associated with the severity of EV-D68. However, these findings have not been completely confirmed by other studies \[[@bib0100], [@bib0210]\]. By contrast, based on our strains sequences analyses, no evidence of increased clinical severity associated with specific molecular signatures in VP1 sequences has emerged from the present study. Further analysis on full-genome sequencing in a large cohort of positive patients with a well-define clinical syndrome could allow to correlated specific molecular signatures with clinical syndromes.

Based on adaptive evolution analysis, the overall dN/dS ratio (dN/dS \<1) of the VP1 gene observed in this study was similar to those reported by Lau et al. \[[@bib0215]\] Although several amino acid changes have been observed in VP1 sequences the meaning of dN/dS ratio reveal a purifying (negative) evolution for the VP1 protein. In addition, in keeping with this report, no codons with a significant signal of positive selective pressure were detected in the VP1 gene. By contrast, the phylogenetic analysis showed how the genetic changes may have driven the emergence of a new lineage of EV-D68 during the last 2 years. In fact, the worldwide circulation of EV-D68 strains belonging to the B2 and B3 lineages might have evolved from a common ancestor originating in 2011--2012 \[[@bib0220]\]. Recently, the evolutionary analyses performed on complete genome sequences have suggested that subclade B3 was not originated by a mutation occurred on subclade B1 strains. Indeed, it has been suggested a distinct evolutionary origin for the new B3 subclade \[[@bib0225]\]. Multiple clades of EV-D68 have worldwide circulated simultaneously between 2005 and 2014 \[[@bib0045]\], while a predominance of one clade (subclade B3) has recently been observed. Additional both epidemiological and molecular studies are needed to better clarify the reasons of this EV-D68 strain selection and widespread diffusion. It is interesting to note that the emergence of new EV-D68 lineages does not represent the circulation of a new virus. Therefore, diagnostic real-time assays, designed in a well conserved region (5′ UTR), are capable to detect all the new EV-D68 variants.

In conclusion, EV-D68 infections have been confirmed as a cause of lower respiratory illness in pediatric age. The circulation of one EV-D68 lineage has been proven in Italy as well as in the European region during 2016. However, further studies are required to investigate whether some strains or lineages may possess a higher affinity for the lower airway or central nervous system. In addition, the continuous global monitoring of the clinical and molecular epidemiology of EV-D68 should be improved to understand the factors determining possible spatiotemporal differences of EV-D68 infections.

Funding {#sec0060}
=======

This work was supported by the Ministero della Salute (Fondazione IRCCS Policlinico San Matteo, Ricerca Corrente grant 80622; Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico Ricerca Corrente grant 2017 850/02).

Competing interests {#sec0065}
===================

The authors declare that they have no potential conflict of interest.

Ethical approval {#sec0070}
================

The Ethics Committees of Fondazione IRCCS Policlinico San Matteo, University of Pavia (Italy) and the Pediatric Highly Intensive Care Unit of the Fondazione IRCCS Ca' Granda Policlinico, University of Milan (Italy) approved the study. Written informed consent was obtained from a parent of a legal guardian of each enrolled child. Children aged ≥8 years were asked for their assent.

[^1]: NA, not available; SaO~2~, oxygen saturation at hospital admission; ALL, acute lymphoblastic leukemia; HSCT, hematopoietic stem cell transplant; hRSV, respiratory syncytial virus; HRV, human rhinovirus; hAdV, human adenovirus; URTI, upper respiratory tract infection; VL, viral load.
