The CHARTER Study is a prospective, observational study of neurocognitive outcomes in PLWH. Ambulatory, PLWH were enrolled at six medical centers in the U.S., as described previously (Heaton et al., 2010). Detailed, structured interviews and comprehensive neurocognitive examinations, as well as laboratory assessments were conducted to collect information on HIV disease and treatment-related factors from participants at baseline and 6-months follow-up visits according to a protocol that was standardized across sites. Details of CHARTER study eligibility and assessment protocols have been published before (Heaton et al., 2010). For participants who consented, CSF samples were also obtained by lumbar puncture. To assess the GDS, participants underwent a comprehensive test battery that involved seven neurocognitive domains and were assigned test scores which were then converted to demographically corrected standard scores (T-scores) (Heaton et al., 2010). A single T-score was calculated for each of the cognitive domains by averaging the T-scores for each of the tests in that domain. The domain impairment is determined when an individual’s T-score is below one standard deviation from the mean (Antinori et al., 2007). A composite GDS was then derived, as a continuous measurement, with deficit scores converted from standard domain T-scores, using a published objective algorithm (Carey et al., 2004; Heaton et al., 2010). Participants’ NCI status was also determined by applying a GDS cutoff of 0.5, with ‘Normal/Not impaired’ defined by a GDS < 0.5 and ‘Impaired’ by a GDS ⩾ 0.5 (Carey et al., 2004; Jia et al., 2017). The GDS incorporated adjustments for practice (or learning) effects to account for prior neurocognitive testing. Neuro-relevant comorbid conditions were evaluated by experienced clinicians. Conditions such as developmental learning disability and major head injury with loss of consciousness were excluded and the rest were categorized as either “incidental” (absent, minimal and non-contributory) or “contributing” (mild-to-moderate) to NCI (Salvatore et al., 2009; Heaton et al., 2010).

Plasma HIV RNA (viral load) was determined by reverse transcriptase PCR (Heaton et al., 2010). The CD4+ nadir was obtained by self-report and confirmed by documented prior measurements in a subset of CHARTER Study participants (Ellis et al., 2011). CSF HP protein levels were quantified in 405 participants using multiplex bead-based suspension array immunoassays (Kallianpur et al., 2019). Detailed methods for quantification of CSF HP and other markers were published previously (Heaton et al., 2010; Ellis et al., 2011; Kallianpur et al., 2019). Genomic DNA was extracted from peripheral blood mononuclear cells collected at the baseline CHARTER visit using PUREGENE (GentraSystems, Inc., Minneapolis, MN). Genotyping was conducted using the Affymetrix Genome-Wide Human SNP Array 6.0 ^TM by the Vanderbilt Technologies for Advanced Genomics (VANTAGE) at Vanderbilt University in two batches: n = 576 samples were genotyped before 2009 and n = 506 (six repeated for QC) were genotyped in 2012, due to funding reasons only (Jia et al., 2017). Nevertheless, only minor changes were observed from explicit testing of batch effects (Jia et al., 2017). Due to limited accuracy for imputing APOE, the APOE genotypes were determined for a subset (n = 401) of CHARTER participants by genotyping of rs7412 and rs429358, using TaqMan predesigned SNP genotyping assays (C_904973_10 and C_30846793_20; Applied Biosystems, Foster City, CA) as described in a previous publication (Morgan et al., 2013).

The QC and basic data cleaning pipeline used for CHARTER genomic data was published previously (Jia et al., 2017); here we describe the additional QC work that was conducted for HP imputation. Since the HP reference panels are separated for AFR and EUR, we performed all the additional QC steps respectively for AFR and EUR participants. The overall genotyping call rate was checked for each study participant and we found the overall call rate was > 95 % for all samples. The whole dataset was then separated into the AFR set and the EUR set, according to ethnicity clusters defined by principal component (PC) ancestry clustering, which, as previously performed by the CHARTER study groups (Samuels et al., 2016; Jia et al., 2017), is an approach of assigning individuals’ ancestry groups by clustering of the continental-ancestry-corresponding PCs obtained from their genomic composition. In each ancestry group, for each SNP/marker, the overall call rate was then checked, respectively. Markers with < 97 % call rate (i.e., ⩾ 3 % missing rate) within each population were pruned from the dataset. Finally, SNPs that failed the Hardy-Weinberg Equilibrium test (with p-value ⩾ 0.001 ) within each population were removed. The entire QC process was conducted using PLINKv1.9 (Chang et al., 2015; Purcell and Chang, 2019).

We adopted published HP CNV imputation references for AFR and EUR (Boettger et al., 2016). An HP marker that was collapsed from all four HP subtype markers was added to the imputation reference with “0” represents HP1 and “1” represents HP2. We performed in silico validation and found the HP genotype we obtained from the HP marker 100% identical to the genotype we obtained from subtype markers. The full lengths of chromosome 16 of the individuals in the HP imputation reference were extracted from 1,000 Genomes (1KG) data and pre-phased using SHAPEITv2 (Delaneau et al., 2013) software. The HP imputation markers were then extracted from the pre-phased chromosome 16 to obtain the phased HP imputation reference panels. Further validation of this imputation strategy and optimization of the IMPUTEv2 (Howie et al., 2009) software settings were conducted using: 1. HP-genotype-removed imputation reference as data input with an accuracy metric and 2. the GTEx sequencing and expression data by comparing the imputed HP genotype and the read count of the exon 4 & 5 junction (unique to HP2) from RNA sequencing. The HP region of all CHARTER Study participants was imputed for HP genotypes, using both of the phased references and samples that were extracted based on the reference population; i.e., EUR participants were extracted from the imputation using European reference and AFR participants were extracted from the imputation using African reference. The HP imputation was conducted using the IMPUTEv2 software.

The imputed dosages were hardcalled using a 0.9 threshold; in other words, dosages ⩾ 0.9 were converted to genotypes and dosages < 0.9 were replaced by ‘Null’ and excluded from hardcall analyses. The hardcalling and data preparation steps were conducted using Python3.7 programming language. Due to variations in HP allele frequency and other genetic complexities, analyses were conducted respectively for AFR and EUR. For continuous GDS as an outcome, analysis of variance (ANOVA) was applied, whereas for GDS defined NCI as the outcome, the χ ²-square test was used. Multivariate linear regression and logistic regression were used for continuous and categorical outcomes, respectively, with adjustments of age (continuous), CD4+ nadir (continuous), plasma HIV RNA (continuous), sex (categorical), comorbid condition (contributing vs. incidental), and CSF HP protein levels (continuous) if specified. Combined analyses of both populations were adjusted additionally for the first 3 PCs. Longitudinal analyses were conducted on repeated measurements of GDS, NCI, cognitive domain T-scores, and domain impairments with generalized estimating equations (GEE) method using age as the time variable. To obtain a more detailed estimation, we first estimated the ages of the records with a minimum increment of 0.5, and then fit into a GEE model with an empirical estimator and first-order auto-regressive (AR1) covariance structure to test if the HP CNV impacts the trajectory of changes of the outcomes. An identity link was used for the continuous outcomes and a logit link was used for the dichotomous outcomes. Moderation effects between variables were tested using interaction terms in models. In other words, a significant A-B interaction effect suggests that the effect of A on the outcome depends on B and vice versa. Additive and dominant genetic effects were tested for genetic variables by applying different coding methods. All hardcall analyses were performed using R statistical language. Specifically, the geepack (Yan, 2002; Yan and Fine, 2004; Halekoh et al., 2006) R library was used for longitudinal data analyses. SNPTEST (Marchini et al., 2007) software was used for frequentist analyses on continuous GDS and NCI using directly imputed dosages/certainties.

Statistics of all the study variables are in Table 1. Due to the complex evolutionary history of HP alleles in different populations, all analyses were conducted in AFR and EUR ancestral groups separately as well as jointly. The HP imputation process was validated by comparing the imputed HP genotypes and exon 4 & 5 junction counts from the Genotype-Tissue Expression (GTEx) RNA sequencing data (Supplementary Figure S1). HP genotypes were then imputed in the CHARTER AFR and EUR participants with high certainty (IMPUTEv2 info metric = 0.831 for AFR, and 0.830 for EUR). A low imputation dosage indicates a low certainty of the imputed genotype, thus, we filtered out the low confidence genotypes (hardcalls) with a stringent threshold, 0.9, to ensure that we had high quality genotypes for further analyses (Supplementary Figure S2). We obtained 371 hardcalls (86% of 432 samples) and 395 hardcalls (80% of 491 samples) in EUR and AFR, respectively. After hardcalling, Hardy-Weinberg Equilibrium was tested in both populations and no significant deviations were found.

Associations between the HP CNV and CSF HP protein levels were determined using an ANOVA F-test among the participants who had CSF HP levels measured (n = 283, after hardcalling) (Kallianpur et al., 2019). In both populations, the presence of additional HP2 alleles decreased the CSF HP protein levels, as shown in Figure 1A, (p = 2.17 × 10^–5 in n = 131 AFR participants, p = 1.70 × 10^–9 in n = 146 EUR participants, and p = 7.85 × 10^–12 when combined).

The distributions of baseline GDS for each HP genotype in AFR and EUR populations are shown in Figure 1B; we observed an almost linear relationship between the mean GDS of different HP genotypes and the number of HP2 alleles in both AFR and EUR. However, neither ANOVA nor linear regression showed a statistically significant association between the HP CNV and the GDS. We then looked at the NCI as a dichotomous trait, defined by GDS ⩾ 0.5 as previously described (Heaton et al., 2010). χ ²-square tests showed a significant relationship between HP CNV and NCI in both populations, with p = 0.027 in AFR and p = 0.013 in EUR.

We adjusted for known risk factors for the dichotomous NCI in PLWH, including age, sex, comorbidity conditions, plasma HIV RNA (viral load), and CD4+ nadir in a logistic regression model as covariates (Heaton et al., 2010). In contrast to the consistent effects on CSF HP levels seen across both ancestry groups, we observed opposite directions of effect of the HP2 allele on NCI in AFR and EUR individuals. As shown in Table 2, in 395 AFR participants, the HP2 allele was associated with decreased odds of NCI (Odds Ratio, OR = 0.584, p = 0.022), while in 371 EUR participants, the HP2 allele was associated with increased odds of NCI (OR = 2.081, p = 0.040).

To account for uncertainty in genotype imputation estimations, we also tested the association between NCI and HP2 using directly imputed dosages rather than hardcalls. In EUR, we were able to detect a significant additive effect of HP2 dosage with OR = 1.605, (p = 0.002), which is similar to what we obtained using hardcalls. No significant associations were found in AFR using directly imputed allele dosages.

We also found that the effect of HP2 on NCI is independent from its association with CSF HP protein levels in both EUR and AFR PLWH. Although the sample size was limited, we performed sensitivity tests with the CSF HP protein levels to see if the effect of the HP CNV is driven by changes in the HP protein levels. Adjusting for the CSF HP levels did not change the significance of the HP2 alleles in either EUR or AFR PLWH. From Supplementary Table S1, we could see that the CSF HP levels is not significantly associated with NCI risk in either AFR or EUR individuals. In fact, this adjustment strengthened the HP association in EUR (Supplementary Table S1).

Given the longitudinal study design of CHARTER, we were also able to examine the effect of the HP CNV on NCI as the study participants age. For the AFR group, the median follow-up time is 6 months and the mean is 1.406 years. While in EUR, the median follow-up time is 6 months and the mean is 1.679 years.

As shown in Table 3, EUR PLWH having HP2 were at 10.276 higher OR (p = 0.037) of developing NCI than PLWH that do not have HP2 (Figure 2A). This OR decreased by 0.961 per year (p = 0.077, Figure 2B). This decrease reflects a cumulative change in risk; for example, a decrease of 0.961 means the individuals probability of NCI equals the probability of the previous year × 0.961. We also noticed that though the risk of NCI for EUR HP1/HP1 individuals was lower than HP2 individuals at the beginning, it increased faster and became greater than HP2 individuals after approximately 58 years of age. In contrast, though only marginally significant, AFR PLWH with HP2 had lower risk (OR = 0.111, p = 0.054) of developing NCI compared to people without HP2 (Figure 2A). The directions of the HP2 effect on NCI were consistent with the baseline models within each population. These effects are visualized in Figure 2. We see a dramatic change in the trend of predicted NCI probability over age in both AFR and EUR between HP2 dominant and HP1/HP1 individuals. These effects were not observed from the analyses that combines the AFR and EUR participants as the opposing effects of the HP alleles in each population cancel each other out (Figure 2). We also detected an additive effect of HP2 on NCI in EUR as shown in Table 4. In EUR, each copy of HP2 allele was associated with a 4.719 increase in the OR (p = 0.009) of developing NCI (Figure 2C). This OR decreased by 0.972 per HP2 allele per year (p = 0.024, Figure 2D). These effects are not statistically significant in the AFR or combined analyses (Table 4).

HIV RNA levels were significantly different by HP genotype (Supplementary Figures S3–S6), suggesting that the risk for NCI may be influenced by altering viral load. To examine the potential mediation effect of the HIV RNA, we re-ran the analyses in the virus suppressed subgroup [⩽ lower limit of quantitation (LLQ)]. Despite changes in significance levels in some tests due to dramatic decreases in sample sizes, the directions of the effects remained consistent in the longitudinal analyses (Supplementary Tables S2, S3). Further analyses adjusting for the ⩽ LLQ status showed that the ⩽ LLQ status did not significantly contribute to NCI risk (Supplementary Tables S4, S5).

We further decomposed the global model of NCI risk by investigating which cognitive domains are significantly affected by the HP CNV by conducting longitudinal analyses using domain impairments as outcomes. The higher risk of NCI in EUR with HP2 is likely driven by an increased risk in the learning domain impairment. This effect is slightly (but significantly) offset by a significant decrease in risk for verbal domain impairment based on an additive model of the HP alleles. These two effects, combined with more modest (non-significant) effects in other domains together create an overall increase in risk for NCI (Figure 3, Supplementary Tables S6, S7). In EUR, PLWH with HP2 had a higher OR of impairment in the learning domain (OR = 14.526, p = 0.019, Supplementary Table S6 and Figure 3A) with the risk reducing over time (by OR = 0.951, p = 0.030 per year of age, Supplementary Table S6). These over-time effects are reflected in Figure 3B where we observed a dramatic difference in slope between HP2 dominant and HP1 individuals. No statistically significant effects were noted in the AFR and combined (AFR + EUR) analyses (Figure 3B). We also found that with each HP2 allele, EUR PLWH have a significant decrease in risk of verbal domain impairment (OR = 0.259, p = 0.019, Supplementary Table S7 and Figure 3E). This effect is smaller than the contrasting effect of HP2 on the learning domain, and over time this risk modestly increases by OR = 1.020 per year (p = 0.023, Supplementary Table S7 and Figure 3F). In AFR, the HP2 allele was associated with a lower risk of impairment in the speed of information processing domain (Figure 3 and Supplementary Table S8) with this risk increasing by OR = 1.072, per year (p = 0.038, Supplementary Table S8 and Figure 3D).

From Figure 3B, similar to the global NCI risk, we also observed an intersection point of the trend lines between 50 and 60 years of age, after which the risk of learning impairment for EUR HP1/HP1 individuals became higher than HP2 individuals. This alteration was also found in the verbal and speed of information processing domains. In both EUR and AFR, HP is associated with some additional changes to cognitive domain T-scores that did not impact impairment (see Supplementary Tables S9, S10).