Low Postseroconversion CD4⁺ T-cell Level Is Associated with Faster Disease Progression and Higher Viral Evolutionary Rate in HIV-2 Infection

The HIV-2 env V1-C3 region was successfully reverse transcribed and amplified from 53 plasma samples that had been collected longitudinally from 16 study participants of a well-described cohort of police officers in Guinea-Bissau (22, 23). Inclusion criteria and clinical characteristics are presented in Materials and Methods (see also Table S1 and S2 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). The median observation time from inclusion until the last registered visit for the 16 included individuals was 19.2 years (interquartile range [IQR], 15.0 to 20.8 years). A median of seven clones from each time point were sequenced, and the median time between the collection dates of the first and last amplified patient-specific samples was 7.9 years (IQR, 5.2 to 14.0 years). Among the 528 sequences, 119 putative recombinant sequences were removed, leaving 409 sequences for evolutionary analyses. Phylogenetic subtype analysis showed that all sequences belonged to HIV-2 group A. Moreover, all sequences from each individual formed distinctive monophyletic clades in the phylogeny, indicating that sequences from the study participants were not subjected to superinfection, coinfection, mix-up during sample handling, or contamination during the laboratory procedures (see Fig. S1 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE).

To assess the link between HIV-2 evolutionary estimates and disease progression, we classified the study participants into groups of faster disease progressors and slower disease progressors based on longitudinal CD4⁺ T-cell dynamics from the entire HIV-2-infected population of the Guinea-Bissau police cohort. To determine if faster disease progression can be differentiated from slower disease progression by CD4⁺ T-cell decline and level in cases of HIV-2 infection also, as had been previously suggested for HIV-1 infection (12, 14), we analyzed the time to the appearance of AIDS from the first recorded percentages of CD4 (CD4%) for all HIV-2-infected individuals in the cohort with two or more CD4% measurements (n = 192). Three stratifications were used, and faster and slower progressors were defined as those patients whose CD4% values were above and below the mean of the values determined for all participants for each stratification, respectively (Fig. 1). The first stratification, referred to as the CD4% decline rate, was based on individual coefficients of regression for CD4%. Eighty-one participants were classified as faster progressors (mean CD4% decline, 2.7% per year [standard deviation {SD}, 2.6]) and 111 as slower progressors (mean CD4% increase, 0.5% per year [SD, 1.6]). The median time to AIDS from the first recorded CD4% value was 11.7 years (95% confidence interval [CI], 7.3 to 16.1 years) for faster progressors with a faster CD4% decline and 16.8 years (CI, 12.3 to 21.3 years) for slower progressors with a slower CD4% decline (P = 0.008 [log rank test]) (Fig. 1A). The second stratification, referred to here as the CD4% level, was determined as the CD4% level at the midpoint in time between the first and last recorded CD4% levels using the regression line generated in the first stratification. Eighty-seven participants were classified as faster progressors (mean CD4% level, 21.0% [SD, 4.3%]) and 105 as slower progressors (mean CD4% level, 35.0% [SD, 5.7%]). The median time to AIDS from the first recorded CD4% value was 9.4 years (95% CI, 6.7 to 12.1 years) for faster progressors with a low CD4% level and 15.5 years (95% CI, 14.3 to 16.6 years) in slower progressors with a high CD4% level (P < 0.001 [log rank test]) (Fig. 1B). In the third stratification, referred to as the combined coefficient, the CD4% decline rate and CD4% level were combined to have equal influences on the stratification (see Materials and Methods). Eighty-five participants were classified as faster progressors (combined coefficient below the mean) and 107 as slower progressors (combined coefficient above the mean). The median time to AIDS from the first recorded CD4% value was 8.6 years (95% CI, 6.5 to 10.8 years) for faster progressors with a low combined coefficient and 18.7 years (95% CI, 13.6 to 23.8 years) for slower progressors with a high combined coefficient (P < 0.001 [log rank test]) (Fig. 1C). Interestingly, only 50% of the participants were consistently classified as faster or slower progressors in the stratifications (see Table S3 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). A Cox proportional-hazards model with CD4% decline and level as discrete covariates (according to the stratifications described above) showed that CD4% decline and level were independently associated with disease progression rate (P = 0.007 for CD4% decline and P < 0.001 for CD4% level [Wald test]; the independence of these covariates was determined by a nonsignificant interaction term between the covariates corresponding to P = 0.85 [Wald test]).

Since all three stratifications described above were significantly linked to the rate of disease progression, we decided to assess differences in virus evolutionary parameters for all three stratifications (Table 1). Model evaluation within the hierarchical phylogenetic modeling (HPM) framework indicated that a strict clock with a constant size demographic model generally fitted our data best for both the nucleotide and codon models (see Materials and Methods and Text S1 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE for details). The overall mean evolutionary rate of the V1-C3 region was 23.5 × 10⁻³ codon substitutions/site/year (95% highest posterior density [HPD] interval, 20.3 × 10⁻³ to 26.6 × 10⁻³) (Table 2). No association was found between evolutionary rate and disease progression when individuals were classified as faster or slower progressors based on CD4% decline rate (Bayes factor [BF] = 0.3) (Table 2). In contrast, when individuals were classified based on either CD4% level or the combined coefficient, the mean evolutionary rate was significantly higher in faster progressors (28.6 × 10⁻³ codon substitutions/site/year; 95% HPD, 24.2 × 10⁻³ to 33.5 × 10⁻³) than in slower progressors (14.9 × 10⁻³ codon substitutions/site/year; 95% HPD, 12.2 × 10⁻³ to 17.6 × 10⁻³) (BF = 20.3) (Table 2). Similar differences were found for the nucleotide substitution rate (BF = 0.3 for association with CD4% decline rate and BF = 17.5 for association with CD4% level or the combined coefficient) (see Table S4 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE).

To determine if the evolutionary rates were similar in different regions of env, we partitioned the data set into four well-defined regions: V1V2, C2, V3, and C3. The evolutionary rate was higher in both the V1V2 and C3 regions than in the C2 and V3 regions, respectively (Table 2) (P < 0.05 [Wilcoxon signed rank test with Bonferroni correction for all pairwise comparisons]). No significant difference was found between faster and slower progressors in any of the four V1-C3 regions when individuals were stratified by CD4% decline rate (BF values ranged from 0.2 to 0.3). In contrast, moderate to strong associations between evolutionary rate and disease progression were found in the V1V2 (BF = 11.8), C2 (BF = 28.4), and C3 (BF = 6.1) regions for individuals grouped according to CD4% level or the combined coefficient. A weak association was also noted in the V3 region (BF = 2.4).

Since the CD4% level and the combined coefficient stratifications resulted in identical groupings (Table 1) and all analyses indicated a strong association between evolutionary rate and CD4% level—in contrast to the absence of associations between evolutionary rate and CD4% decline rate—only results based on the CD4% level stratification are presented for the subsequent analyses.

Next, we hypothesized that the differences in evolutionary rate between the progressor groups could have originated in parameters influencing selection pressure or viral replication. We therefore estimated the ratio of nonsynonymous substitution rates to synonymous substitution rates (the dN/dS rate ratio), which indicates whether a gene or site had been subjected to positive selection (dN/dS > 1) or to negative selection (dN/dS < 1) or had evolved neutrally (dN/dS = 1). We found global negative selection over the entire V1-C3 region (dN/dS rate ratio = 0.56; 95% HPD interval, 0.49 to 0.63), with no statistical differences apparent between faster and slower progressors (BF = 1.3) (see Table S5 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). The region-specific analyses showed that the V1V2 and C3 regions evolved neutrally (dN/dS rate ratios, 1.27 [95% HPD interval, 0.84 to 1.80] and 0.94 [95% HPD interval, 0.70 to 1.22], respectively) whereas the C2 and V3 regions had been subjected to strong negative selection (dN/dS rate ratios, 0.22 [95% HPD interval, 0.18 to 0.27] and 0.41 [95% HPD interval, 0.27 to 0.57], respectively). No significant differences were found in dN/dS rate ratios between the progressor groups in the V1V2, V3, or C3 regions (BF values of <3 for all regions), whereas a moderately higher level of negative selection was found in the C2 region among faster progressors than among the slower progressors (BF = 3.7) (see Table S5 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE).

Whereas the dN/dS rate ratio provides an estimate of the overall selective pressure, analyses of expected nonsynonymous (E[N]) and expected synonymous (E[S]) substitution rates can add detailed information concerning replication and selection pressure (17). The E[N] and E[S] divergences were estimated and plotted as accumulated divergence rates over time from the first analyzed sample (Fig. 2; see also Table S6 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). The HIV-2 divergence rates increased in a relatively linear manner during the asymptomatic phase of infection, and viral populations in faster progressors accumulated both nonsynonymous and synonymous substitutions faster than those in slower progressors. The median E[N] rate ranged from 3.0 × 10⁻³ to 6.6 × 10⁻³ nucleotide [nt] substitutions/site/year for faster progressors and from 1.4 × 10⁻³ to 4.6 × 10⁻³ nt substitutions/site/year for slower progressors (P = 0.005 [two-tailed Mann-Whitney U test {M-W}]). The median E[S] rate ranged from 4.5 × 10⁻³ to 11.2 × 10⁻³ nt substitutions/site/year for faster progressors and from 1.1 × 10⁻³ to 4.0 × 10⁻³ nt substitutions/site/year for slower progressors (P < 0.001 [M-W]). These results are in line with the higher overall viral evolutionary rates for faster progressors than for slower progressors seen in the HPM analyses and clarify why no statistically significant differences were found in dN/dS rate ratios.

Although the analysis described above indicated the presence of general negative selection across the entire V1-C3 region, we hypothesized that a number of specific residues could have been subject to positive selection. To assess this hypothesis, we used a Renaissance counting procedure to estimate the dN/dS rate ratio at each codon site (Fig. 3). Only a low proportion of sites showed signatures of positive selection, and the proportions of positively selected sites were not significantly different between the viral populations of faster and slower progressors across the V1-C3 region (7% and 11%, respectively; P = 0.503 [M-W]) or in the V1-C3 regions in separate analyses (see Table S7 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). While the numbers of sites under positive selection were similar in the faster and slower progressors in the V1V2, V3, and C3 regions, slower progressors had more sites under positive selection in the C2 region than faster progressors (35 versus 20 sites) (P = 0.026 [two-tailed Fisher’s exact test {FET}]) (Fig. 3). Since it has been shown that the HIV-2 C2 and C3 regions are exposed similarly to the corresponding HIV-1 regions and are under negative selection, we hypothesized that the effects of positive selection on amino acids subjected to structural and functional constraints would negatively impact viral fitness (24). We therefore defined the amino acids critical to viral fitness as those amino acids that were conserved between HIV-2 and the simian immunodeficiency virus (SIV) from which HIV-2 originated through cross-species transmission (i.e., SIVsm, the SIV infecting sooty mangabeys [Cercocebus atys]) and compared the ratios of positive selection of viral lineages between slower and faster progressors (25). Among the 246 amino acids in the V1-C3 regions, 84 were conserved between SIVsm and HIV-2 (see Table S8 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). Comparisons of the numbers of conserved HIV-2/SIVsm sites under positive selection showed that slower progressors had more conserved sites under positive selection than faster progressors (20 versus 5 sites (P = 0.002 [FET]) in the V1-C3 region. The differences between the slower and faster progressors with respect to positive selection at conserved sites were most highly pronounced in the C2 region (12 versus 3 sites, respectively) (P = 0.021 [FET]), while such differences were not observed in the V1V2, V3, and C3 regions. To confirm previous observations that C2 is well exposed on the HIV-2 env gene, we used the published structural data of HIV-2 gp125 to visualize the amino acids in the V1-C3 region (24, 26). This analysis indicated that the majority (15 of 22) amino acids associated with positive selection mapped to exposed surfaces on HIV-2 gp125 (see Fig. S2 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). Moreover, the surface accessibility of amino acids within the V1-C3 regions showed that the positively selected sites were associated with residues with solvent-exposed surfaces in the C2 region (29 of 68 amino acids) (P = 0.040 [FET]).

This study included individuals from a large cohort of police officers in Guinea-Bissau, West Africa, which was formed in 1990 (22, 23). At inclusion, and at follow-up visits scheduled with an interval of 12 to 18 months, individuals were examined and a plasma sample was collected. The civil war in 1998 to 1999 temporarily (from June 1998 until the end of 2002) ended inclusion, but annual visits of previously included individuals were resumed in July 2000. The cohort was followed routinely until February 2011, when the cohort was closed. In September 2013, however, selected individuals from the cohort were asked to participate in a special sampling round, including a clinical examination and collection of a plasma sample. In early 2006, the police cohort was included in the national antiretroviral therapy (ART) program which was introduced into Guinea-Bissau in 2005. HIV testing was performed at the National Public Health Laboratory (LNSP), Bissau, as previously described (22). The CD4⁺ T-cell count and CD4% were determined at each follow-up visit after infection (14, 22, 23). In the absence of commercial HIV-2 RNA assays, HIV-2-infected patients are generally monitored by CD4⁺ T-cell levels, particularly in areas of endemicity. Consequently, viral load measurements have not been included as a standard procedure in Guinea-Bissau, preventing us from a comprehensive analysis of viral load data in this study.

The cohort includes 438 seroincident and seroprevalent HIV-2-infected individuals, 83 of whom had an estimated date of seroconversion, defined as the midpoint between the last HIV-2 seronegative sample collection date and the first seropositive sample collection date. Individuals with long follow-up series, including both CD4% measurements and available plasma samples, were considered for inclusion in the study. Our initial goal was to conduct this study by including only individuals with estimated dates of infection. Thus, amplification of viral RNA was attempted on plasma samples from seroincident individuals from whom three or more longitudinal plasma samples were available (n = 19). As expected, due to the generally low viral loads in HIV-2 infections, amplification was successful for only a minority of samples and seven individuals with two or more successfully amplified longitudinal samples could be included in the study. We therefore decided to also include seroprevalent individuals (i.e., individuals who were HIV-2 infected already at enrollment) in the study. Due to the anticipated difficulties in amplifying HIV-2 RNA, only individuals with more than six available plasma samples were considered for inclusion. Amplification was attempted on samples from 19 individuals, and successful amplification of two or more longitudinal samples was achieved for 9 individuals. Taking the data together, amplification was attempted on samples from 38 individuals, where 16 individuals fulfilled the inclusion criterion of the availability of two or more amplified longitudinal samples (total, 53 samples). The majority of the plasma samples were collected from the 16 participants included in this study during the asymptomatic phase of infection, defined by a CD4⁺ T-cell count of >200 cells/μl, CD4% of >14, and a lack of clinical AIDS symptoms (WHO stage 4 and CDC stage C [47, 48]). However, seven samples from six individuals were collected after the individuals had developed AIDS, and two samples from two individuals were collected after the individuals had initiated antiretroviral therapy (see Table S2 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). Baseline characteristics of the 16 individuals are presented in Table S1 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE, and successfully amplified samples are listed in Table S2 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE.

The 16 included individuals were subsequently stratified as faster or slower progressors using three different parameters: CD4% decline rate, CD4% level, and a combined coefficient. The combined coefficient was achieved by combining the CD4% decline rate and the CD4% level. Since these values were quantified on different scales and since the CD4% decline rate could be either positive or negative, we transformed the values before combining them. The two variables were transformed and rescaled to have equal levels of influence on the combined coefficient. More specifically, to account for negative CD4% decline rates, all rates were transformed to positive rates with the corresponding relative difference. Next, CD4% decline rates and CD4% levels were rescaled to include the same extrema. Individual combined coefficients were then determined for each study participants by multiplying the rescaled CD4% decline rates and CD4% levels. There were no differences in age between the groups at the time of the first HIV-2-positive sample (data not shown).

Viral RNA was extracted from patient plasma samples using an miRNeasy microkit (Qiagen, Stockholm, Sweden) with minor modifications to the manufacturer’s instructions. Briefly, 200 μl of plasma was disrupted in 2,000 μl QIAzol and loaded onto an RNeasy MinElute Spin column in the presence of 15 μg carrier RNA (Qiagen, Stockholm, Sweden). DNA was removed using an on-column DNase treatment (Qiagen, Stockholm, Sweden), and purified RNA was eluted in 22 μl RNase-free H₂O. An approximately 935-bp fragment that included the complete V1-C3 region of env (nt 6986 to 7920 in the BEN reference sequence]; GenBank accession number M30502) was amplified using 9.5 μl eluted RNA in a SuperScript III One-Step reverse transcription-PCR (RT-PCR) system with Platinum Taq DNA polymerase followed by a seminested PCR approach using Platinum Taq High Fidelity (Invitrogen, Copenhagen, Denmark). Primers KH2_OF (5′-GAGACATCAATAAAACCATGTGTC-3′) and TH2_OR (5′-TTCTGCCACCTCTGCACTAAAGG-3′) were used for One-Step PCR, and primers KH2_OF and KH2_OR (5′-ACCCAATTGAGGAACCAAGTCA-3′) were used for nested PCR (5, 42). Following the initial cDNA synthesis performed for 30 min at 50°C, the PCR conditions were identical for One-Step PCR and the nested PCR: initial denaturation for 2 min at 94°C; 40 cycles of 15 s at 94°C, 30 s at 50°C, and 1 min at 68°C; and a final elongation step for 5 min at 68°C. The sensitivity of the PCR was found to be 12 RNA copies/ml plasma as determined by a dilution series performed with an electron microscopy-counted HIV-2 particle (Advanced Biotechnologies, Eldersburg, MD, USA). Molecular cloning of the amplified fragments using a pCR2.1 TOPO cloning system (Invitrogen) was performed by BaseClear BV (Leiden, The Netherlands), and 12 individual clones were routinely picked for subsequent sequencing on both strands using conventional M13 primers.

Sequences were manually edited using CodonCode Aligner v1.5.2 (CodonCode Corporation, Dedham, MA, USA) and aligned in MEGA5 using the Clustal algorithm (49). Sequences of poor quality and sequences containing stop codons were removed from the analysis. Parts of the sequence that were difficult to align were removed in full codons to preserve an open reading frame, resulting in a final alignment length of 774 bp. The alignment spanned the last 30 bp in the 3′ end of the C1 region, the entire V1-C3 region, and the first 6 bp in the 5′ end of the V4 region. To analyze different subregions, the fragment was partitioned as follows: V1V2, bp positions 31 to 201 (positions 7040 to 7318 in the BEN reference sequence); C2, positions 202 to 498 (7319 to 7615); V3, positions 499 to 600 (7616 to 7717); C3, positions 601 to 768 (7718 to 7891).

Kaplan-Meier analyses were performed for progression time to AIDS. Cases that did not reach AIDS during follow-up were right censored at their last clinical examination date. Statistical comparisons were performed by the log rank test. A Cox proportional-hazards model was applied to determine independency of covariates (as defined by nonsignificant interaction terms between covariates).

Putative intrapatient recombinant sequences were identified by the pairwise homoplasy index (PHI) test using an exhaustive and iterative search algorithm and were then removed from the data set (50) (the Perl script for the iterative search is available from the authors upon request). Maximum likelihood (ML) phylogenetic trees were reconstructed using the inferred model, GTR +I + G, with Garli v2.0 (51). Statistical support for internal branches was determined by ML-based approximate likelihood ratio test (aLRT) Shimodaira-Hasegawa (SH)-like branch support, as implemented in PhyML 3.0 (52). SH values of >0.9 were considered statistically significant (53). For subtype analysis, our data set was aligned with reference sequences of the major HIV-2 subtypes (downloaded from the Los Alamos Sequence Database [54]) in MEGA5 using the Clustal algorithm (49), followed by phylogenetic analysis.

Analyses of HIV-2 evolutionary rates were performed in BEAST v1.7.5 (55) by reconstructing Bayesian rooted and time-measured phylogenetic trees. Unless otherwise stated, all analyses were performed by running a Markov chain Monte Carlo (MCMC) analysis for 50 × 10⁶ generations, with sampling performed every 2,500 to 5,000 generations. Convergence was determined by calculation of effective sample sizes (ESS) of >100 and inspection of traces, as assessed in Tracer v1.6 (available from http://beast.bio.ed.ac.uk/software/tracer/), following removal of 10% after burn-in. In explorative analyses, we evaluated different models to find the one that best fit our data (see Table S9 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE). For each individual, the nucleotide substitution rates were estimated using two different clock models (strict and uncorrelated lognormal relaxed clocks), two different demographic models (Bayesian skyline plot and constant size), and either partitioned (1st plus 2nd and 3rd) or nonpartitioned codon positions. Nucleotide substitution rates were estimated using the Hasegawa, Kishino, and Yano (HKY) substitution model (56) with gamma-distributed rates. Exploratory analyses were performed by running a single analysis of 50 × 10⁶ MCMC generations as described above. All subsequent analysis were performed in duplicate, and the results were combined in LogCombiner v1.7.5 after the removal of 10% after burn-in (55). Phylogenetic trees were visualized in Figtree v1.40 (available from http://tree.bio.ed.ac.uk/software/figtree/).

Hierarchical phylogenetic modeling (HPM) with fixed effects, as implemented in BEAST v1.7.5, was used to compare differences in evolutionary rate between progressor groups (55, 57). HPM allows simultaneous analysis of sequence data from multiple individuals. Information of evolutionary parameters is pooled across populations or individuals through hierarchical prior specification, resulting in a shrinkage effect of the variation when data are sparse for study participants. In addition, fixed effects across a group of individuals (in this case, the faster and slower progressor groups) can be included to test and quantify differences between them. Analyses were performed using a strict or uncorrelated lognormal relaxed clock model with a constant population size model as the tree prior. Both nucleotide substitution rates (HKY) and codon substitution rates (determined using the GY94 codon model) were estimated (56, 58). In a region-specific analysis, the data set was partitioned into the env V1V2, C2, V3 and C3 regions and the evolutionary rate was estimated for all regions simultaneously by HPM (strict clock, constant size, and GY94 codon model). Evolutionary rates were compared between groups as assessed by Bayes factors (BFs). A BF value of >3 was considered to represent a significant association (59).

The ratio of nonsynonymous and synonymous substitution rates (dN/dS rate ratio) has been widely used as an indicator of selection and molecular adaptation. However, the dN/dS rate ratio cannot be used to detect simultaneous increases or simultaneous decreases in nonsynonymous and synonymous rates (60). To address this issue and to further dissect and explore the molecular adaptation process, the rate of substitution in every branch in a tree can be divided into expected nonsynonymous (E[N]) and expected synonymous (E[S]) substitution rates. These absolute rate estimates reflect the respective contributions of E[N] and E[S] substitution rates to the overall substitution rate for a particular branch (17). These rates are uncorrected for the number of possible nonsynonymous and synonymous alterations (i.e., the number of possible nonsynonymous alterations is higher than the number of possible synonymous alterations). Consequently, direct comparisons between E[N] and E[S] estimates may be difficult to interpret. However, relative differences between, e.g., patient groups in E[N] or E[S] estimates can still be explored. Analyses of E[N] and E[S] and divergence plots were performed as described by Lemey et al. (17). Briefly, 200 random trees from the HPM analysis (the model settings included HPM performed with a relaxed clock, constant population and constant size, and the GY94 codon model [58]) were selected for each individual to determine E[N] and E[S] rates in HyPhy 2.2.0 (61). The accumulated divergence was estimated and plotted over time using the weighted average rate. The analysis was performed using a relaxed clock model because a strict clock model, by definition, would assume a linear relationship.

A Renaissance counting approach, as implemented in BEAST v1.8.1(55), was employed to estimate the ratio of nonsynonymous and synonymous substitutions at each codon site, allowing identification of sites that were under positive or negative selection (62). A dN/dS value of 1 suggests neutral selection, a dN/dS value of >1 suggests positive selection, and a dN/dS value of <1 suggests negative selection. The analysis was performed using a strict clock model with a constant population size model as the tree prior and a triply partitioned (first, second, and third) HKY nucleotide substitution model (56). Significant positive or negative selection was considered to have occurred when the HPD interval did not include a value of 1.

The env reference sequences of HIV-2 and SIVsm (n = 61 and n = 54, respectively) were downloaded from the Los Alamos Database and aligned using the Clustal algorithm. A strict consensus sequence was obtained using Consensus Maker (available at https://www.hiv.lanl.gov/). Amino acid positions with 100% amino acid conservation were considered to represent HIV-2/SIVsm conserved positions (see Table S8 at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE).

The solvent accessibility surface area of the HIV-2 gp125 core (PDB accession number 5CAY) was quantified using the Web server running DSSP (http://www.cmbi.ru.nl/xssp/) (63). The relative solvent accessibility (RSA) was calculated as the accessible area divided by the maximum accessible area of the amino acid in extended tripeptide confirmation (Gly-X-Gly), as previously described (64). Residues with an RSA value of ≥5% were considered to be surface exposed (65). The HIV-2 Env protein (PDB accession number 5CAY) was visualized and surface exposed, and positively selected amino acids were edited using UCSF Chimera (https://www.cgl.ucsf.edu/chimera/).

Statistical analyses were performed in IBM SPSS Statistics 21 using the two-tailed Mann-Whitney U test (M-W) for continuous variables, the two-tailed Fisher’s exact test (FET) for categorical variables, Spearman’s rho for correlation analyses, and Friedman’s test for comparisons of multiple groups. Results from the Friedman’s test showing statistical significance were then analyzed by the Bonferroni-corrected Wilcoxon signed rank test for pairwise comparisons of groups.

The study was approved by the ethics committees of the government of Guinea-Bissau; the University of Lund, Sweden; and the Karolinska Institute, Stockholm. Study participants were counselled and provided informed oral consent.

Supplementary information (including supplementary text, figures, and tables) has been deposited at https://dochub.com/patrik-medstrand/zqp8z8/supplemntaryinformation_palm_mbio_ver2_180822?dt=Nx2HKD__sa_SKWC-XzVE.

The data sets generated and analyzed during the current study are available in the GenBank repository (accession numbers: KM390990 to KM391398).