Genetic diversity of the enteroviruses detected from cerebrospinal fluid (CSF) samples of patients with suspected aseptic meningitis in northern West Bank, Palestine in 2017

Background Human enterovirus genus showed a wide range of genetic diversity. Objectives To investigate the genetic diversity of the enteroviruses isolated in 2017 in northern West Bank, Palestine. Study design 249 CSF samples from aseptic meningitis cases were investigated for HEV using two RT-PCR protocols targeting the 5’ NCR and the VP1 region of the HEV genome. The phylogenetic characterization of the sequenced VP1 region of Echovirus18 (E18) and Coxsackievirus B5 (CVB5) isolated in Palestine along with 27 E18 and 27 CVB5 sequences available from the Genbank were described. Results E18 and CVB5 account for 50% and 35% of the successfully HEV types, respectively. Phylogenetic tree of E18 and CVB5 showed three main clusters, with all Palestinian isolates uniquely clustering together with those from China and from different countries, respectively. Cluster I of E18, with 13 Palestinian and 6 Chinese isolates, showed the lowest haplotype-to-sequence ratio (0.6:1), haplotype diversity (Hd), nucleotide diversity (π), and number of segregating sites (S) compared to clusters II and III. Furthermore, cluster I showed negative Tajima’s D and Fu-Li’sF tests with statistically significant departure from neutrality (P<0.01). In both E18 and CVB5 populations, high haplotype diversity, but low genetic diversity was evident. Inter-population pairwise genetic distance (Fst) and gene flow (Nm) showed that the Palestinian E18 and CVB5 clusters were highly differentiated from the other clusters. Conclusions The study divulged close genetic relationship between Palestinian HEV strains as confirmed by population genetics and phylogenetic analyses.


Results
E18 and CVB5 account for 50% and 35% of the successfully HEV types, respectively. Phylogenetic tree of E18 and CVB5 showed three main clusters, with all Palestinian isolates uniquely clustering together with those from China and from different countries, respectively. Cluster I of E18, with 13 Palestinian and 6 Chinese isolates, showed the lowest haplotypeto-sequence ratio (0.6:1), haplotype diversity (Hd), nucleotide diversity (π), and number of segregating sites (S) compared to clusters II and III. Furthermore, cluster I showed negative Tajima's D and Fu-Li'sF tests with statistically significant departure from neutrality (P<0.01). In both E18 and CVB5 populations, high haplotype diversity, but low genetic diversity was PLOS  meningitis showing headache, neck stiffness, low grade fever, and CSF with or without pleocytosis in the absence of an acute bacterial pathogen [23]. All CSF samples were proven negative for classical bacterial pathogens by the hospital microbiology laboratory. Patients' demographic data and clinical history including age, sex, sign and symptoms, place of residence, date of onset of symptoms, and CSF laboratory test results were retrieved from the patients' files. The study was approved by the Ministry of Health in Palestine under reference number ATM/125/2013 that included the approval of using the CSF samples and the clinical history from anonymous children. Undocumented verbal consent was taken in the presence of the physician and parents/guardians. This included informing all parties that this study was approved by Ministry of Health and only leftover samples following routine laboratory diagnosis will be used, which would normally be discarded, in order to achieve study objectives. Physician gave the order to provide leftover samples for the study. No additional samples will be requested for the study and all remaining samples will be discarded immediately after the study according to protocol. The data was analyzed anonymously.

Viral RNA extraction
The HEV RNA was extracted from 200 μl of CSF using a QIAamp viral RNA Mini Kit (QIA-GEN, Valencia, CA, USA) according to the manufacturer's instructions.

RT-PCR
RT-PCR for HEV diagnosis. To identify the enterovirus associated aseptic meningitis cases, the HEV RNA was amplified using two primer sets targeting the 5'NCR of the HEV genome as described and modified previously [12,24]. Briefly, one step RT-PCR test was carried out in 25μl-reaction mixture containing 5μl viral RNA extraction, 10 U Reverse transcriptase (AMV), 10 pmol of the outer primers (P1; 5'-GTA CCT TTG TGC GCC TGT T-3' and P4; 5'-GAT GGC CAA TCC AAT AGC TA-3') and 12.5μl of PCR Reddy master mix (Thermo Scientific). Two microliters from the first round were further amplified in 25μl reaction mixture containing 10 pmol of the inner primers (P2; 5'-TGG CTG CGT TGG CGG CCT G-3' and P3; 5'-ACA CGG ACA CCC AAA GTA GTC GG) and 12.5μl of the PCR master mix (Thermo Scientific). In each PCR run, negative and positive controls were included. Five microliters of PCR product were analyzed by electrophoresis on a 2% agarose gel containing ethidium bromide and were visualized using the Gel Doc System 2000 (Bio-Rad Laboratories-Segrate, Milan, Italy). A band of 203 bp indicated a positive result.
RT-PCR for HEV genotyping. For HEV genotyping, an RT-PCR with two primer sets was used to amplify the 5' half of the VP1region of the viral genome as described previously [12,25]. Briefly, in the first round, a one-step RT-PCR was performed in 25 μl reaction mixture containing 20 pmol of the outer primers (EV-F; 5'-GYDGARACNGGVCACACRTC-3' and EV-R; 5' CTMATGAAHGGDATNGAYATBC-3'), 10 U AMV and 12.5 μl PCR master mix (Thermo Scientific). Two microliters from the first round was further amplified in 25μl reaction mixture containing 20 pmol of the inner two primers (ENTNES-F; 5'-GAYACWATG-CARACVMGRCAYGT-3' and ENTNES-R; 5'-GRGCAYTVCCYTCTGTCCA-3') and 12.5 μl PCR master mix (Thermo Scientific). Negative and positive controls were used in each run. A band of 400 bp visualized on agarose gel electrophoresis indicated a positive result. PCR amplicons were purified using the NucleoSpin Gel and PCR Clean-up from Machery Nagel (Germany) before sending for commercial sequencing.

Enterovirus typing and phylogenetic analysis
The HEV identity was investigated by searching the partial VP1 sequences obtained in this study in comparison with the HEV prototypes and other HEVs sequences available in GenBank using basic local alignment (BLAST, http://www.ncbi.nlm.nih.gov/BLAST). The HEV that showed more than 75% nucleotide similarity was assigned to be the same genotype. Enterovirus Genotyping Tool version 1.0 was used to type positive samples (http://www.rivm.nl/mpf/enterovirus/ typingtool/) [26]. Phylogenetic analyses of the RNA sequences of E18 and CVB5 genotypes were conducted based on the neighbor-joining (NJ) methods with Tamura-Nei and Kimura 2-parameter models implemented in MEGA 6 with 1000 bootstraps replicates for branch confidence in clades in each tree (Tamura et al., 2013). Poliovirus was used as an out-group.

Genetic diversity analysis and neutrality tests
To account for the limited number of genotype HEV cases (26) and the intermediate length of VP1 sequence (400bp) and the diversity that might be caused external factors such as amplifications and sequencing, several diversity indices were used [36]. The genetic diversity of the 13 E18 and the 11 CVB5 (9 in the present study and the 2 previously detected in 2013) detected in Palestine were analyzed based on VP1 gene. The 27 E18 from 10 different countries reported in the period of 1999-2017 and 25 CVB5 from 13 different countries reported in the period between 1971 and 2017 were retrieved from the Genebank and included in the analysis. EV18 and CVB5 RNA sequences were separately aligned using MEGA 6 [37]. Then, the genetic diversity of the VP1 gene in both E18 and CVB5 genotypes were calculated based on parameters including haplotype diversity (H d ), nucleotide diversity (π). Haplotype diversity (or gene diversity) refers to the number of haplotypes in the population. Nucleotide diversity is the average number of nucleotide differences per site in pairwise comparison between RNA sequences [38]. To detect departure from neutral theory of evolution (random mutation) at a constant population size, two neutrality tests were used. The first was Tajima's D test which compares the differences between the numbers of segregating sites (S) and the average number of nucleotide differences between two randomly chosen sequences from within in the population (K) [39]. The second test was Fu and Li's F test which compares differences between the number of singleton mutations and the average number of nucleotide differences between pairs of sequences [40]. DnaSP 5.1 software was used for the calculation at default settings. Although the number of isolates (sequences) is not fairly high, yet positive results were drawn with signs of genetic diversity.

Statistical analysis
Epidemiologic data were analyzed using EpiInfo statistical package. Analysis included distribution, 2x2 tables, and frequency tables. Fisher's exact test and Chi square with 95% confidence interval were calculated. The level of statistical significance was P<0.05.

Demography of circulating EVs
A total of 249 CSF samples were collected from hospitalized patients suspected of having aseptic meningitis in the period between March and October, 2017. Overall, 54/249 (22%) yielded positive results for HEV using RT-PCR targeting the 5'UTR. The general characteristic and the clinical history of the study samples are shown in Table 1.
Twenty-six (48%) were successfully genotyped by sequencing the amplicon of the amplified 5' VP1 region. Four different types of HEVs were detected. All of them belong to HEV-B species. Thirteen of the detected HEVs (50%) were E18, 9 (34.5%) were coxsackievirus B5 (CVB5), 3 (11.5%) were E25 and 1 (3.8%) was CVB2. Demographic data, the clinical history and the gene Bank accession number of the genotypes HEV cases are shown in Table 2.

Phylogenetic analysis of the partial VP1 gene of the E18 and CVB5
Phylogenetic analysis of the partial VP1 gene was conducted using the thirteen E18 and eleven CVB5 strains from Palestine (9 in 2017 and 2 in 2013). Twenty-seven E18 and 25 CVB5 sequences were retrieved from GenBank were included for comparison. The phylogenetic tree for E18 showed three main clusters with all Palestinian isolates uniquely clustering together along with those from China. Similarly, The CVB5 isolates were distributed into three clusters with Palestinian isolates in 2017 clustering together, along with isolates from different areas of the world, whereas, the two Palestinian isolates in 2013 grouped within the Predominantly-Chinese cluster (Fig 1)

Recombination analysis
The recombination analysis did not find any statistically significant evidence for recombination events between the aligned sequences, both on the Palestinian strains level (13 E18 strains and 11 CVB5) and international level that included Palestinian strains pooled with strains from the Gene bank (40E18 strains and 36CVB5). However, minimum number of recombination events was detected by DnaSP 5.1 software for E18 and CVB5. RDP4 software (PHI statistic) revealed good chance of recombination, but at a very low P-value (0.00001).

Diversity indices
Population nucleotide diversity indices and neutrality tests were calculated for the partial VP1 gene for E18 and CVB5, separately, based on phylogenetic clusters (Tables 3 and 4). The totally haplotype diversity (Hd) for the 40 E18 sequences was 0.98± 0.01 and 0.99± 0.006 for the 36CVB5. At the same time, the total genetic diversity (π) for E18 was 0.12± 0.02 and 0.17± 0.01 for CVB5. The average number of nucleotide differences between any two sequences (k) for E18 and CVB5 were 35.5 and 41, respectively. Cluster I of E18 showed peculiar results compared to the other two clusters. In this group, we detected 12 haplotypes in 19 sequences with    (Tables 3 and 4). Both E18 and CVB5 showed a combination of high haplotype diversity (Hd), but low genetic diversity (π). Also, both showed overall negative values of neutrality tests; Tajima's D and Fu-Li's F. Inter-population pairwise genetic distance (Fst) in the three E18 populations ranged from 0.39 to 0.63 with Nm value from 0.29 to 0.78 (Tables 5 and 6). Fst for cluster I containing all Palestinian strains compared to clusters II and III were high (0.51 and 0.63) (positive Tajima's D) indicating population differentiation with low gene flow, Nm (0.48 and 0.29). However, genetic differentiation between clusters II and III is low (Fst = 0.39) with high Nm (0.78), but still reflected genetic differentiation. This is supported by the negative value of Gst (-0.006).
As for CVB5, cluster III which contained most of the Palestinian strains is differentiated from the other two clusters, I and II, as reflected by the high Fst and low Nm values (Fst = 0.54, Nm = 0.43) and (Fst = 0.49 and 0.52), respectively. At the same time, the genetic differentiation between clusters I and II is relatively low as supported by very low Gst value (Gst = 0.00005) and relatively low Fst (0.32) with high gene flow (Nm = 1.06).

Discussion
So far, only one report in Palestine confirmed the isolation of 7 different enteroviruses genotypes including echovirus 13,14,9,30,16,6 and coxsackievirus B5 from sporadic cases of aseptic meningitis and/or sepsis like illness [12]. In the present study, all of the HEV positive cases occurred in children less than 7 years old and 50% of them < 1-year-old. Highest rate of enterovirus-positive cases was reported previously in children less than 1 year of age in Palestine, United States, and Korea [12,42]. On the contrary, other studies, reported that HEV aseptic meningitis occurs most frequently in patients with age range of 3-12 years' old [15][16][17]43]. The discrepancy in the HEV age groups could be due to the source of cases whether sporadic or from an outbreak.  The present study showed that E18 and CVB5 were the most predominant genotypes, representing, 50% and 34.6% respectively. Several recent studies in Germany, Taiwan, Korea, and China reported that both E18 and CVB5 were the most predominant HEVs reported from sporadic and outbreak cases of aseptic meningitis [15,[17][18]44].
All E18 sequences included in the phylogenetic analysis grouped into three major genetic clusters. The E18 isolates from Palestine clustered together along with those from China reported in 2015 in cluster I. Clusters II and III contain E18 from Australia, France, USA, Sweden, Russian, Germany and Netherlands reported in the periods 2005-2012 and 1999-2012, respectively. Other minor clusters included cluster IV and V which included the prototype Metacalf and few HEV reported from India and China in the period 2000-2008. The low number of E18 isolates in clusters IV and V is due to low number of the E18 included in the phylogenetic analysis. Similar results reported a genetic divergence in the complete VP1 gene of the E18 that resulted in the formation of five phylogenetic clusters [44]. These viruses were reported recently from China, India, South Korea, Australia, Netherlands, Germany, Sweden, Russia and France [44]. Phylogenetic analysis of partial or complete VP1 gene of E11, E30, E13, and E6 of the HEV-B species revealed high genetic diversity showing several clusters and sub-clusters [44]. Accordingly, such data indicate that the same HEV genotype may circulate in different geographies at different times. Therefore, and due to lack of reporting HEVs in Palestine, the possibility of knowing whether E18 isolates in cluster I are recently imported or have already been circulating before remains a dilemma.
The CVB5 sequences grouped into three genetic clusters (Fig 1). Cluster III contained the 9 Palestinian strains isolated in 2017 along with CVB5 from Brazil, France, Cyprus, Germany and Australia reported in the period from 2006 to 2017. The 2 Palestinian isolates in 2013 grouped with a more recent CVB5 in cluster II (2003-2014) from China, Denmark, France and Netherland. Cluster I contained the prototype Faulkner and isolates from France, Germany, Japan, Korea, China, Hong Kong, India, Finland and Sweden. Recently, few studies compared the partial or the complete VP1 gene revealing genetic diversity of CVB5 in 2-5 clusters [45][46][47]. Accordingly, the findings in the present study reaffirm that two or more clusters may co-circulate and co-evolve in the same region as a result of the genetic diversity forces such as mutation or recombination.
The recombination for both species E18 and CVB5 was minimal, if at all present between the study sequences which could be explained by the infrequency of recombination within the capsid region, VP1, low level of co-infection, different tissue tropism, low level of viral prevalence, low geographic differences in the occurrences of virus, high nucleotide similarity in viral genome, presence of transcriptional factors and proteins like Cis-acting elements [21][22]48]. However, it is worth mentioning that recombination analysis was based on partial VP1 sequence (400 nt) and not on complete sequence of the gene. This puts forward mutation rates as the main cause of genetic diversity. The significant departure from neutrality as confirmed by negative Tajima's D and Fu-Li's F accompanied by high haplotype (Hd) and low nucleotide diversity (π) for cluster I of E18 (Table 3) may suggest recent rapid population expansion phase following bottleneck or genetic hitchhiking (genetic draft or gene sweep) that bring excess number of rare alleles. This was supported by the negative values Tajima's D and Fu-Li's F tests in the individual populations and the overall value for the three populations. A study in Taiwan showed that partial VP1 genes of E18 have low genetic diversity with high similarity between regions [16]. The clustering of the Palestinian and Chinese isolates in cluster I could have been brought about by the activity of Palestinian traders to China and back. Cluster II and III had similar diversity indices, but without any significance. Cluster III of CVB5 showed lower genetic diversity than the others (Table 4), yet insignificant; suggesting that cluster III may have undergone a neutral or contraction period or may be due to small subpopulations with limited number of sequences to yield sufficient statistical power. A study from Korea revealed very high similarity in the partial VP1 sequences among Korean samples (intra-population similarity) and between Chinese strains (inter-population similarity) [18]. The genetic similarity of the Palestinian, Korean and Taiwanese E18 and CVB5 strains to those from China, may hint to their dispersal from China to the rest of the world [16,18]. However, a more extensive study of several endemic areas should be conducted to confirm the ancestral origin of E18 and CVB5, the route of spread and the distribution.
The estimations of inter-population comparison indices (Fst, Nm, Kxy, Dxy, Gst and Da) ( Tables 5 and 6) support high level of genetic differentiation between the three main clusters of both E18 and CVB5. These results are supportive of the cluster in phylogenetic analyses (Fig  1).
In conclusion, the present study unravels three main clusters for each of the E18 and CVB5 with both showing high haplotype diversity compared to lower genetic diversity. The Palestinian isolates grouped mainly into one cluster for each viral species. Finally, present study supports close genetic relationship between Palestinian HEV species as confirmed by population genetics and phylogenetic analyses.