Reconstructing the phylodynamic history and geographic spread of the CRF01_AE-predominant HIV-1 epidemic in the Philippines from PR/RT sequences sampled from 2008 to 2018

Abstract

The Philippines has had a rapidly growing human immunodeficiency virus (HIV) epidemic with a shift in the prevalent subtype from B to CRF01_AE. However, the phylodynamic history of CRF01_AE in the Philippines has yet to be reconstructed. We conducted a descriptive retrospective study reconstructing the history of HIV-1 CRF01_AE transmissions in the Philippines through molecular epidemiology. Partial polymerase sequences (n = 1144) collected between 2008 and 2018 from three island groups were collated from the Research Institute for Tropical Medicine drug resistance genotyping database. Estimation of the time to the most recent common ancestor (tMRCA), effective reproductive number (R_e), effective viral population size (N_e), relative migration rates, and geographic spread of CRF01_AE was performed with BEAST. R_e and N_e were compared between CRF01_AE and B. Most CRF01_AE sequences formed a single clade with a tMRCA of June 1996 [95 per cent highest posterior density (HPD): December 1991, October 1999]. An increasing CRF01_AE N_e was observed from the tMRCA to 2013. The CRF01_AE R_e reached peaks of 2.46 [95 per cent HPD: 1.76, 3.27] in 2007 and 2.52 [95 per cent HPD: 1.83, 3.34] in 2015. A decrease of CRF01_AE R_e occurred in the intervening years of 2007 to 2011, reaching as low as 1.43 [95 per cent HPD: 1.06, 1.90] in 2011, followed by a rebound. The CRF01_AE epidemic most likely started in Luzon and then spread to the other island groups of the country. Both CRF01_AE and Subtype B exhibited similar patterns of R_e fluctuation over time. These results characterize the subtype-specific phylodynamic history of the largest CRF01_AE cluster in the Philippines, which contextualizes and may inform past, present, and future public health measures toward controlling the HIV epidemic in the Philippines.

Keywords: CRF01_AE; HIV-1; Philippines; Subtype B; epidemic; molecular epidemiology; phylodynamics.

Figure 1.

(A) The maximum likelihood phylogenetic tree of Philippine CRF01_AE PR/RT sequences relative to international sequences obtained from the LANL database reconstructed with IQ-TREE 2. The arrow indicates the node, with SH-aLRT and UFBoot branch support values, of the most recent common ancestor of the large MC to which the majority of CRF01_AE sequences belong. This clade is also emphasized with a rectangular highlight. The bottom scale bar shows a reference branch length for 0.02 substitutions/site. Country abbreviations: AT: Austria; AU: Australia; CM: Cameroon; CN: China; DK: Denmark; ES: Spain; FI: Finland; FR: France; GB: United Kingdom; HK: Hong Kong; ID: Indonesia; JP: Japan; KR: South Korea; KW: Kuwait; MM: Myanmar; MY: Malaysia; PH: Philippines; PK: Pakistan; SE: Sweden; SG: Singapore; TH: Thailand; TW: Taiwan; US: USA; VN: Vietnam. (B) Root-to-tip plot generated with TempEst showing a positive correlation between time and divergence or accumulating mutations among sequences, indicating suitability of sequence data for time-scaled phylogenetic and phylodynamic analysis. (C) Posterior probability estimate of the tMRCA of the large Philippine CRF01_AE clade with corresponding MCC tree generated using BEAST2 from the full dataset of 1144 Philippine CRF01_AE sequences.

Figure 2.

Phylodynamics of the largest HIV CRF01_AE cluster in the Philippines between the mid-to-late 1990s and 2018 reconstructed by BEAST2 analysis of the full 1144 sequence dataset of Philippine CRF01_AE PR/RT sequences. (A) Time-scaled MCC tree from BDSKY analysis summarized with TreeAnnotator. Common ancestor internal nodes are highlighted as dots. (B) The change in the effective reproductive number (R_e) over time, with the median represented as a solid line and the 95 per cent HPD as a shaded interval, obtained from BDSKY analysis. A dashed line indicates the value of R_e equal to 1. (C) The change in CRF01_AE effective population size (N_e) over time in log scale, with the median represented as a solid line and the 95 per cent HPD as a shaded interval, obtained from BSP analysis.

Figure 3.

Analysis of geographic spread and relative migration rates of HIV CRF01_AE across Philippine island groups. (A) MCC tree of the largest Philippine CRF01_AE cluster from the full dataset of 1042 sequences with sample collection date and location metadata generated with BEAST under the ‘phylogeo’ model and summarized with TreeAnnotator. All nodes achieved a posterior probability support value of 1.0 except for the nodes labeled with a star symbol, corresponding to a support value greater than or equal to 0.95 and less than 1.0, and a circle symbol, corresponding to a support value greater than 0.66 and less than 0.95. The branch color is set according to the most likely geographic location of a branch at the level of Philippine island groups: Luzon, Visayas, Mindanao. (B) The forest plot with median and 95 per cent HPD estimate of the asymmetric relative migration rates of CRF01_AE between all pairs of Philippine island groups. A dashed vertical line indicates a relative migration rate equal to 1.0, or no greater or less than the mean of migration rates. (C) The phylogeographic spread of CRF01_AE between Luzon, Visayas, and Mindanao island groups over time visualized with SPREAD4. Edges represent transitions between island groups. The size of the polygons over island groups corresponds to the intensity of localized virus transmission at the specified location and time.

Figure 4.

Comparison of phylodynamics between Philippine CRF01_AE and Subtype B. (A) Time-scaled MCC tree from BDSKY analysis of the full set of 192 Subtype B PR/RT sequences summarized with TreeAnnotator. Common ancestor nodes are highlighted as dots. (B) The change in the Subtype B effective reproductive number (R_e) over time obtained from BDSKY analysis, superimposed over that of CRF01_AE, with the medians represented as solid lines and the 95 per cent HPDs as shaded intervals. A dashed line indicates the value of R_e equal to 1. (C) The change in Subtype B effective population size (N_e) over time in log scale obtained from BSP analysis, superimposed over that of CRF01_AE, with the medians represented as solid lines and the 95 per cent HPDs as shaded intervals.