Increasing airline travel may facilitate co-circulation of multiple dengue virus serotypes in Asia

The incidence of dengue has grown dramatically in recent decades worldwide, especially in Southeast Asia and the Americas with substantial transmission in 2014–2015. Yet the mechanisms underlying the spatio-temporal circulation of dengue virus (DENV) serotypes at large geographical scales remain elusive. Here we investigate the co-circulation in Asia of DENV serotypes 1–3 from 1956 to 2015, using a statistical framework that jointly estimates migration history and quantifies potential predictors of viral spatial diffusion, including socio-economic, air transportation and maritime mobility data. We find that the spread of DENV-1, -2 and -3 lineages in Asia is significantly associated with air traffic. Our analyses suggest the network centrality of air traffic hubs such as Thailand and India contribute to seeding dengue epidemics, whilst China, Cambodia, Indonesia, and Singapore may establish viral diffusion links with multiple countries in Asia. Phylogeographic reconstructions help to explain how growing air transportation networks could influence the dynamics of DENV circulation.


Author summary
In the past half century the incidence of dengue fever worldwide has increased 30-fold, with an estimated~390 million infections per year. The years 2014 and 2015 were characterized by large dengue outbreaks worldwide, which are a threat to public health, especially in Asian countries. Here we use a phylogeographic approach to reconstruct historical virus movement and evaluate multiple potential predictors of the spatial spread of dengue virus (DENV) among Asian countries. We show that air traffic has played a more important role than maritime transport in shaping large-scale geographic dispersal of DENV in Asia. Our study suggests that the patterns of DENV spread result from the interplay between intensity and structure of human mobility through the air transportation network. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Introduction
Dengue virus (DENV) is a growing threat to public health, with nearly 390 million infections every year worldwide, of which~96 million are symptomatic [1,2]. An estimated 2.5 billion people are at risk of dengue infection [3]. Dengue is regarded as the world's most important mosquito-borne viral disease and is endemic in more than 100 countries [4], with most disease burden limited to tropical and subtropical regions [5]. However, in 2014 an outbreak of dengue occurred in Japan for the first time in over 70 years. This occurred despite the country's temperate climate, and viral phylogenetic analysis suggests that the Japanese outbreak resulted from international travel from Southeast Asia [6]. Furthermore, the number of reported dengue cases and outbreaks continues to increase. In 2014, a dengue outbreak affected several Asian countries, including China, Thailand, Vietnam and Japan [7,8].
The majority of DENV infections are asymptomatic or cause a mild febrile disease known as dengue fever, while the more severe forms of dengue infection-dengue hemorrhagic fever (DHF) and dengue shock syndrome (DSS)-may be life threating with >20% mortality [9]. Epidemiological studies have indicated that DHF often occurs when a dengue-immune person acquires a second infection with a different DENV serotype [10,11], and it has been hypothesized that DHF/DSS may result from a process of antibody-dependent enhancement [12,13]. The geographical areas in which transmission of multiple dengue serotypes occurs has grown in recent years and the pattern of co-circulation is conspicuously different from that which prevailed decades ago [14,15]. Furthermore, a growing number of DHF/DSS cases in the last 50 years, especially in Asia, demonstrates the need for a better understanding of how DENV genetic diversity and transmission jointly shape dengue epidemics.
The growing scale of human mobility, particularly through air transportation, underscores an increase of pathogen introductions into geographic areas suitable for transmission, potentially contributing to the emergence and the re-emergence of epidemics [16,17]. Recent examples include the SARS epidemic, novel influenza A virus strains, the Middle East respiratory syndrome (MERS) coronavirus, and Zika virus [18][19][20][21]. At local scales, endemic circulation of DENV strains can be driven by viral lineage introduction events that eventually lead to lineage replacement [22] and viral introductions are expected to increase with human mobility. Together, these observations suggest that increasing human mobility through airtraffic networks may impact the spread of pathogens at large geographical scales. Here we combine Asian socio-economic data, air and maritime transportation network data and viral genetic data to identify key drivers of the spread of DENV serotypes 1, 2 and 3 in Asia, and to reconstruct the patterns of virus circulation among Asian countries over the past half century.

Compilation of genetic datasets
Maximum phylogenetic information would be obtained by studying whole DENV genomes. However, available DENV whole genomes from Asia have insufficient coverage through time and space for reliable analysis, and most available sequences comprise partial or complete coding E gene sequences. In order to generate a data set with both acceptable phylogenetic diversity and spatiotemporal sampling, we used DENV E gene sequences in subsequent analyses. DENV (DENV-1 to DENV-3) envelope (E) gene sequences with known collection dates and locations of sampling in Asia were collected from GenBank. DENV-4 was not included in this study because too few samples were available (only 64 sequences from 11 countries; Fig 1). The remaining strains comprised a total of 2,202 sequences sampled between 1956 and 2015, from 20 distinct countries or geographic regions (Fig 1). Sequences were grouped by serotype and aligned separately using MAFFT [23]. Recombination was inspected using the methods implemented in RDP3 and SimPlot [24]. After removing duplicate and recombinant strains, the final data set contained 1,272 DENV-1 sequences, 628 DENV-2 sequences and 302 DENV-3 sequences. In the complete sequence dataset, some countries, such as Vietnam, Cambodia, Thailand, and Singapore, were over-represented. In order to control for possible bias from uneven sampling, we randomly subsampled the complete sequence datasets by location and sampling time. At most 10 sequences were sampled per country and per year in order to create a more equitable spatio-temporal sampling distribution. After sub-sampling, the total number of sequences analyzed here was 327 for DENV-1, 357 for DENV-2, and 202 for DENV-3, sampled over a total of 59 years (S1 Fig). Details of the sequences in each data set, including information on the year of isolation, sampling location, and accession numbers, are provided in Supplementary Information (S1 Table). Time-scaled phylogenetic tree reconstruction For each serotype subsampled data set, we first estimated the correlation between root-to-tip genetic divergence and sequence sampling dates, using TempEst [25]. This preliminary analysis indicated a good temporal signal for all serotypes (S2 Fig). To reconstruct past population dynamics, we used a coalescent-based Gaussian Markov random field (GMRF) method with the time-aware smoothing parameter [26], as implemented in BEAST v1.8.2 [27]. A GTR+I+Γ nucleotide substitution model and an uncorrelated lognormal relaxed molecular clock model were used, with a prior distribution for the evolutionary rate parameter set to a Γ distribution with shape = 0.001 and scale = 1000. The BEAGLE library was used to accelerate computation [28]. For each serotype, three independent analyses of 150 million generations were performed, sampling parameters and trees every 15,000 generations. Analyses were combined after the removal of a burn-in of 10-20% of the samples and were checked visually in Tracer v.1.5.

Estimating viral migration through discrete geographic locations
To reconstruct the spatial dynamics of DENV-1-3, we used a Bayesian Markov chain Monte Carlo (MCMC) phylogeographic discrete approach [29,30] that estimates the ancestral locations along each branch of a viral phylogeny, as implemented in BEAST v1.8.2 [27]. A GTR+ I +Γ nucleotide substitution model was used in this analysis. Three independent MCMC chains were run for 150 million states, sampling every 15,000 states after the removal of 10% burn-in. Maximum clade credibility (MCC) trees were summarized using TreeAnnotator and visualized using SPREAD [31].
To provide a minimal set of location state changes that provide an adequate description of viral spread, we used Bayesian stochastic search variable selection (BSSVS) [30,32]. BSSVS uses Bayes factors (BF) and binary indicator variables (I) to identify statistically supported viral lineage movement routes. Values of BF > 6 and I > 0.5 were considered as denoting a significant migration pathway, where BF > 1,000 indicates decisive statistical support, 100 BF < 1,000 indicates very strong support, 30 BF < 100 indicates strong support, 10 BF < 30 indicates substantial support and 6 BF < 10 indicates support [32]. We next estimated the number of expected transitions among location-pairs ("Markov jum" counts) along the phylogeny branches [33], which provided a quantitative measure of successful viral introductions among countries [34].

Air transportation network in Asia
We investigated the air transportation network in Asia using mobility data from ICAO (International Civil Aviation Organization; http://www.icao.int/Pages/default.aspx). This database contains the number of passengers traveling among 373 airports in Asia during the years 1982, and 1992-2012, and included all scheduled flights both for large and small aircraft. Countrylevel movement of passengers was obtained by aggregating airport-level movement for each country.
To obtain insight into temporal changes across the air transportation network in Asia, directed and weighted air flow networks were constructed, with countries as the network nodes. Network edges were weighted by the number of air passengers connecting pairs of countries. Hubs within the Asian air transportation network were also identified using two network topology properties: degree centrality and betweenness centrality. Degree centrality measures the number of edges connected to a node. Here, the degree of a given country refers to the number of airlines linking to it in the airline network. Betweenness centrality depends on the proportion of shortest paths between all pairs of vertices that pass through a given node. Thus the betweenness centrality of a given country is a measure of the extent to which a country lies on routes between other countries in the airline network and is calculated as: where σ st is the total number of shortest paths from node s to node t and σ st (k) is the number of those paths that pass through k.

Maritime transportation network in Asia
To represent the shipping connectivity between pairs of countries, the liner shipping bilateral connectivity index (LSBCI) in Asia from 2008 to 2014 was obtained from United Nations Conference on Trade and Development (UNCTAD; http://unctadstat.unctad.org/EN/Index.html), which captures the amount of goods transported between two countries. Container port throughput (CPT) for each Asian country was also obtained from UNCTAD.

Identifying potential determinants of DENV lineage movement
To investigate the potential predictors driving DENV-1, DENV-2 and DENV-3 spatial spread, we used a generalized linear model (GLM) extension to Bayesian phylogeographic inference [20,35]. The phylogeographic GLM method parameterizes the continuous-time Markov chain (CTMC) matrix of among-location lineage migration parameters as a log-linear function of several potential predictors. Each function includes a coefficient β (in log space) and a binary indicator variable δ. BSSVS [30] was used to estimate the posterior probability that each predictor is included or excluded from the model. We considered several potential predictors of DENV diffusion. These included log-transformed and standardized measures of demographic and economic data, geographical distances among countries, absolute latitude, DENV sample sizes. Finally, predictors of human mobility from shipping connectivity and air transportation movement were included (averaged values were estimated using data in 1982, 2000, and 2012). Specific details of potential predictors are as follows: (i) average distance between two countries was estimated as the average of pairwise distances between all pairs of airports in the two countries, (ii) absolute latitudes for each country were calculated as the latitudes of their geometric center, (iii) passenger flow represents the number of passengers on fights between each pair of countries, (iv) LSBCI (linear shipping bilateral connectivity index) represents goods transported by the sea between pairs of countries, (v) population sizes and densities for each country in 2000 were obtained from the UN World Population Prospects database (http:// www.un.org/en/development/desa/population/), (vi) gross domestic product (GDP) for each country in 2000 was collected from the UN Statistics Division (https://unstats.un.org) and (vii) average temperature and rainfall for each Asian country was extracted from the data set provided by WorldClim [36]. In addition, we included sample sizes (number of dengue sequences per country) for origin and destination locations as potential predictor variables, in order to test the impact of heterogeneous sampling. Fig 1A depicts the geographic locations of the sequences used in this study, from a total of 20 distinct countries or geographic regions in Asia; most sequences were sampled in Southeast Asia. To reduce potential bias in the reconstruction of spatial spread that may arise from oversampling particular locations, we subsampled the original data (see Methods). It is important to note that, in recent years, particularly since the 1980s, genetic evidence of multiple dengue virus sequences from a given location has increased (Fig 1B) in Asian countries (S3 Fig). Evolutionary analysis of 327 DENV-1 E gene sequences showed that Asian sequences fell into five distinct lineages, genotypes I-V. The 357 DENV-2 E gene sequences were classified into five genotypes: Asia I, Asia II, America/Asia, Cosmopolitan, and the America genotype.

Asian airline network growth
Air transportation has historically exhibited significant growth [37]. We found that Asian air transportation grew substantially from 1980 onwards (Fig 2A). We observe that several countries act as hubs in airline network, e.g. India, Singapore, and Thailand, while other countries exhibit growing network centrality over time, such as China and Malaysia (Fig 2B). It is interesting to note that the sequences obtained from dengue patients in 2014 in Japan shared a very high identity with a sequence from China, which shares one of the busiest Asian flight routes with Japan ( Fig 2B). These observations prompt us to question the role of air passenger transportation in DENV spread in Asia.

Air travel may drive the spatial spread of DENV in Asia
To infer the contribution of candidate factors driving DENV diffusion in Asia we used a generalized linear model that simultaneously estimates ancestral geographic reconstruction and identifies the contribution of potential predictors of spatial spread [20,35]. Our results show that most candidate predictors of virus spread, such as geographic distance and demographic factors, are not significantly associated with viral spread (Fig 3). However, we do find that air passenger flow is a dominant driver of DENV lineage movement and the inclusion of this factor in the model is supported for all three DENV serotypes. GDP at the location of origin and CPT (container port throughput) at the destination location are negatively associated with DENV-1 lineage movement in Asia, but not for the other two genotypes (Fig 3). Sample sizes for each location are not associated with viral lineage movement, which suggests that our conclusions are not driven by sampling biases.

Viral migration through discrete geographic locations
To understand the spatial circulation of DENV in Asia, we reconstructed the past spatial transmission patterns of serotypes DENV-1-3 over the study period, for each country from which sequences were sampled (Fig 4). A Bayesian stochastic search variable selection procedure was employed to infer a minimum set of location exchange events, while a Markov jump (MJ) analysis was used to quantity viral lineage movement among pairs of locations (see Methods). The phylogeographic analysis in Fig 5 indicates several significant migration links (with BF support > 1000). For DENV-1 these are between Cambodia and Vietnam, Thailand and Laos, and Singapore and China. For DENV-2, the well supported links are between Thailand and Cambodia, China and Indonesia, and Indonesia and Singapore. Finally, for DENV-3 the strongly supported links are between India and Sri Lanka (BF > 800). We also inferred a number of other well supported movements among countries; a full list is provided in S2 Table. We further find that estimated number of viral lineage migrations (including both importations and exportations) for each country is associated with measures of centrality in the air transportation network (Pearson correlation: r = 0.45, P = 0.05, for degree centrality and state transitions; and r = 0.73, P < 0.01, for betweenness centrality and state transitions; Fig 6A, S4  Fig). However, countries do not all contribute equally to viral lineage dissemination. An analysis of the number of virus lineage movements and the air transportation network analysis for each country indicates increased lineage virus movements from Thailand, India, and Indonesia to other locations [38,39], and a trend towards viral lineage importation from other locations for Vietnam and China. Lineage import and export is approximately equal for Cambodia and Singapore (Fig 6B and S5 Fig). Given the betweenness centrality of Thailand and India in the regional airline network, these countries may act as net sources of viral lineages, while China, Cambodia, Indonesia and Singapore may establish strong links with multiple countries within the network.
There has been an increasing frequency of DENV lineage co-occurrence in China over time, concomitant with its increasing centrality in the Asian air transportation network (Figs 2 & 5). In contrast, Vietnam, a dengue-endemic country, has a low frequency of viral lineage import and export with other locations (Fig 6B), likely due to its lower network centrality. However, we find that some non-endemic countries that do not contribute significantly to viral lineage movement across Asia may still maintain high airline passenger flows, such as Japan, South Korea, and countries in Western Asia. This may be because such areas lack a climate suitable for widespread dengue virus transmission. Further we are not able to infer transmission between countries if no viral sequences are available for those locations. The results here do not appear to be driven by heterogeneous sampling; similar results were obtained using downsampled datasets with a maximum of 5 sequences per country per year (S6 Fig).

Discussion
Our study provides a temporal description of the patterns of DENV spread across Asia. The spatial dynamics of dengue virus in Asia inferred here suggests that DENV genetic diversity in Asia is dynamic yet spatially structured, with frequent virus lineage movement among countries, and frequent co-circulation of dengue virus lineages. Our analysis suggests that the air transportation network has contributed to the spatial distribution of DENV serotypes in Asia. This can be explained by the movement of viraemic individuals; if the destination and timing of virus introductions coincides with a climatically suitable period then an epidemic can become established in a susceptible recipient population. Further, it is possible that the mobility of viraemic mosquitoes through air transportation networks is non-negligible [40]. For example, one study estimated that 8-20 Anopheline mosquitoes were imported into France per flight in one 3-week period in 1994 [41]. Previous studies have also shown that mosquito species, including Aedes albopictus, can survive long-haul flights [42,43]. Implementing disinfection is proving effective in vector control and disease prevention [40,44], highlighting the importance of air transportation in pathogen importation.
The positive coefficient β for the air traffic predictor in our phylogenetic GLM model indicates that DENV lineage dispersal is associated with air traffic; this result is consistent with previous studies of DENV in Brazil, which concluded that air travel of humans and/or mosquitoes is associated with virus lineage movement [22]. Notably there is strong statistical support Spread of dengue viruses in Asia for the inclusion of the air travel predictor variable in the GLM models for all three DENV serotypes in our study. Further, we find that origin GDP and destination CPT (container port throughput) are negatively related to DENV diffusion for DENV-1; this result is possibly linked to better public health infrastructure and vector control programmes in wealthier countries [45]. However, there was no support for the inclusion of those two predictors in the analyses of DENV-2-3. This may stem from the larger number of sequences for DENV-1, or from differences in sampling distribution among countries, despite the fact we have used a subsampling procedure in an attempt to mitigate potential bias. Alternatively, observed among-genotype differences may be driven solely by stochastic variation in the process of spatial dissemination; for example, the well-supported migration links inferred using DENV-1-3 sequences were not identical (Fig 5).

Spread of dengue viruses in Asia
Analysis of the air transportation network indicates differences in the roles of air transport nodes in the regional-scale dissemination of DENV lineages. Our data suggest that large transportation hubs account for most of the inferred virus lineage movements within their respective networks. For example, Thailand and India seem to act as central hubs of DENV lineage movement in Asia, demonstrating both inward/outward viral migrations with other countries in the region. In contrast, our results suggests that the low level of network centrality of Vietnam may explain its limited role in regional DENV geographical dissemination. Some countries, such as China and Malaysia, may have had an increasing role in DENV spread through time. We note that dengue is still considered to be an imported disease in mainland China [46]. It might be expected that some DENV endemic countries, with climatic suitability for dengue transmission and high betweenness-centrality in air transportation network, contribute to seeding Asian dengue epidemics. However, it has been shown that the risk factors for the persistence and transmission of dengue are complex [47], once the socioeconomic conditions were taken into account (as shown in Fig 3).
Previous studies based on simulation models have shown that global warming is likely to increase the area of land with a climate suitable for dengue virus transmission by altering the distribution of Aedes aegypti, the main mosquito vector of dengue [48,49]. This may result in more countries being at risk of dengue virus via the air transportation network, resulting in the potential for spread across even greater geographic scales. Temperature increases may lead to lengthened mosquito lifespans and shortened extrinsic incubation periods, which may result in more infected mosquitoes for a longer period of time [50,51].
Our study has several limitations. First, we cannot discern the relative contribution of infected humans versus infected mosquitoes in the spread of DENV, which requires more detailed epidemiological data and pathogen genome sequences from both mosquitoes and patients in future studies. However the information provided by the genetic analyses here can be useful to parameterize future spatially-structured transmission models. For example, the structure and travel flux of the airline transport network, as well as the among-country rates of Spread of dengue viruses in Asia lineage movement from the phylogeographic model, could be used to inform future simulations [20]. Second, our findings are based on DENV E genes for Asian countries, as this is the only genome region for which sufficient numbers of sequences are available, yet complete virus genomes would provide better phylogenetic resolution [52,53]. Further, our analyses were limited to DENV in Asia and do not consider transmission to or from other regions. Finally, the uneven sampling of DENV sequences in Asia, especially in a few countries after 2000 (e.g. Vietnam and Singapore) necessitated the use of sub-sampling analyses that accounted for the number of samples per location. Results obtained from a more representative data set indicates that our main conclusions are robust to sample sizes.
Future trends in global mobility could potentially accelerate the appearance and diffusion of DENV worldwide. Prevention and control of dengue epidemics requires a better understanding of its mode of geographical dissemination, especially for countries in the tropics. Our study highlights the importance of developing a multivalent DENV vaccine in order to cope with increasing frequency of DENV serotype co-occurrence. The potential impacts of vaccination on dengue epidemics in Asia should be considered [54]. Viral spatial dissemination, together with host age-structure and host-vector interactions are required in mathematical models to inform future vaccine deployment decisions.