Evolutionary study and phylodynamic pattern of human influenza A/H3N2 virus in Indonesia from 2008 to 2010.

Influenza viruses are by nature unstable with high levels of mutations. The sequential accumulation of mutations in the surface glycoproteins allows the virus to evade the neutralizing antibodies. The consideration of the tropics as the influenza reservoir where viral genetic and antigenic diversity are continually generated and reintroduced into temperate countries makes the study of influenza virus evolution in Indonesia essential. A total of 100 complete coding sequences (CDS) of Hemagglutinin (HA) and Neuraminidase (NA) genes of H3N2 virus were obtained from archived samples of Influenza-Like Illness (ILI) surveillance collected from 2008 to 2010. Our evolutionary and phylogenetic analyses provide insight into the dynamic changes of Indonesian H3N2 virus from 2008 to 2010. Obvious antigenic drift with typical 'ladder-like' phylogeny was observed with multiple lineages found in each year, suggesting co-circulation of H3N2 strains at different time periods. The mutational pattern of the Indonesian H3N2 virus was not geographically related as relatively low levels of mutations with similar pattern of relative genetic diversity were observed in various geographical origins. This study reaffirms that the existence of a particular lineage is most likely the result of adaptation or competitive exclusion among different host populations and combination of stochastic ecological factors, rather than its geographical origin alone.

Introduction

Influenza A viruses cause acute respiratory disease in humans and are responsible for annual epidemic with 3–5 million hospitalization and approximately 500,000 deaths globally [1]. In temperate countries, influenza viruses activities typically occur in the winter season, characterized by the increasing number of influenza cases and associated deaths above a seasonal baseline, known as seasonal influenza pandemic [2]. In contrast, in the tropics, influenza occurs throughout the year without a well-defined pattern, although there is some evidence of higher frequency of respiratory viral infections in rainy season [3].

Two major surface glycoproteins of influenza A viruses—the Haemagglutinin (HA) and Neuraminidase (NA), are the main targets of neutralizing antibodies. Therefore, genes encoding these glycoproteins undergo most frequent mutations. The gradual accumulation of amino acid leads to changes over time within the antigenic sites of HA and NA, referred to as antigenic drift, allowing the virus to escape neutralization by pre-existing antibody [4]. Antigenic drift has contributed to the persistence of epidemics in the community due to inefficient immune clearance by hosts, the occurrence of antiviral resistance [5, 6], implying the need for periodical update of vaccines [7].

Influenza A/H3N2 and Influenza A/H1N1 subtypes together with influenza B have circulated in human population and become the source of annual epidemics in temperate countries. Among those viruses, influenza A/H3N2 subtype has caused frequent outbreaks with associated severe illness [7, 8], and also shows the strongest antigenic drift [4, 9]. The composition of WHO recommended vaccine has been updated several times since 1999 due to recurrent changes of influenza A/H3N2 virus [7].

Indonesia with approximately five million square kilometres of territorial area is the largest country in Southeast Asia. It has more than 17,000 islands and is home to around 250 million people with more than 300 distinct ethnic groups [10, 11]. The bigger western islands of Java, Bali, Sumatra, Kalimantan (Borneo), and Sulawesi are the most populated with faster socio-economic development than the small islands of eastern part of Indonesia. In addition, there is unique ecological and biological arrangement across geographical regions of the country [12]. This diversity is associated with varied resistance and susceptibility to diseases, which could be reflected in the diversity of pathogens due to host-agent interaction among different ethnic populations in Indonesia [13, 14].

The endemicity of poultry H5N1 and high number of human H5N1 cases [15] make Indonesia as an important country associated with the dynamic changes of influenza A viruses. The proximity of domestic animal and human in Indonesia raises concern for the generation of new influenza viruses through possible reassortment of H5N1 with other circulating influenza viruses [16]. The newly generated virus might further obtain beneficial multiple mutations allowing it to adapt and spread among humans. It has been suggested that tropical zone may function as the permanent mixing pool for viruses from around the world [17].

Most of the data concerning seasonal appearance and epidemic of influenza A in human come from temperate countries, and very limited data have been reported on seasonality of influenza in tropical region [18]. With influenza A/H3N2 virus infection occurring throughout the year without well-defined pattern [2], the above situation necessitates studying the dynamics changes and evolution of this virus across years and geographical regions in Indonesia. The present study describes the evolutionary analyses of HA and NA genes of Influenza A/H3N2 virus in Indonesia from 2008 to 2010.

Materials and methods

Population study and sampling strategy

The study samples were 242 archived clinical specimens of nasal swabs with confirmed H3N2-positivity obtained during the Influenza-Like Illness (ILI) Surveillance, which was conducted by the National Institute of Health Research and Development (NIHRD), Indonesia, from 2008 to 2010. All ILI samples were collected from outpatients that showed ILI symptoms according to WHO protocols. The confirmation was performed by NIHRD laboratory during the surveillance period following the WHO protocols [19]. All samples were suspended in Hank Balance Salt Solution (HBSS) transport medium and stored properly in -80°C. This study was approved by the Ethics Committee of NIHRD with Federal Wide Assurance (FWA) No. IORG0002751/IRB00003331. All samples were fully anonymized and could not be traced back to individual participants of the study.

RNA isolation, RT-PCR and direct sequencing

Viral RNA was extracted directly from 140 μL of sample using QIAmp Viral RNA Minikit (Qiagen, Hilden, Germany) according to the manufacturer’s instruction. RT-PCR and direct sequencing to obtain complete coding sequence (CDS) of HA and NA genes were conducted according to previous study [20].

Dataset preparation

The dataset consisted of sequences generated in this study and sequences retrieved from GenBank that were significant and relevant to be included in the analyses. As illustrated in S1 Fig, the selection of the sequence database was carried out through the following steps:

First, previously published sequences of HA and NA genes of H3N2 from other countries were obtained from GenBank in ‘GenBank-sequence’ format as per July 2, 2013 (10,972 HA and 5,530 NA sequences). The retrieved sequences were then parsed to yield HA-fasta and NA-fasta files containing sequence data with labels indicating isolate accession number, location, and date of isolation as fractioned years.

Second, the sequences from first step were then compared with Indonesian sequences using Usearch, which is a new algorithm for sequence database search that can perform high-scoring of local and global alignment [21]. This process resulted in UC files for HA and NA sequences. The sequences that had 90% similarity were included for the subsequent bioinformatics analysis together with Indonesian samples.

In the third step, the sequences from step two and the Indonesian sequences were aligned using MUSCLE [22]. Maximum likelihood tree of Indonesian and sequences database were constructed by FastTree [23, 24] to find the sequences that are close to Indonesian sequences. Database sequences that did not belong to the same branch as the Indonesian sequences and the non-significance database sequences were removed using patristic analysis [25]. This step resulted in datasets for HA and NA genes that were further analyzed by BEAST and MrBayes methods.

Phylogenetic analysis, bayesian skyline analysis, TMRCA, and substitution rate estimations

Bayesian Markov Chain Monte Carlo (MCMC) method in BEAST package version 1.6.2 was used to estimate the rate of nucleotide substitutions per site per year, The Most Recent Common Ancestor (TMRCA) and relative genetic diversity (expressed as Nеτ, where Ne is the effective population size and τ denotes the generation time of host-to-host infection) [26]. The general time reversible GTR+I+Г4 model of nucleotide substitutions determined by jModelTest prior BEAST analysis [27], relaxed molecular clock model, and Bayesian skyline based on coalescent approach were used for the analysis [28]. We ran MCMC for 100 million generations with sampling every 1000 generations to produce at least 100,000 trees to ensure adequate sample size of the posterior, prior, nucleotide substitution rates and likelihoods (effective sample size >200). The mean substitution rates, TMRCA, and maximum clade credibility (MCC) phylogenetic tree were calculated after 10% removal as burn-in following visual inspection with TRACER version 1.5. The remaining 90.000 sampled trees were summarized using TreeAnnotator to infer the MCC trees, followed by visualization with FigureTree version 1.3.1 where the posterior probability values shown on the major branches [26]. Vaccine strains from 1968 to 2010 and other sequences from Northern Hemisphere, Southern Hemisphere, and tropical countries were included and analysed together with the Indonesian samples. The sequence data were divided into nine taxon groups based on the density and movement of the people and the access for transportation between the respective origins as follows: a) Non Indonesia, b) Indonesia, c) Balikpapan, d) Banjarmasin, e) Batam, Medan and Aceh, f) Java and Lampung, g) Jayapura, h) Makassar, and i) Merauke, as shown in Fig 1.

Fig 1

Indonesian map with color-coded geographical origin.

Selection analysis

Selection analysis was undertaken using MrBayes version 3.1.2. using a GTR evolution model coupled with gamma rates and codon model, with one million iterations and subsampling of every 100 iterations. The analysis of selective pressure on HA and NA glycoproteins was estimated using MrBayes by looking at the ratio of non-synonymous (dN) to synonymous (dS) substitutions per site, referred as ω = dN/dS. The mean probability of 0.5 was set as threshold to determine positive selection pressure [29].

Results

Sequence dataset and phylogenetic analysis

One hundred of the 242 samples were successfully sequenced for the CDS of both HA and NA genes. Based on the originating places, the sequences were listed according to the main biogeographical (geographical, ethnological, and zoological) regions of the Indonesian archipelago (Table 1) [12]. The sequences and origins of the samples are shown in S1 Table and have been deposited in Global Initiative on Sharing All Influenza Data (GISAID) with accession number EPI465537-EPI465700, and EPI465712-EPI465748.

Table 1

Sample distribution based on biogeographical origins.

Year	Sumatra (Banda Aceh, Batam, Lampung, Medan)	Java	Kalimantan (Banjarmasin, Balikpapan)	East Timor and Bali (Kupang, Denpasar)	Sulawesi (Makassar)	Papua (Merauke, Jayapura)	Total
2008	15 (8)	35 (17)	6 (5)	-	34 (12)	21 (9)	111 (51)
2009	10 (5)	4 (3)	8 (6)	4 (2)	7 (3)	1 (1)	34 (20)
2010	33 (7)	25 (4)	8 (6)	5 (1)	28 (6)	8 (5)	107 (29)
Total	58 (20)	64 (24)	22 (17)	9 (3)	69 (21)	30 (15)	252 (100)

Number in bracket represents the number of the sequenced samples

The sequences were analysed together with 320 and 396 sequences (S2 Table), which were selected from 2,364 and 2,495 non-redundant sequences of HA and NA, respectively, of the published sequences in the corresponding years (see Materials and Methods). The phylogenetic tree with maximum clade credibility (MCC) of HA and NA are shown in Figs 2 and 3, respectively. All sequences were grouped into three major lineages: I, II and III, each with multiple sub-lineages. Lineage I consisted of sequences mostly collected in 2007 and 2008 together with the Indonesian sequences of 2008. All strains from this lineage did not continue their lines further than 2008, with the exception of one sample (Batam260). Lineage II mostly consisted of 2009 sequences, while lineage III was composed of 2010 sequences. The topology of the three linages showed the typical ‘ladder-like’ phylogeny with replacement of old strains by newer ones, with each lineage leaving a trunk that became the ancestor of the next generation.

Fig 2

MCC tree from Bayesian analysis of HA gene of H3N2 viruses originated from Indonesia and other countries.

Indonesian samples are illustrated in different colored branches according its geographical origin. The numbers I, II and III represent major lineages of the virus. Insert figure is the color-coded geographical origin.

Fig 3

MCC tree from Bayesian analysis of NA gene of H3N2 viruses originated from Indonesia and other countries.

Indonesian samples are illustrated in different colored branches according its geographical origin. The numbers I, II and III represent major lineages of the virus. Insert figure is the color-coded geographical origin.

There were two sub-lineages in Lineage III: (a) sub-lineage 1, consisting of two sequences (Makasar917 and Cirebon1828) collected in 2008, and four sequences (Banjarmasin079, Banjarmasin097, Banjarmasin064 and Kupang105) collected in 2009, as well as some from 2010; (b) sub-lineage 2, which mainly consisted of Indonesian sequences from 2010 and one from 2009 (DKI94). Some variants that were distant from the trunk in an evolutionary terminal or dead-end clade were also found.

Bayesian skyline analysis

The Bayesian skyline plot, which illustrates the changes in genetic diversity through time [30], showed the changing level of genetic diversity of HA and NA genes of Indonesian H3N2 viruses sampled from 2008 to 2010 (Fig 4). There were variations of genetic diversity across years showing three distinct increases with different levels of relative diversity, both in HA and NA genes. The first increase of HA Bayesian skyline plot was at the end of 2008, whereas the second increase was at the end of 2009, followed by rapid decrease at the beginning of 2010. The third increase of genetic diversity occurred from March to June 2010. In NA genes, the first increase of genetic diversity occurred at the end of 2008. The second increase was identified as a small peak from September to October 2009, slightly earlier compared with HA Bayesian skyline plot. The third increase was observed from May to June 2010.

Fig 4

The Bayesian Skyline plot of HA (A) and (B) NA genes. The x-axis represents the time (mm/yy) and the y-axis represents a measure of relative genetic diversity (Net, where the Ne is the effective population size and t the generation time from infected host to infected host). Black line represents the mean value while the 95% confidence limits shown in blue line.

To determine the detailed Bayesian skyline plots according to geographical origin, each sample group was grouped based on the frequencies of people movement within the area and the proximity of one area to the others (S2 Fig). Analysis of population dynamics revealed that the relative genetic diversity both in HA and NA of Indonesian sequences had similar patterns regardless of geographical origin, although Balikpapan and Merauke have slightly different patterns.

Rates of nucleotide substitution, TMRCA, and selection pressures of the HA and NA genes

The rates of nucleotide substitution and TMRCA are shown in Table 2, while the positively selected sites and dN/dS values in the selection pressure for HA and NA genes are depicted in Table 3. For HA gene, 44 sites were identified to be positively selected. Specifically, 21 positively selected sites were found in HA1 domain, in which eight sites were located at B-cell epitope (antigenic sites A to E), one site at T-cell epitope, two at the receptor binding sites, and one at the glycosylation site, while 23 positively selected sites were found in HA2 domain, with two sites were located at T-cell epitope (S3 Table). For NA gene, 88 positively selected sites were found (S4 Table), with 12 sites were located at B-cell epitope and four sites at T-cell epitope. No positively selected sites were found at the catalytic sites as well as at the framework sites.

Table 2

Rates of nucleotide substitution and TMRCA of the HA and NA genes.

Gene	Mean rate of nucleotide substitution(x 103 substitution/per/year)		TMRCA(years)
	Mean	95% HPD*	Mean	95% HPD*
HA	3.27	2.9–3.7	5.8	5.4–6.3
NA	4.64	4.0–5.4	5.3	5.3–5.4

*HPD: Highest Probability Density

TMRCA: Time of the most recent ancestor.

Table 3

Positively selected sites and dN/dS values for HA and NA genes.

Gene		Positively selected sites* (n)			dN/dS
		Total	Antigenic sites
			B-cell epitope	T-cell epitope
HA					3.3
	HA1	21	8	1
	HA2	23		2
NA		88	12	4	3.0

*using 0.5 of mean probability as threshold

Discussion

Influenza viruses are unique in that their ecological and evolutionary dynamics occur on the same timescale. Consequently, populations of these fast evolving viruses can accumulate detectable genetic differences in just a few days and can adapt to their environment swiftly [31]. This specific characteristic is examined by investigation of Influenza A/H3N2 virus that continuously changes across different time and location, with a consequence of drift in its phylogeny. To the best of our knowledge, this study is the first in Indonesia that focuses on the evolutionary pattern of the genes encoding the surface glycoproteins influenza A/H3N2 virus across time and geography. In this study, the virus was collected throughout the year with increasing number of samples in rainy seasons (from December to February) when more influenza cases would be more frequent [3]. The ladder-like phylogenies in Figs 3 and 4 represent the antigenic drift of H3N2 as the result of gradual evolution of the virus with multiple mutations occurred cumulatively across years [32]. The central single trunk of each phylogeny characterizes the ancestor of the successful strain with advantageous mutations fixed by natural selection through time [17].

This study observed the presence of positively selected sites which were higher in NA than in HA genes. (Table 3, S3 and S4 Tables). This result was not in accordance with a previous study that reported higher number of positively selected sites in HA compared with those in in NA [33]. However, when analysed more closely, the proportion of positively selected sites found to be in antigenic sites were higher (25%) in HA compared with those observed in NA gene (18%), suggesting that immunological events occurred more frequently in HA than in NA gene. The positively selected sites may vary by the dataset applied, method used, and the significance levels for positive selection determination [18].

Phylogenetic trees both of HA and NA genes revealed that multiple lineages co-circulated and some lineages persisted locally, indicated by the long-line appearance of some strains through the years in the phylogenetic tree. In contrast, the phylogenetic trees also showed that the common ancestor of Indonesian virus existed one or two years before virus sampling (mid-2007). This finding suggests that some lineages became ancestors and were seeded into Northern or Southern Hemisphere in the following seasons or years, and that these viruses might have prevailed and circulated in Indonesian host populations before the dates of collection. Previous study reported that Indonesian H3N2 isolates were detected earlier than the WHO-designated vaccine strain for Northern Hemisphere [34]. This result confirms the long term persistence of virus strains in the tropics, exhibited by co-circulating multiple lineages and often with common ancestors that have prevailed earlier. Consistent with the ‘source-sink’ theory as the evolution model of influenza A virus, the A/H3N2 virus may have persisted locally in the tropics as a reservoir (‘source’ population) before being exported to temperate countries (‘sink’ population) [9].

Phylogenetic analysis also showed a tendency for the viruses isolated from same years to cluster together, regardless geographical origins (Figs 2 and 3); most of the Indonesian samples were intermixed even when the originating islands were separated by oceans. This suggests that these viruses have circulated in various places that are geographically different, with frequent mixing with other viruses from Northern or Southern Hemispheres [35]. The increasing activity of air travel worldwide, with extensive trade and travel between islands or continents connecting millions of people provide a considerable conduit to pathogen transmission. The rapid people movement emphasizes the potential transmission of diseases globally, including the airborne diseases such as influenza [36].

Bayesian skyline plot depicted relatively low genetic diversity within the Indonesian samples through the long and slow seasons in tropical regions [30]. This in contrast with that observed in Northern and Southern Hemispheres with fluctuating levels of relative genetic diversity, which represents stronger natural selection over the short and fast seasons in temperate regions [37]. In relation to neutral evolution, fluctuations in genetic diversity also reflect the underlying changes in the number of infected host [38]. When analysed by geographical origin, a similar pattern of the low genetic diversity of the virus across the Indonesian archipelago was observed, with exception that from Balikpapan and Merauke (S2 Fig). This could be due to the smaller sample size or lower virus population variation in the two places [37].

One limitation of the study was that the collection of samples was carried out during a Ministry of Health geography-based project of three years, which could not enable to study longer period of the antigenic drift. Another limitation was the relatively small number of samples restricted the likelihood of identifying differences between the changes of antigenic sites in each respective year, affecting the ability to interpret significant results. Studies with larger numbers of influenza samples collected over longer time are needed to provide a more distinct evolutionary pattern of Influenza A/H3N2 virus in Indonesia. However, the present study does reveal the evidence of antigenic drift over years and the correspondence of viral strains across geographical regions of the archipelago.

Conclusions

This study describes the evolutionary pattern of the Influenza A/H3N2 virus in Indonesia from 2008 to 2010 based on the analyses of HA and NA genes. Obvious antigenic drift with typical ‘ladder-like’ phylogeny was observed with multiple lineages found in each year, suggesting co-circulation of H3N2 strains at different time periods despite of the various biogeographical regions in Indonesia. This study reconfirms the long term persistence of virus strains in the tropics before being transferred to temperate regions, suggesting the importance to consistently investigate the characteristics of A/H3N2 strains circulating in the tropics, which may contribute important information regarding further spread of the virus. The existence of a particular lineage is most likely the result of adaptation or competitive exclusion among different host populations and combination of stochastic ecological factors, rather than its geographical origin alone.

List of Indonesian samples and the accession number.

(DOCX)

Click here for additional data file.

List of HA and NA sequences from GenBank database included in evolutionary analyses.

(DOCX)

Click here for additional data file.

Positively selected sites within Indonesian HA gene.

(DOCX)

Click here for additional data file.

Positively selected sited within Indonesian NA gene.

(DOCX)

Click here for additional data file.

Schematic diagram of bioinformatics analyses.

(PDF)

Click here for additional data file.

Bayesian skyline plots of HA and NA genes of Indonesian H3N2 virus differentiate according its geographical origin.

Datasets were divided based on the samples origin as follow: Balikpapan (a), Banjarmasin (b), Batam, Aceh and Medan (c), Java and Lampung (d), Jayapura (e), Makasar (f), and Merauke (g). The Bayesian skyline plot at the left panel was genererated from HA gene while at the right panel from NA gene. The x-axis represents the time (mm/yy) and the y-axis represents a measure of relative genetic diveristy (Net, where the Ne is the effective population size and t the generation time from infected host to infected host). The black line represents the mean value while the 95% confidence limits shown in blue line.

(PDF)

Click here for additional data file.

HA MCC tree with sequence identifier.

The Indoneisan sequences were denoted with red labels, while sequences from Southern, Northern Hemispheres, and Tropics were denoted with black labels. Banch color scheme, scale and symbols are similar as those in Fig 2.

(PDF)

Click here for additional data file.

NA MCC tree with sequence identifier.

The Indoneisan sequences were denoted with red labels, while sequences from Southern, Northern Hemispheres, and Tropics were denoted with black labels. Banch color scheme, scale and symbols are similar as those in Fig 3.

(PDF)

Click here for additional data file.