Genetic Diversity and the Spatio-Temporal Analyses of Hantaviruses in Shandong Province, China

Hemorrhagic fever with renal syndrome (HFRS) is a serious public health problem in Shandong Province, China. We conducted an epizootiologic investigation and phylogeographic and phylodynamic analyses to infer the phylogenetic relationships of hantaviruses in space and time, and gain further insights into their evolutionary dynamics in Shandong Province. Our data indicated that the Seoul virus (SEOV) is distributed throughout Shandong, whereas Hantaan virus (HTNV) co-circulates with SEOV in the eastern and southern areas of Shandong. Their distribution showed strong geographic clustering. In addition, our analyses indicated multiple evolutionary paths, long-distance transmission, and demographic expansion events for SEOV in some areas. Selection pressure analyses revealed that negative selection on hantaviruses acted as the principal evolutionary force, whereas a little evidence of positive selection exists. We found that several positively selected sites were located within major functional regions and indicated the importance of these residues for adaptive evolution of hantaviruses.

Introduction

Hantaviruses (genus Orthohantavirus, family Hantaviridae, Order Bunyavirales) are among the most widely distributed groups of zoonotic viruses and have been detected in Asia, Europe, the Americas, and Africa (Mir, 2010). They are the causative agents of hemorrhagic fever with renal syndrome (HFRS) in Eurasia, and hantavirus cardiopulmonary syndrome (HCPS) in the Americas. Hantaviruses are negative-stranded RNA viruses with a tripartite genome, consisting of large (L), medium (M), and small (S) segments, encoding RNA-dependent-RNA polymerase (L), membrane glycoproteins (Gn and Gc), and nucleocapsid protein (N), respectively. Each hantavirus species is primarily associated with one or a few closely related rodent, insectivore, or bat species in nature (Zeier et al., 2005; Jonsson et al., 2010). Humans usually acquire hantavirus infection by contact or inhalation of aerosols and secretions from infected rodent hosts.

HFRS is endemic in China and remains a serious public health problem with about 10,000 human cases annually in recent years (Zhang et al., 2010, 2014; Jiang et al., 2016), even though many measures such as environmental management, host surveillance, and HFRS vaccine implementation have been carried out to control the disease. In China, the Hantaan virus (HTNV) and Seoul virus (SEOV), mainly carried by Apodemus agrarius and Rattus norvegicus, respectively, are the predominant causative agents of HFRS. The common symptoms of HFRS include headache, myalgia, abdominal and back pain, nausea, vomiting, and diarrhea. HTNV usually causes a more severe form of HFRS than SEOV does.

Shandong Province is located in eastern China and since 1980s, has been one of the areas most seriously affected by HFRS in China (Fang et al., 2010; Huang et al., 2012). Historically, it has had the largest HFRS burden in China–the cumulative number of human cases in Shandong Province accounts for one third of the national total (Kang et al., 2001; Fang et al., 2010). Although the number of human cases has been rapidly reduced within recent decades, it has remained at about 1,000–2,000 cases annually in Shandong, accounting for nearly one fifth of the total number of cases in China. An epidemiologic investigation suggested the presence of two pathogenic species of hantaviruses, HTNV and SEOV in Shandong Province (Li et al., 2007). However, little sequence information on hantaviruses from Shandong has been available to date, and there is limited knowledge on the spatio-temporal epidemiology and genetic characteristics of hantaviruses in Shandong Province. In this study, we conducted an epizootiologic survey, and phylogeographic and phylodynamic analyses to infer the viral phylogenetic relationships in space and time, and gain further insights into the evolutionary dynamics of hantaviruses in Shandong Province, China.

Materials and Methods

Ethics Statement

This study was reviewed and approved by the Ethical Review Board, Science and Technology Supervisory Committee at the Beijing Institute of Microbiology and Epidemiology. The animal work described here adhered to the guidelines of the Animal Subjects Research Review Boards in the Beijing Institute of Microbiology and Epidemiology.

Sample Collection and Processing

Rodents were captured in farm lands, grass lands, orchards, and residential areas from 2012 to 2015. Sample collection sites were indicated on a digital map of Shandong Province, which was provided by the Chinese Academy of Surveying & Mapping, using the ArcGIS 10.2 software (ESRI Inc., Redlands, CA, USA) (Figure 1). Cages with a treadle release mechanism were used for live trapping, according to previously described protocols (Mills et al., 1995). All animals were anesthetized with ether, and then sacrificed via cervical dislocation in the field owing to difficulties in transporting them. After the rodent species were identified by morphological assessment, they were necropsied on-site and lung tissues were collected. Instruments that had been cleaned with ethanol and heat-sterilized were used for each animal. Precautions for working with animals potentially infected with dangerous pathogens were strictly followed, and all processes were conducted in accordance with the National Monitoring Program for Hemorrhagic Fever with Renal Syndrome (Trial). Lung tissues were stored either in RNAfixer (Bio Teke Co., China) or liquid nitrogen. After they were transported to the laboratory, tissues were stored at −80°C freezer until further analysis.

FIGURE 1

Figure 1. Sample collection sites. Sample collection sites are represented by dots colored according to the disease etiology (red, Seoul virus; blue, Hantaan virus; gray, negative). Map was obtained from the Chinese Academy of Surveying & Mapping. Geographic coordinates of each sample collection site were obtained using the Baidu Coordinates Pick-up System, and were indicated on the map using the ArcGIS 10.2 software.

Lung tissue samples were screened using the reverse transcription-polymerase chain reaction (RT-PCR). Total RNA, which was extracted from 20 to 50 mg of lung tissue using the TRIzol reagent (Invitrogen, USA), was reverse-transcribed using M-MLV reverse transcriptase (Invitrogen, USA) and the primer P14: 5′-TAGTAGTAGACTCC-3′ (Wang et al., 2000). The PCR was performed to amplify a partial L sequence, using primer pairs as previously described (Klempa et al., 2006). The primers for S and M sequences were designed from consensus regions of SEOV and HTNV sequences (Table 1). The amplicons were sequenced using the ABI 3730XL Genetic Analyzer. The complete M and S sequences were generated from overlapping fragments using the Lasergene package (DNAstar, Inc., USA) and manually edited.

TABLE 1

Table 1. Primers used for nested and heminested PCR.

Sequence Comparison and Phylogenetic Analyses

Sequences were compared by the Lasergene package (DNAstar, Inc., USA) to determine the divergence of nucleotide and amino acid sequences.

The sequences were aligned using the Clustalx 1.8 software (Thompson et al., 1994). Phylogenetic trees were generated by the Bayesian Markov Chain Monte Carlo (MCMC) tree-sampling methods, implemented by the MrBayes 3.1 software (Ronquist and Huelsenbeck, 2003) using the GTR evolutionary model, with gamma-distributed rate variation across sites and a proportion of invariable sites. The run was stopped when the standard deviation of split frequencies was below 0.01. The tree was visualized using the FigTree v.1.4.1 software (http://tree.bio.ed.ac.uk/software/figtree/).

To gain a better understanding of the characteristics of the hantaviruses in Shandong Province, we reconstructed the phylogenetic relationships of SEOV and HTNV, including only those strains obtained in the province, based on M and S segment open reading frame (ORF) sequences using two different methods. Maximum clade credibility (MCC) trees were constructed by the Beast v 2.4.6 software (Bouckaert et al., 2014), using the model derived from the Modeltest 0.1 program (Posada, 2008) with the strict clock, relaxed clock log normal, and random local clock models. The convergence was assessed on the basis of the effective sampling size (ESS) using the Tracer v1.6 software (http://tree.bio.ed.ac.uk/software/tracer/). Only log-likelihoods with ESS values of > 200 were accepted. Models were compared by AICM values (Baele et al., 2012), and the lower AICM value was selected as the better model fit. The MCC phylogenetic tree was summarized by the TreeAnnotator program and visualized using the FigTree v1.4.1 program. Phylogenetic network analysis was performed using the median-joining (MJ) method with mutation codons. It was implemented by the Network 5.0.0 software (Bandelt et al., 1999) (http://www.fluxus-engineering.com/sharenet.htm) and invariant sites were removed from the dataset.

Phylogeny-Trait Association Analyses

The effect of geographic origin on hantavirus populations was evaluated by phylogeny-trait association analysis using the BaTS 2.0 software (Parker et al., 2008), based on M and S segment ORF sequences of SEOV and HTNV. Each BaTS analysis was performed using 180 Bayesian phylogenies, which were obtained using the LogCombiner program. The BaTS analysis of SEOV was run using three states (East, central and Jiaxiang). As for HTNV, we used four states (Northeast, Yiyuan, Jiaozhou and Huangdao). These states were determined according to the phylogenetic trees and geographic distribution of the strains. The BaTS results were reported with three scores: the association index (AI), parsimony score (PS), and monophyletic clade (MC). The AI and PS scores are two methods of testing the overall structure of all traits tested to tree topologies, whereas the MC scores indicate the association of specific traits with the tree topologies. Any score < 0.05 indicated an association for that particular test and score < 0.01 was indicative of strong association.

Molecular Diversity and Demographic Analyses

The nucleotide diversity (π) and Tajima's D test (Tajima, 1989) of the M and S segment ORF sequences for each geographical population were evaluated using the DnaSp 5.0 software (Librado and Rozas, 2009). A value of Tajima's D < 0 means that the population size is not at equilibrium, but expanding, due to an over representation of infrequent polymorphisms, which is typical after a bottleneck or a selective sweep and/or purifying selection.

The Star Contraction (SC) method (Forster et al., 2001) was used to calculate the SC network on the basis of the MJ network for SEOV. The SC network graphs were then compared to those of the MJ network, as this method can facilitate the analysis of historical demographic expansion events. We did not calculate the SC network for HTNV because of its relatively small dataset.

The coding regions of the M and S segments of each strain were concatenated in the order M–S. The population dynamics of hantaviruses were examined using the Bayesian skyline plot (BSP) (Drummond et al., 2005), as implemented in the BEAST program (Bouckaert et al., 2014) with the concatenated sequences. This method is a nonparametric coalescence analysis that uses standard MCMC sampling procedures to estimate past population dynamics from the posterior distribution of the effective population size (Neτ) (Drummond et al., 2002). The model estimation and selection were similar to that of the MCC tree construction.

Natural Selection Analyses

The Codon-based Z test of selection was performed using the MEGA 7.0 software (Kumar et al., 2016) to test the hypothesis the following about H0: (a) dN = dS (test of neutrality); (b) dN > dS (positive selection); and (c) dN < dS (purifying selection). Site-specific selection pressure was assessed using five varying methods, including single likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), internal fixed effects likelihood (IFEL), random effects likelihood (REL), mixed effects model evolution (MEME), and fast unconstrained Bayesian approximation (FUBAR), using the web-based interface Datamonkey (http://www.datamonkey.org/). Significance levels were set to p < 0.05, Bayes Factor (BF) > 100, and Posterior Probability (PP) > 0.95.

Results

Rodent Screening and Sequence Comparison

A total of 1,798 small animals representing 10 species were captured at 58 different collection sites in Shandong Province (Figure 1). They included R. norvegicus (n = 908), A. agrarius (n = 269), Mus musculus (n = 373), Sorex araneus (n = 73), Cricetulus triton de Winton (n = 31), Cricetulus Barabensis (n = 3) Niviventer confucianus (n = 25), Rattus flavipectus (n = 23), Cricetulus spp. (n = 92) and Microtus arvalis (n = 1). Among them, 111 samples generated the PCR products of expected size for the partial L sequence. Prevalence in individual counties ranged from 0 to 14.3%, and the highest prevalence was found in Qingzhou County (Table 2). Sequence comparison of the PCR products indicated that 93 strains from R. norvegicus and one from M. musculus were most related to SEOV, with 81.2–98.8% nucleotide similarity. Among them, 93 strains were similar to each other with more than 95.7% nucleotide similarity, and only one strain from R. norvegicus in Jiaxiang County (JX20141175) was divergent from most of the SEOV samples obtained in this study with < 83.3% nucleotide similarity. Seventeen samples from A. agrarius were similar to the HTNV with 80.9–85.6% nucleotide similarity.

TABLE 2

Table 2. RT-PCR detection of hantaviruses RNAs in rodents of Shandong Province.

A total of 76 complete M segment, and 71 complete S segment nucleotide sequences of SEOV were also determined. They were all obtained from R. norvegicus, with the exception of one strain from M. musculus. The M sequences of SEOV were 3,644–3,656 nt long. They all consisted of an ORF of 3,402 nt that encoded a putative 1,134 amino acid (aa) glycoprotein precursor with 95.3–99.3% nucleotide similarity to other SEOV strains. The complete S sequences of SEOV were 1,764–1,775 nt long. They all contained the single ORF that encoded the N protein of 430 aa with 95.7–99.4% nucleotide similarity to other SEOV strains. No differences were observed between the SEOV from R. norvegicus and that from M. musculus. We did not obtain the M and S sequences from the sample JX20141175.

A total of 14 complete M segment, and 13 complete S segment nucleotide sequences of HTNV were determined. The complete M sequences of HTNV were 3,605–3,616 nt long, containing a single ORF that encoded a putative 1,136 aa glycoprotein precursor. They were similar to other HTNV strains with 84.0–96.7% sequence similarity. The complete S sequences of HTNV were 1,701–1,703 nt long. All contained a single ORF that encoded a putative nucleocapsid (N) protein of 430 aa. It showed 85.8–97.0% nucleotide similarity to other HTNV strains.

Phylogenetic Analyses

Phylogenetic trees were constructed, including partial L segments (n = 94) (accession no. KY639712-KY639799), ORF sequences of M (n = 90) (accession no. KY639536-KY639625), and S segments (n = 84) (accession no. KY639626-KY639711) obtained in this study, and other representative sequences from different areas in East Asia downloaded from GenBank. Combined with the geographic origin of the strains, the results showed the distribution of SEOV throughout Shandong Province, whereas HTNV seemed to be co-circulating with SEOV only in the eastern areas (Zhaoyuan, Laixi, Pingdu, Jiaozhou, and Huangdao) and southern areas (Junan and Yiyuan) (Figures 1, 2). The phylogenetic trees indicated that most SEOV strains from Shandong clustered together with the strains from other areas in China, and the genetic variants from a single collection site clustered together (Figure 2). Only one strain from Jiaxiang (JX20141175) clustered together with the strains from Qingdao and formed a new lineage of SEOV in the L tree. They seemed to be a new subtype of SEOV. The HTNV strains from Shandong were clustered together with other HTNV strains and formed a new sister lineage (Figure 2). The phylogenetic trees also indicated geographic clustering of both HTNV and SEOV, especially for the M and S segments. This might be due to the increased length of the M and S sequences in the present study, which could have provided more information.

FIGURE 2

Figure 2. Phylogenetic trees estimated by the Bayesian method using the MrBayes software based on the 329-nucleotide alignment of L-segment sequences, and ORF sequences of the M segment and S segment. Phylogenetic trees were constructed using the GTR evolutionary model, with Dabieshan virus NC167 strain as the outgroup, gamma-distributed rate variation across sites and a proportion of invariable sites. The run was stopped when the standard deviation of split frequencies reached below 0.01. The L (accession no. KY639712–KY639799), M (accession no. KY639536–KY639625), and S (accession no. KY639626–KY639711) sequences of hantaviruses obtained in this study are shown in different colors according to their geographic origin. The numbers at each node are posterior probabilities (PP) based on 6,000,000–12,000,000 trees. Only PP values > 0.6 are shown. The scale bar indicates the nucleotide substitutions per site. SEOV, Seoul virus; HTNV, Hantaan virus; DBSV, Dabieshan virus; ASV, Amur-Soochong virus.

The Bayesian MCC trees of SEOV indicated that the strains from Shandong fell into four different lineages for both M and S segments (Figure 3), which were roughly consistent with their geographic origin. Lineage 1 included the strains from Jiaxiang County. Lineage 2 included the strains from the central areas (Linzi, Zichuan, Qingzhou, Yiyuan, and Anqiu), and one strain (1584325) obtained in Huangdao. Lineage 3 included strains from the eastern areas (Huangdao and Pingdu). Lineage 4 included only one strain from Huangdao. The MCC trees of HTNV strains in Shandong showed that they fell into four different lineages according to their geographic distribution (Figure 4), specifically, the strains from Huangdao, Yiyuan, Jiaozhou, and strains from the northeastern areas (Zhaoyuan, Laixi, and Pingdu). Interestingly, the strain ZY55 fell within the Huangdao cluster for the M-segment, and within the Northeastern cluster for the S segment. However, the notion that the ZY55 strain was a reassortment was not strongly supported by the present data, as the PP value for the stem was very low (0.44) in the M phylogenetic tree.

FIGURE 3

Figure 3. Relative phylogenetic relationships of M and S gene sequences for Seoul virus defined in MCC trees. Different colors represent the geographic origin of the samples. Numbers at each node are posterior probabilities (PP) based on 6,000,000–12,000,000 trees. Only PP values > 0.6 are shown. The scale bar indicates the nucleotide substitutions per site.

FIGURE 4

Figure 4. Relative phylogenetic relationships of M and S segment sequences for Hantaan virus defined in MCC trees. Different colors represent the geographic origin of the samples. Numbers at each node are posterior probabilities (PP) based on 1,000,000–5,000,000 trees. Only PP values > 0.6 are shown. The scale bar indicates the nucleotide substitutions per site.

The MJ network, constructed from 46 alleles for the M segment, demonstrated that three haplotypes (A, B, and C) were highly common and shared by various strains (Figure 5, SEO-M-MJ). Haplotype A included strains from the central areas, and was connected to several low-frequency tip haplotypes from those areas. This suggested that most strains in these areas shared the same ancestor. Haplotype B mainly included strains from Jiaxiang. Haplotype C included strains from Qingzhou County. The MJ network, constructed from 30 alleles for S segment of SEOV, indicated that only two haplotypes (A and B) were highly common and similar to those in SEO-M-MJ (Figure 5, SEO-S-MJ). The MJ network of HTNV showed no obvious haplotype trait.

FIGURE 5

Figure 5. Median-joining (MJ) and Star Contraction (SC) phylogenetic networks of Seoul virus and Hantaan virus. Each pie chart was colored according to the geographic origin of the samples and sized relative to its frequency in the dataset. The numbers in each pie chart represent the haplotypes of the strains. Small white circles represent median vectors (roughly equivalent to hypothetical unsampled haplotypes). Length of the branch is proportional to the number of mutational changes between haplotypes.

Phylogeny-Trait Association Analysis

Phylogeny-trait association analysis using the BaTS 2.0 software detected a very strong association between the geographic traits and phylogeny for both SEOV and HTNV (p = 0 for AI and PS statistics). When the MC statistic was used to test the extent of the phylogenetic clustering of SEOV in three regions (East, Central, and Jiaxiang), the results showed high significance (p < 0.01) among these geographic populations (Table 3). The MC statistic for HTNV (four regions) showed no significance in some populations (p > 0.05) (Table 3). Nevertheless, it should be noted that only a limited number of HTNV sequences were available in the present study, which might have led to false evaluation.

TABLE 3

Table 3. Region-phylogeny association parameters of hantaviruses in Shandong Province (three states for SEOV and four states for HTNV).

Molecular Diversity and Demographic Analysis

Neutrality tests performed on the M and S datasets produced negative Tajima's D values for most of the geographic populations (Table 4). Particularly in the central areas, the Tajima's D test value was significant for SEOV, which might be the main reason that these test results were significant for “all areas.” These findings support a model of the population expansion of SEOV in the central areas. In addition, when Tajima's D test was performed, multiple evolutionary paths were detected for the M segment of SEOV in the eastern areas (Huangdao and Pingdu), and for both M and S segments in the central areas.

TABLE 4

Table 4. Genetic diversity indices and demographic history parameters of hantaviruses in Shandong Province.

Historical demographic expansion events are characterized by star-like clusters of nodes around a founder population node. The SC algorithm can identify such clusters and shrink the nodes of a cluster back toward the founder node. We constructed the SC network and compared it with the MJ network to facilitate the analysis of the historical demographic expansion events. The results indicated two founder nodes (sco1 and sco2) for the M segment (Figure 5, SEO-M-SC). The sco1 node was formed by contraction of the strains from Linzi. The sco2 node was formed by strains from Qingzhou which were contracted within haplotype C. In addition, some strains from the central areas were contracted within haplotype A, and some strains from Jiaxiang were contracted within haplotype B (Figure 5, SEO-M-SC). Although no founder nodes were detected in the SC network of the S segment, some strains from the central areas were contracted within haplotype A (Figure 5, SEO-S-SC). These results suggest that there might be historical demographic expansion events for SEOV in the central areas. These findings are consistent with Tajima's D test results.

To explore the demographic history of the sampled population over the study period, the relative effective population size over time was inferred using the Bayesian skyline model, which essentially plots effective population size (Neτ) as a function of time. The results indicated that the effective population size of SEOV increased before the end of 2012, reached a plateau (2013–2015), and then rapidly decreased (Figure 6, SEO). The HTNV population remained roughly steady (Figure 6, HTN).

FIGURE 6

Figure 6. Population dynamics of SEOV and HTNV in Bayesian skyline reconstruction plots. The demographic inference of SEOV and HTNV genetic diversity is depicted using the Bayesian skyline reconstruction plot method. The plot was referred to the concatenated sequences in the order M–S. Solid lines indicate the median estimates, and the purple area displays the 95% highest posterior density (HPD). SEOV, Seoul virus; HTNV, Hantaan virus.

Selection Pressure Analyses

The selection pressure analyses using the Z test and codon-specific analyses revealed that negative selection was acting as the principal evolutionary force, with very low dN/dS values (far lower than 1), and an abundance of negatively selected sites (Table 5). The mean dN/dS values for Gc was 0.0393 and the dN/dS values for Gn was 0.0669. Codon-specific analyses also indicated little evidence of positive selection (Table 5). One positively selected site in the Gn/Gc protein and one in the N protein of SEOV showed significance (p < 0.05), based on the MEME method. Another positively selected site was found in the N protein of SEOV using the FUBAR method yielding a significance value (p < 0.05). Regarding the Gn/Gc protein of HTNV, three positively selected sites based on the REL method, and another positively selected site based on the FUBAR method all yielded a significance value (p < 0.05). No positively selected sites were detected in the N protein of HTNV.

TABLE 5

Table 5. Summary of selection pressures acting on Gn/Gc protein and N protein of SEOV and HTNV obtained in this study.

Discussion

In this study, we carried out an epizootiological investigation by collecting a large number of samples in 58 collection sites, covering the most severe affected HFRS endemic areas in Shandong Province, China (Fang et al., 2010; Wei et al., 2011; Cui et al., 2013; Wang et al., 2016). We found that both SEOV and HTNV were co-circulating in eastern and southern Shandong. However, the dominant hantavirus in the western and northern parts of Shandong was SEOV. Generally, HTNV is mainly associated with A. agrarius, which remains predominant in rural and forested areas, whereas SEOV is typically carried by R. norvegicus, which is distributed in almost all areas, especially the urban areas. Our results suggest that prioritization of control efforts should be focused on peridomestic rodents, such as R. norvegicus in residential areas throughout Shandong Province. M. musculus should also be monitored, as they could be affected by spillover infections with SEOV. However, the virus seems to be far less prevalent in M. musculus than in R. norvegicus. In eastern and southern Shandong, both the sylvatic and peridomestic rodents, R. norvegicus and A. agrarius should be monitored.

Phylogenetic analysis indicated that one strain from Jiaxiang (JX20141175) clustered together with two strains from Qingdao and formed a new lineage of SEOV in the L tree. This seemed to be a new subtype of SEOV that was widely distributed in Shandong. Unfortunately, we did not obtain the M and S sequences of the sample because of the high level of nucleotide divergence and limited availability of the tissue. Furthermore, no information on the M and S sequences of these strains was available in GenBank. Therefore, further epidemiological and virological studies in Shandong Province are necessary to evaluate the issue further.

Phylogeny-trait association analyses indicated very strong clustering of SEOV by geographic origin. This suggests that SEOV could be transmitted over a limited distance under normal conditions, possibly because its spread was accompanied by dispersion of the reservoir, which usually moves within a limited range. Therefore, the dispersal of SEOV should be slow and steady under wild conditions. However, when Tajima's D test was performed, some codons with multiple evolutionary paths were detected in the central areas. In addition, significant negative Tajima's D test value supported a model of population expansion for SEOV in the central areas. The SC analyses also suggested historical demographic expansion events associated with SEOV in the central areas. Furthermore, Bayesian skyline plot suggested an increase of the effective population size of SEOV. These results all suggest that SEOV in the central areas might have diverse origins and they might have experienced rapid and wide spread. This phenomenon was possibly related to its animal host, R. Norvegicus, which is spatially in closer proximity to human beings than any other hosts of hantaviruses. These hosts could migrate to other places further away by means of transportation. Therefore, SEOV could be transmitted over a relatively long distance by infected R. Norvegicus aboard motor vehicles (Plyusnin and Morzunov, 2001; Zuo et al., 2011). In addition, SEOV could be spread over a long distance via the fecal matter of infected R. norvegicus on cargo. These factors could facilitate the long-distance transmission of SEOV, which might be responsible for the appearance of SEOV in Europe (Plyusnina et al., 2012; Dupinay et al., 2014; Reynes et al., 2017). These factors also created favorable conditions for SEOV from different origins to enter the same area. A Combination of the highly frequent cargo trade in Shandong Province, with the construction of several highways likely presents suitable conditions for the transmission and dispersal of, and interaction between, SEOV strains, which is considered to be ongoing. Therefore, both SEOV and its hosts, R. Norvegicus, should be rigorously monitored, in order to control the disease in Shandong Province, particularly in freight terminals, wharfages and warehouses.

Our study indicates that hantaviruses in Huangdao show the greatest diversity, with two species of hantaviruses, and the highest nucleotide diversity for SEOV (0.0141 for the M segment and 0.0107 for the S segment) in Shandong. The present study also indicates that some codons with multiple evolutionary paths were detected for the M segment of SEOV in Huangdao when Tajima's D test was performed. Interestingly, the strain 1584325 from Huangdao shared the same ancestor with strains from the central areas. This finding suggests that this strain might have been transmitted from the central areas and could be regarded as an example of long-distance transmission of SEOV. Hantavirus diversity and complexity in the Huangdao district might be related to its geographical position. Huangdao is one of the emerging economic development zones located in Qingdao City, a port city with heavy traffic congestion. The city is affected by modes of transport over land and sea that move imports and exports on a daily basis. In addition, hantaviruses in this area seem to have undergone a relatively long evolutionary history, as a previous study indicated that Huangdao was one of the first locations at which hantaviruses appeared in Shandong (Fang et al., 2010). All of these factors increased the chances of the transmission of, and interactions among, the various strains of hantaviruses.

Understanding the driving forces acting on viral evolution and its associated mechanisms is one of central aspects of evolutionary biology. In addition, the study of these forces would facilitate the development of efficient strategies for the management of viral diseases. In this study, several methods of selection pressure analyses showed that the ratio of non-synonymous to synonymous substitutions for Gn/Gc and the N protein gene was very low, indicating that most amino acid changes were deleterious polymorphisms removed by negative selection. In comparison to the Gc gene sequence, the greater non-synonymous Gn sequence variation was consistent with membrane-distal localization and supported the notion that Gn was subjected to selective pressure of the humoral immune response (Li et al., 2016). Previous studies have also suggested that negative selection acts as the principal evolutionary force on SEOV and HTNV (Hughes and Friedman, 2000; Lin et al., 2012; Castel et al., 2014). Negative selection has important implications for the evolution of RNA viruses, because a large number of rare non-synonymous mutations resulting from viral replication are often purified from the population (Hughes and Hughes, 2007).

Based on our analyses, in addition to negative selection pressure, episodic diversifying selection could be another evolutionary force working on the hantavirus genomes, even though several sites were positively selected. We identified one positively selected site at aa 827 of the Gn/Gc protein, and two positively selected sites at aa 37 and aa 367 of the N protein for SEOV. The aa 35–38 of the N protein reportedly shows strong cross reactivity with, and is important for the recognition of, the SEOV antigen (Lindkvist et al., 2008). We identified three positively selected sites for the Gn/Gc protein of HTNV at aa 289, aa 373, and aa 889. The aa 882–896 is reportedly a B-cell epitope that displays neutralizing activity to HTNV infection (Yan et al., 2012). The higher proportion of positively selected positions located in the variable Gn/Gc envelope glycoproteins was consistent with their functional roles in the viral escape from immunological responses (Holmes, 2003). We noticed that none of these positively selected sites were identified by more than one method. Although the accurate functions of all mutated sites remain unknown, based on the literature, some evidently play important functional roles. Other positively selected sites might also have important functional roles. However, these findings do not necessarily mean that identification by one method was sufficient to conclude whether a site was indeed under positive selection. The natural selection profiles of the various sites should be further tested. We presumed that the mutated sites might provide a selective advantage for the virus to evade host immunity, and indicate the importance of the associated residues for the adaptive evolution of hantaviruses. These findings might be beneficial to the development of vaccines.

Author Contributions

S-QZ, X-JL, Z-QW, and W-CC conceived and designed the experiments. S-QZ, X-JL, Z-QW, and J-FJ collected the samples. S-QZ, X-JL, J-SZ, and Q-MZ performed the experiments. S-QZ, L-QF, and W-HZ analyzed the data. S-QZ, X-JL, and W-CC wrote the paper. All authors reviewed the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 81473025, 81673238, 30972522), the Basic Work on Special Program for Science & Technology Research (No. 2013FY114600), the Special Program for Prevention and Control of Infectious Diseases in China (No. 2013ZX10004218), the Natural Science Foundation of Shandong Province (No. ZR2016HM75), the National Key Research and Development Program (2016YFC1201304-04) and State Key Research Development Program of China (2016YFC1201902).

References

Baele, G., Lemey, P., Bedford, T., Rambaut, A., Suchard, M. A., and Alekseyenko, A. V. (2012). Improving the accuracy of demographic and molecular model clock model comparison while accommodating phylogenetic uncertainty. Mol. Biol. Evol. 29, 2157–2167. doi: 10.1093/molbev/mss084

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandelt, H. J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C. H., Xie, D., et al. (2014). BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10:e1003537. doi: 10.1371/journal.pcbi.1003537

PubMed Abstract | CrossRef Full Text | Google Scholar

Castel, G., Razzauti, M., Jousselin, E., Kergoat, G. J., and Cosson, J. F. (2014). Changes in diversification patterns and signatures of selection during the evolution of murinae-associated hantaviruses. Viruses 6, 1112–1134. doi: 10.3390/v6031112

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, F., Wang, T., Wang, L., Yang, S., Zhang, L., Cao, H., et al. (2013). Spatial analysis of hemorrhagic fever with renal syndrome in zibo City, China, 2009–2012. PLoS ONE 8:e67490. doi: 10.1371/journal.pone.0067490

PubMed Abstract | CrossRef Full Text | Google Scholar

Drummond, A. J., Nicholls, G. K., Rodrigo, A. G., and Solomon, W. (2002). Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data. Genetics 161, 1307–1320.

PubMed Abstract | Google Scholar

Drummond, A. J., Rambaut, A., Shapiro, B., and Pybus, O. G. (2005). Bayesian coalescent inference of past population dynamics from molecular sequences. Mol. Biol. Evol. 22, 1185–1192. doi: 10.1093/molbev/msi103

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupinay, T., Pounder, K. C., Ayral, F., Laaberki, M. H., Marston, D. A., Lacôte, S., et al. (2014). Detection and genetic characterization of seoul virus from commensal brown rats in France. Virol. J. 11:32. doi: 10.1186/1743-422X-11-32

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, L. Q., Wang, X. J., Liang, S., Li, Y. L., Song, S. X., Zhang, W. Y., et al. (2010). Spatiotemporal trends and climatic factors hemorrhagic fever with renal syndrome epidemic in Shandong Province, China. PLoS Negl. Trop. Dis. 4:e789. doi: 10.1371/journal.pntd.0000789

PubMed Abstract | CrossRef Full Text | Google Scholar

Forster, P., Torroni, A., Renfrew, C., and Röhl, A. (2001). Phylogenetic star contraction applied to Asian and Papuan mtDNA evolution. Mol. Biol. Evol. 18, 1864–1881. doi: 10.1093/oxfordjournals.molbev.a003728

PubMed Abstract | CrossRef Full Text | Google Scholar

Holmes, E. C. (2003). Error thresholds and the constraints to RNA virus evolution. Trends Microbiol. 11, 543–546. doi: 10.1016/j.tim.2003.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Yin, H., Yan, L., Wang, X., and Wang, S. (2012). Epidemiologic characteristics of haemorrhagic fever with renal syndrome in Mainland China from 2006 to 2010. WPSAR 3, 12–18. doi: 10.5365/wpsar.2011.2.2.007

PubMed Abstract | CrossRef Full Text

Hughes, A. L., and Friedman, R. (2000). Evolutionary diversification of protein-coding genes of hantaviruses. Mol. Biol. Evol. 17, 1558–1568. doi: 10.1093/oxfordjournals.molbev.a026254

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, A. L., and Hughes, M. A. (2007). More effective purifying selection on RNA viruses than in DNA viruses. Gene 404, 117–125. doi: 10.1016/j.gene.2007.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H., Du, H., Wang, L. M., Wang, P. Z., and Bai, X. F. (2016). Hemorrhagic fever with renal syndrome: pathogenesis and clinical picture. Front. Cell. Infect. Microbiol. 6:1. doi: 10.3389/fcimb.2016.00001

CrossRef Full Text | Google Scholar

Jonsson, C. B., Figueiredo, L. T., and Vapalahti, O. (2010). A global perspective on hantavirus ecology, epidemiology, and disease. Clin. Microbiol. Rev. 23, 412–441. doi: 10.1128/CMR.00062-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, D. M., Ruan, Y. H., Fu, J. H., Zhang, Z. B., Zhang, X. L., and Wand, K. A. (2001). Epidemic characteristic and its changing trend of hemorrhagic fever with renal syndrome in Shandong Province from 1990 to 1998. Chin. J. Epidemiol. 22, 475–476.

Klempa, B., Fichet-Calvet, E., Lecompte, E., Auste, B., Aniskin, V., Meisel, H., et al. (2006). Hantavirus in African wood mouse, Guinea. Emerg. Infect. Dis. 12, 838–840. doi: 10.3201/eid1205.051487

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Stecher, G., and Tamura, K. (2016). MEGA7: molecular evolutionary genetics is version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874. doi: 10.1093/molbev/msw054

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Zhao, Z., Wang, Z., Liu, Y., and Hu, M. (2007). Nucleotide sequence characterization and phylogenetic is of hantaviruses isolated in Shandong Province, China. Chin. Med. J. 120, 825–830.

PubMed Abstract | Google Scholar

Li, S., Rissanen, I., Zeltina, A., Hepojoki, J., Raghwani, J., Harlos, K., et al. (2016). A molecular-level account of the antigenic hantaviral surface. Cell. Rep. 15, 959–967. doi: 10.1016/j.celrep.2016.03.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Librado, P., and Rozas, J. (2009). DnaSP v5: a software for comprehensive is of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi: 10.1093/bioinformatics/btp187

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, X. D., Wang, W., Guo, W. P., Zhang, X. H., Xing, J. G., Chen, S. Z., et al. (2012). Cross-species transmission in the speciation of the currently known murinae-associated hantaviruses. J. Virol. 86, 11171–11182. doi: 10.1128/JVI.00021-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindkvist, M., Näslund, J., Ahlmb, C., and Buchta, G. (2008). Cross-reactive and serospecific epitopes of nucleocapsid proteins of three hantaviruses: prospects for new diagnostic tools. Virus. Res. 137, 97–105. doi: 10.1016/j.virusres.2008.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Mills, J. N., Childs, J. E., Ksiazek, T. G., Peters, C. J., and Velleca, W. M. (1995). Methods for Trapping and Sampling Small Mammals for Virologic Testing. Atlanta, GA: Centers for Disease Control and Prevention.

Google Scholar

Mir, M. A. (2010). Hantaviruses. Clin. Lab. Med. 30, 67–91. doi: 10.1016/j.cll.2010.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Parker, J., Rambaut, A., and Pybus, O. G. (2008). Correlating viral phenotypes with phylogeny: accounting for phylogenetic uncertainty. Infect. Genet. Evol. 8, 239–246. doi: 10.1016/j.meegid.2007.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Plyusnin, A., and Morzunov, S. P. (2001). Virus evolution and genetic diversity of hantaviruses and their rodent hosts. Curr. Top. Microbiol. 256, 47–75. doi: 10.1007/978-3-642-56753-7_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Plyusnina, A., Heyman, P., Baert, K., Stuyck, J., Cochez, C., and Plyusnin, A. (2012). Genetic characterization of Seoul hantavirus originated from Norway rats (Rattus norvegicus) captured in Belgium. J. Med. Virol. 84, 1298–1303. doi: 10.1002/jmv.23321

PubMed Abstract | CrossRef Full Text | Google Scholar

Posada, D. (2008). jModelTest: phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256. doi: 10.1093/molbev/msn083

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynes, J. M., Carli, D., Bour, J. B., Boudjeltia, S., Dewilde, A., Gerbier, G., et al. (2017). Seoul virus infection in Humans, France, 2014–2016. Emerg. Infect. Dis. 23, 973–977. doi: 10.3201/eid2306.160927

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronquist, F., and Huelsenbeck, J. P. (2003). MrBayes 3: bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. doi: 10.1093/bioinformatics/btg180

PubMed Abstract | CrossRef Full Text | Google Scholar

Tajima, F. (1989). The effect of change in population size on DNA polymorphism. Genetics 123, 597–601.

PubMed Abstract | Google Scholar

Thompson, J. D., Higgins, D. G., and Gibson, T. J. (1994). CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680. doi: 10.1093/nar/22.22.4673

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Yoshimatsu, K., Ebihara, H., Ogino, M., Araki, K., Kariwa, H., et al. (2000). Genetic diversity of hantaviruses isolated in china and characterization of novel hantaviruses isolated from niviventer confucianus and Rattus rattus. Virology 278, 332–345. doi: 10.1006/viro.2000.0630

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Wang, T., Cui, F., Zhai, S. Y., Zhang, L., Yang, S. X., et al. (2016). Hemorrhagic fever with renal syndrome, Zibo City, China, 2006–2014. Emerg. Infect. Dis. 22, 274–276. doi: 10.3201/eid2202.151516

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, L., Qian, Q., Wang, Z. Q., Glass, G. E., Song, S. X., Zhang, W. Y., et al. (2011). Using geographic information system-based ecologic niche models to forecast the risk of hantavirus infection in Shandong Province, China. Am. J. Trop. Med. Hyg. 84, 497–503. doi: 10.4269/ajtmh.2011.10-0314

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, G., Zhang, Y., Ma, Y., Yi, J., Liu, B., Xu, Z., et al. (2012). Identification of a novel B-cell epitope of hantaan virus glycoprotein recognized by neutralizing 3D8 monoclonal antibody. J. Gen. Virol. 93, 2595–2600. doi: 10.1099/vir.0.045302-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeier, M., Handermann, M., Bahr, U., Rensch, B., Müller, S., Kehm, R., et al. (2005). New ecological aspects of hantavirus infection: a change of a paradigm and a challenge of prevention–a review. Virus Genes 30, 157–180. doi: 10.1007/s11262-004-5625-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Wang, S., Yin, W., Liang, M., Li, J., Zhang, Q., et al. (2014). Epidemic characteristics of hemorrhagic fever with renal syndrome in China, 2006-2012. BMC Infect. Dis. 14:384. doi: 10.1186/1471-2334-14-384

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. Z., Zou, Y., Fu, Z. F., and Plyusnin, A. (2010). Hantavirus infections in humans and animals, China. Emerg. Infect. Dis. 16, 1195–1203. doi: 10.3201/eid1608.090470

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuo, S. Q., Fang, L. Q., Zhan, L., Zhang, P. H., Jiang, J. F., Wang, L. P., et al. (2011). Geo-spatial hotspots of hemorrhagic fever with renal syndrome and genetic characterization of seoul variants in Beijing, China. PLoS Negl. Trop. Dis. 5:e945. doi: 10.1371/journal.pntd.0000945

PubMed Abstract | CrossRef Full Text | Google Scholar