Host genetic variation modifying HIV-1 acquisition risk can inform development of HIV-1 prevention strategies. However, associations between rare or intermediate-frequency variants and HIV-1 acquisition are not well studied. We tested for the association between variation in genic regions and extreme HIV-1 acquisition phenotypes in 100 sub-Saharan Africans with whole genome sequencing data. Missense variants in immunoglobulin-like regions of CD101 and, among women, one missense/5’ UTR variant in UBE2V1, were associated with increased HIV-1 acquisition risk (p = 1. 9x10-4 and p = 3. 7x10-3, respectively, for replication). Both of these genes are known to impact host inflammatory pathways. Effect sizes increased with exposure to HIV-1 after adjusting for the independent effect of increasing exposure on acquisition risk. Trial registration: ClinicalTrials. gov NCT00194519; NCT00557245
The discovery of the protective deletion variant, CCR5-delta32, in the chemokine receptor 5 gene, encoding a HIV-1 co-receptor [1–3], generated great enthusiasm to search for additional host genetic variants and pathways associated with HIV-1 acquisition as a means of identifying targets for new HIV-1 prevention and treatment strategies. This enthusiasm was further enhanced by documentation of HIV-1 exposed seronegative (HESN) individuals who had very high exposure and lacked CCR5-delta32, [4–7] suggesting existence of additional genetic factors that alter the risk of sexually transmitted HIV-1 infection (Online Mendelian Inheritance in Man [OMIM] phenotype #609423) [8]. At least one in vitro experiment supports the hypothesis of a strong genetic component to infection risk, reporting 50% heritability in cellular susceptibility to HIV-1 infection [9]. Nevertheless, genome-wide association studies (GWAS) to date searching for such genetic risk factors for HIV-1 infection risk have met with limited success [10–17]. Most GWAS have had moderate (~80%) average power to detect common variants, very low power to detect variants with minor allele frequency (MAF) = 5% (approximately 1% power for an OR = 2) and even less power to detect associated rare variants (RVs) (MAF≤1%). Additionally, susceptibility to HIV-1 can only be assessed among individuals who are exposed to the virus. While some assessment of HIV-1 exposure was used for most HIV-1 GWAS [10–17], exposure measurement error and/or exposure misclassification, including that related to lack of information about the HIV-1 infected partners’ plasma HIV-1 RNA level (s), can result in lower statistical power than anticipated. Furthermore, out of necessity, some studies have been forced to attempt replication across different ancestral/racial groups [10]. However, risk variants might differ between such groups, thereby lessening the power for replication. Hence, major gaps in HIV-1 genetic association studies still exist, and focus on power to detect rare associated variants employing high accuracy in HIV-1 exposure measurements is warranted. The contribution of rare variants to risk of HIV-1 acquisition is of particular interest because effect size is generally inversely correlated with MAF when an association does exist [18] and large effects can be expected to translate to strong interventional impact (e. g. , HMG-CoA reductase inhibitors and familial hypercholesterolemia [19]). With these issues in mind, we undertook an association study of HIV-1 acquisition using whole genome sequencing (WGS) of extreme phenotypes sampled from two large clinical trials and one observational study of African HIV-1 serodiscordant couples (stable heterosexual couples with one partner HIV-1-infected and the other partner HIV-1-seronegative at enrollment) (n = 8,593 couples). These studies (S1 Table) had similar clinical follow-up, including quarterly risk assessments, PCR-verification of HIV-1 infection, reports of protected and unprotected sexual activity from both partners, measurement of the infected partner’s plasma HIV-1 RNA level and molecular confirmation of transmission linkage through viral sequencing [20–22].
The 100 Discovery stage genomes were sequenced by Complete Genomics Inc. (CGI) with high quality results (S2 Table). The RVT1 test (“rare variant test 1”) [24], a statistical test designed specifically for rare variant association studies, was used to test the difference between extremes in functional variant burden (defined as the total number of minor alleles (aggregated variant scores) by gene comparing cases and controls; see Materials and methods) for each of 18,354 genic regions, including 284,632 functional variants in the tests. The regions with the two lowest RVT1 p-values (S2A Fig) had an estimated 83% probability that at least one of these was a true positive, based on a False Discovery Rate (FDR) analysis [25]. QQ-plots of the RVT1 results showed good adherence to expected behavior and no evidence of confounding by major ancestry nor by spatially-isolated pockets of ancestry [26] (S2B and S2C Fig). These two regions are transcribed regions of CD101 (NCBI Gene ID: 9398) and UBE2V1 (NCBI Gene ID: 7335). For both of these genes, individuals with greater numbers of polymorphic functional sites had increased risk of HIV-1 acquisition: CD101 odds ratio (OR) = 2. 7 (per functional site with at least one minor allele), 95% CI = [1. 6–4. 8], p = 3. 6x10-5, and UBE2V1 OR = 3. 7,95% CI = [1. 8–7. 5], p = 4. 7x10-5 (Table 2). Within these two genes in the Discovery stage, eight variants were novel at the time of identification, and all eight validated by Sanger sequencing (S3 Table). Based on the FDR results showing >80% chance of a true positive, the CD101 and UBE2V1 regions were moved forward to the Replication stage for testing in a longitudinal analysis. Test regions from the Discovery stage almost certainly include variants significantly associated, as well as, variants not significantly associated (i. e. , noise) with outcome. The benefits, and perhaps even necessity, of using biological knowledge to increase signal-to-noise ratio for rare variant replication is well recognized [27,28]. To increase the signal-to-noise ratio and increase power, we prioritized and grouped variants (see Materials and methods) in CD101 and UBE2V1 for Replication stage testing. Fourteen functional variants in CD101 were designated as “primary replication variants” (PRVs) based on a direction of effect consistent with that for the Discovery result (S4 Table) and their degree of significance in by-variant tests (S3A Fig). These variants were subdivided into four sub-groups for replication testing (see Materials and methods): (1) five missense variants in regions encoding extracellular CD101 immunoglobulin-like (Ig-like) protein domains [29], (2) five missense variants in the CD101 cytoplasmic domain, (3) two 3’-UTR variants and (4) two splice site variants (Fig 1A, S3A Fig, S4 Table). Creating four separate a priori replication test groups increases the multiple-testing penalty, but is expected to further increase the signal-to-noise ratio within some of these variant test groups. Similarly, 11 of 15 predicted functional variants in UBE2V1 were designated as PRVs, and these were divided into two groups: (1) six 5’-UTR and (2) five 3’-UTR variants (Fig 1B, S3B Fig, S5 Table). One of the 5’-UTR variants, rs6095771, is identified as a missense variant in the canonical transcript, and as a 5’-UTR for other transcripts. It was grouped with the 5’-UTR variants due to low predicted power to replicate a single variant with MAF = 0. 02. Three of these variants were novel and another six are not found among Kenyans in the 1kG database [30] (S5 Table). Furthermore, two of the functional UBE2V1 Discovery variants that are indexed in the ESP [31] or ExAC [32] databases are extremely rare and found in ESP only among individuals of African descent (rs187204768, MAF = 4/4046 and rs6095771, MAF = 20/4400) (S5 Table). These results indicate that the UBE2V1 association could not be detected at the variant level using GWA-type methods and can only be detected in a population with substantial recent African ancestry, such as the present study. Across these three HIV-1 serodiscordant couples cohorts, the pool of participants who were available for the Replication stage had varying amounts of reported sexual exposure. Because inclusion of individuals with no or little HIV-1 exposure could reduce power, we identified and excluded these individuals using a Protected-sex Index (PI). PI is defined as the proportion of study visits for which only abstinence or 100% condom use was reported; we considered this our measure of baseline behavior tendency. PI predicted overall HIV-1 seroconversion rates among the three HIV-1 cohorts (Fig 2), indicating reasonably low measurement error and a high signal-to-noise ratio for PI. Simulation studies predicted that power would be maximized when only individuals with PI ≤0. 6 were included in the Replication analysis, despite a larger sample size when individuals with lower exposure are included (see Materials and methods). Accordingly, the Replication cohort was restricted to the 261 individuals HIV-1 uninfected at baseline, who each had PI ≤ 0. 6 and were not selected for the Discovery stage analysis. Genotyping of variants in CD101 and UBE2V1 for the 261 Replication stage individuals was completed using molecular inversion probe sequencing (MIPs) [33] technology (see Materials and methods). Given the cost efficiency of using MIPs, an auxiliary sample of 968 participants selected by more common methods rather than sexual exposure also was sequenced and used for verification of simulated power studies (see Materials and methods) for a total of 1229 individuals with MIPs data for CD101 and UBE2V1 (S6 Table). Neither CD101 splice variant from the Discovery stage was found in the Replication sample (S7A Table), reducing the number of CD101 primary test groups from four to three. In total, 83 CD101 SNVs were detected among the 1229 individuals (S8 Table). For UBE2V1, three of five 3’-UTR PRVs were found in the Replication sample but only one of the six 5’-UTR PRVs was present (S7B & S9 Tables). This is consistent both with the rarity of these UBE2V1 variants and enrichment for these variants in the Discovery sample. To test for further evidence of a true positive association between HIV-1 acquisition risk in women and Ig-like CD101 variants or UBE2V1 rs6095771, we augmented the Replication data with variant data from the auxiliary sample and tested for an interaction between the reported PI score and each PRV score to ascertain whether a dose-response relationship was present. “Dose” here is HIV-1 virus exposure quantified in terms of frequency of unprotected sex with a partner with the average plasma HIV-1 RNA in the Replication sample partners. For UBE2V1 rs6095771, this increased the number of women with 5’-UTR PRVs in the analysis from 4 to 28 (S7B Table). We found strong dose-response relationships for both genes, indicating that the association between variant scores and HIV-1 risk is positively related to the frequency of unprotected sex: p = 5. 0x10-8 for the CD101 Ig-like dose-response and p = 8. 2x10-6 for the UBE2V1 rs6095771 dose-response, with exposure assessed through the PI score and modeled on a log-scale (Fig 4; S7 Fig). Because the model is adjusted for PI, these significant associations with increasing dose are in addition to the effect of decreasing PI. Overall model significance levels including the variant and dose-response variables were marked at p = 5. 1x10-12 and p = 1. 7x10-12, respectively, for CD101 and UBE2V1 in explaining variation in HIV-1 seroconversion risk (Tables 4 and 5). The interaction models also provide estimated HRs for the risk groups under the assumptions of no protection and the average frequency of heterosexual intercourse among those in the Replication Stage (CD101 HR = 5. 4, p = 1. 4x10-3,95% CI = [1. 9,15. 2]; and UBE2V1 HR = 16. 2, p = 2. 8x10-4,95% CI = [3. 6,72. 6]) for the Ig-like PRV score and UBE2V1 rs6095771, respectively. Again, these HRs are adjusted for PI, which means that the increase in effect size is above that explained by the PI variable. When we included all 1229 samples genotyped after the Discovery stage (Replication sample and auxiliary sample) in a model that does not account for sexual exposure/PI we found highly attenuated HRs (Tables 4 and 5) and non-significant p-values. This phenomenon occurs because the estimated HRs are averages over a group that contains individuals who have no risk due to no sexual exposure, “diluting” the effect we detected in those with higher exposure. To illustrate these effects as a function of the exposure to HIV-1 viral quantity, HRs were estimated with step-wise decreases in sample size, excluding exposed individuals at each step (Fig 4, S6 & S7 Figs). Levels of cytokines in blood are a useful measure for immunologic function and may indicate presence of a generalized host pro-inflammatory state [36]. We evaluated whether the three CD101 Ig-like PRVs with FDR < 0. 05 in the replication stage were also associated with altered plasma cytokine levels compared to individuals without any of these variants (see Materials and methods). Among the 163 individuals in this subset for whom CD101 genotypes were determined and for whom measurements of 25 plasma cytokines were available [37], carriers of these CD101 risk alleles had significantly lower serum levels of IL1R1 (right-shifted distribution) compared to those without any CD101 risk alleles (OR = 0. 19 for achieving the 75th percentile IL1R1 value, 95% CI = [0. 07,0. 54], p = 1. 7x10-3; adjusted p = 0. 04) (S8 and S9 Figs). There was also a tendency toward lower levels of sCD40L (p = 0. 0049; S10 Table). Frequencies of the UBE2V1 5’-UTR PRVs were too low in the individuals with cytokine measurements to allow for assessment of any association with cytokine levels. None of these variants were associated with HIV-1 plasma RNA set point among HIV-1 seroconverters (S11 Table), indicating that the associations we have found are unlikely to act through altered viral replication with these gene products functioning as intracellular viral restriction factors [38].
Our findings demonstrate that aggregates of host genetic variants, including variants with MAF<10%, can have a strong and replicated association with HIV-1 acquisition risk. Our replication of CD101 Ig-like variants showed an HR for HIV-1 infection of 4. 3 (95% CI = [2. 1–8. 9], p = 6. 4x10-5); and replication of association of UBE2V1 rs6095771 with HIV-1 acquisition risk showed an HR of 6. 4 in women (95% CI = [2. 1,19. 1], p = 9. 5x10-4). Two CD101 missense/regulatory variants reached individual statistical significance after adjustment for multiple testing, four had individual FDRs < 0. 05, and a strong dose-response relationship with the frequency of unprotected sex collectively strengthened evidence that association with HIV-1 acquisition is real. While many of these variants were individually rare or infrequent, nearly 18% of Kenyans evaluated in the 1K Genomes Project [30] had one or more of the three most common CD101 Ig-like HIV-1 risk variants (rs17235773, rs3754112, rs12093834). Hence, these variants or others in CD101 and UBE2V1 could have a substantial impact on the population-based HIV-1 infection risk. CD101 and UBE2V1 are both biologically plausible candidates for influencing HIV-1 acquisition. CD101 is expressed on CD4+ and CD8+ T-cells, dendritic cells and monocytes, [39] and appears to alter CD4+/CD25+/FOXP3+ T regulatory cell (Treg) function based on both a murine graft-versus-host disease model [40], and through IL10 secretion from human dendritic cells [41]. Monoclonal antibody ligation of CD101 reduces T cell proliferation through a Ca2+ and tyrosine kinase-dependent pathway possibly by preventing translocation of nuclear factor of activated T cells (NFAT) and IL-2 production [42]. Recent experiments using a murine model of chronic colitis demonstrate that adoptive transfer of CD101-/- Tregs is associated with Th17 cell proliferation and more severe colitis [43]. CD101 expression is also strongly associated with the immune suppression function of Tregs in humans. [44] Reduced expression of CD101 on mucosal CD8+ T cells has been associated with increased tissue inflammation in studies of human intestinal [45], and pulmonary mucosa. [46] Given that local inflammation and CD4+ [47] and CD8+ [48,49] T cell immune activation have been associated with increased HIV-1 acquisition risk, and reduced immune activation [50,51] or immune quiescence [52,53] have been associated with natural resistance to HIV-1 in HESN, our results support the idea that CD101 gene variants may modify HIV-1 heterosexual acquisition risk through altered levels of genital mucosal inflammation. The finding that CD101 risk variants are associated with lower plasma IL1R1 levels, but not with HIV-1 RNA set point, suggests that these variants have systemic immunological effects in the seronegative partner while not directly acting on HIV-1 replication. Recent studies in mouse models indicate that IL1 may inhibit Treg and enhance Th17 differentiation [54]. However, it is unclear at this point how variants in CD101, specifically those identified in CD101 Ig-like domains, might modify either IL1 or IL1R1 levels. The IL1R1 rs2234650 genotype has been reported to be associated with HIV-1 acquisition in infants of HIV-1 infected mothers with modification of risk by IL1 gene family haplotypes [55]. In addition to CD101 being associated with increased Treg function, IL1R1 expression has been associated with increased Treg and anti-inflammatory IL10 secretion [44]. Combining these prior data with our findings, we hypothesize that CD101 Ig-like variants reduce Treg function with associated reductions in IL1R1 levels and an enhanced pro-inflammatory environment, which leads to increased risk of HIV-1 acquisition. While efforts to test this hypothesis and to develop a more detailed understanding of CD101 function are underway, our results suggest that targeting CD101 activity could be a novel approach to host-directed, HIV-1 prevention. UBE2V1 associates with TRIM5-α, a host restriction factor involved with HIV-1 capsid uncoating [56]; however, rare and uncommon UBE2V1 variants have not previously been studied in association with HIV-1 acquisition risk. Previous GWAS without rare variant burden/aggregation tests cannot adequately assess the UBE2V1 association observed here because of the rarity of the variants found to be associated with HIV-1 acquisition risk in this study. UBE2V1 forms an ubiquitin-conjugating complex generating unattached polyubiquitin chains that may stimulate NF-κB activation and consequent pro-inflammatory cytokine production [56,57]. This complements reports of reduced systemic immune activation associated with resistance to HIV-1 acquisition in Kenyan sex workers [50,53,58]. Of note, UBE2V1 is one of only 23 genes differentially down-regulated by HIV-1 trans-activator of transcription (TAT), an HIV-1 protein that is required for efficient replication of the HIV-1 virus and potential escape from the host immune system [59]. Distinguishing between a group carrying variants protecting against HIV-1 acquisition versus a group carrying risk-increasing variants requires that both groups are exposed to HIV-1, i. e. to assess susceptibility to HIV-1 in any given individual, that individual must be exposed to HIV-1. This translates directly to a mathematical proof that statistical power to detect a significant association with HIV-1 acquisition risk is increased by selectively identifying individuals with sustained high levels of HIV-1 exposure. Indeed, the success of this study was dependent on our ability to quantify HIV-1 exposure with relatively high accuracy to identify exposed individuals. Given that the overall per-contact probability of heterosexual HIV-1 transmission is intrinsically low (estimated at ~1/1000 for vaginal intercourse among the study population) [60], our findings suggest that inclusion of individuals with lesser HIV-1 exposure may explain, in part, why it has been difficult to identify and/or validate genetic risk factors for HIV-1 acquisition. The use of the extreme phenotypes design is another strength of this study. The discovery of CCR5-delta32 was based on the observation that extreme resistance phenotypes existed [61,62] leading to case-control studies to identify the CCR5-delta32 variant among candidate genes [1–3], including extreme hemophiliac controls [1] and the others employing just two to four extreme resistant controls [2,3]. Identifying individuals with phenotypes that represent the extremes of risk of HIV-1 acquisition has been challenging because few studies are able to assess both the behavioral and biologic dimensions (e. g. frequency of unprotected sex, and plasma HIV-1 RNA level in the HIV-1 transmitting partner) that contribute to exposure. Our use of data from both partners added accuracy to the phenotypes we used in this study—e. g. , plasma HIV-1 RNA level contributed by the HIV-1 infected partner, reported frequency of unprotected sex from both partners, and epidemiologic data from the HIV-1 uninfected partner (e. g. , male circumcision status) [23]. Further improvements in accurate exposure quantification could increase power even more and perhaps very extreme phenotypes can be identified that provide superb power with small samples but lead to generalizable treatments. Recent examples that underscore the value of the extreme phenotype analysis approach include discovery using a sample size of a dozen individuals of a human antibody to a malaria protein that prevents death and provides a promising new malaria vaccine target [63], and discovery of the protective effect of PKC9 loss of function variants against cardiovascular disease in a small group of extreme individuals [64], leading to the development of PKC9 inhibitors for lowering LDL cholesterol. Our analysis identified variants in CD101 and UBE2V1 associated with increased risk of HIV-1 acquisition, but we did not identify variants associated with reduced risk of HIV-1 acquisition. Statistical power is lower to identify protective variants than it is to identify risk variants (of the same magnitude but inverse) when the outcome has low incidence. Discovery of protective variants for HIV-1 infection in a population with average sexual contact and some use of protection against transmission would require extended observation time to detect differences in the survivor rates (“rates of non-seroconversion, ”) while differences in rates of seroconversion of the same magnitude can be detected statistically within a shorter time period. A limitation of these results is that confounding cannot be ruled out with certainty as the source of the associations. This is true for any observational study, though observational studies remain a fundamental part in building a step-wise scientific case for causal association for many exposures that cannot be tested experimentally in humans (including exposure to genetic variants, the epidemiological exposure being tested here), with smoking as a cause of lung cancer being a prime example. [65] We have guarded against confounding to the extent possible by evaluating for known potential confounders, including ancestry group (which likely includes HLA), BV, age, cohort, ethnic group affiliation and spatially isolated ancestry. The next general step in establishing causal association is replication by different groups and functional studies. We are currently engaged in the latter and encourage validation of our results in genetic association studies with HIV-1 exposure measurement as well as functional studies by others of the variants/genes discovered and replicated here. In summary, we used quantitative measures of HIV-exposure to select individuals with extreme HIV-1 acquisition phenotypes and thereby optimize our power to detect genes associated with risk of HIV-1 acquisition. We identified variants, including rare variants, in CD101 and at least one in UBE2V1 that are significantly associated with increased HIV-1 acquisition risk. More detailed dissection of the molecular basis for this association may identify unique interventions that use these pathways to improve public health prevention of HIV-1.
We identified individuals for this study from HIV-1 serodiscordant couples recruited into three cohorts of African heterosexual HIV-1 serodiscordant couples: the Partners in Prevention HSV/HIV Transmission Study [20] (ClinicalTrials. gov number, NCT00194519), the Couples Observational Study [14], and the Partners PrEP study [22] (ClinicalTrials. gov number, NCT00557245) (S1 Table). Detailed procedures have been reported elsewhere for each of these studies [14,20,22]. Briefly, routine follow-up visits with both partners were scheduled at least every 3 months, with clinical, behavioral and demographic data collected. HIV-1 seroconversion (SC) was assessed by HIV-1 rapid test at the study clinic; positive rapid tests were confirmed by HIV-1 ELISA at the site laboratory, and by Western Blot in batch at the University of Washington (UW). Plasma virus sequencing performed on both partners for each couple associated with SC was used to confirm transmission linkage [66]. All participants provided written informed consent for participation in the clinical study, and samples for this genotyping study were selected from among those participants recruited at 14 sites across all three cohorts who additionally consented to host genetic studies. Relevant study documents went through ethical review and approval by the following committees: An extreme phenotypes case-control design was employed for the Discovery stage. The extreme phenotypes design provides the greatest statistical power for a fixed Discovery sample size of 100 individuals (with power increasing as the percentiles of phenotype become more extreme in the two arms). Hence, individual phenotype was a primary consideration in Discovery stage participant selection. Extreme cases here comprise individuals with relatively low estimated exposure who seroconverted during the study, especially those who converted early in the study. Extreme controls comprised individuals with high estimated risk who remained seronegative over the full observation period with follow-up for at least nine months. Highly-exposed HIV-1 exposed seronegative control individuals with longer follow-up time were considered more extreme based on cumulative exposure scores across all study visits. Cumulative exposure score to rank extremes was calculated as previously described using plasma HIV-1 RNA level of the infected partner, frequency of unprotected sex and male circumcision status [23], with the modification that all participants must have reported unprotected sex by at least one partner in the couple to be eligible for the Discovery stage sample. It is possible for an individual to have a fairly high exposure score in this model even if no unprotected sex is reported, because the original risk score is based on empirical estimates of risk of seroconversion given the variable values, and some participants who reported no unprotected sex did seroconvert with plasma HIV-1 genomes matching those of their infected partners. Hence, the risk among the group that reported “never unprotected sex” is not zero. Nevertheless, the latter risk is smaller than that for those who report unprotected sex, all else equal, and the risk rises as the proportion of reports with unprotected sex rises (S6 Fig). Seronegative individuals from couples that did not report unprotected sex were excluded for statistical power reasons: including unexposed individuals lowers the statistical power to the point where a p-value of 6. 3x10-5 in the Replication stage (N = 261) (Table 3) becomes p = 0. 03 (N = 1229) (Table 4; replication stage methods below). Given these exposure scores/conditions defining the extremeness of phenotype within the potential cases and controls, we then incorporated both gender balance within-group and pairwise ethnicity/sex matching between-groups to this design to avoid happenstance confounding by differing proportions of sex or ancestry in the relatively small samples. The potential impact of ancestry-confounding in this subpopulation was largely unknown when the Discovery stage was designed, and we opted to take this precaution against confounding. Specifically, seroconverter (SC) “cases” were selected from the Partners in Prevention HSV/HIV Transmission Study and COS cohorts among couples with laboratory confirmed linked HIV-1 transmission [67], who were HIV-1 polymerase chain reaction (PCR) negative at enrollment and had the lowest exposure scores, conditional on relatively equal numbers of males and females. Individuals from the Partners PrEP cohort (S1 Table) were not available for Discovery stage sampling because the trial had not come to complete closure with data available for ancillary studies at that time. For selection of Discovery stage controls, we excluded from consideration: SC individuals (with either linked or unlinked transmissions), couples in which HIV-1 infected partners reported use of any antiretroviral therapy (ART), couples with no reported unprotected sex during the study and couples with less than nine months follow-up time. After these exclusions, for each selected case, all HESN of matching sex and self-reported ethnicity were identified, and from these matched individuals, the HESN individual with the highest cumulative exposure score was selected as the matching control. A total of 65 case-control pairs were identified in this manner, with identification of low exposure (extreme) cases being the highest selection priority, followed by ethnicity/sex matching for controls, followed by criteria for high exposure among the matching controls. We quality controlled these case-control pairs for gender check, cryptic relatedness and genetic heterogeneity across multiple longitudinal whole blood DNA samples, using a custom Illumina Goldengate chip with 384 single nucleotide polymorphisms (SNPs). These test SNPs were selected as being the most predictive of ancestry clusters and individual identity from a principal component analysis (PCA) on data from a previous genome-wide association study [14] that included samples from these same cohorts. Specifically, we first performed a PCA on a pruned set of 133,991 SNPs from that GWAS that had low linkage disequilibrium. The first five PCs for this analysis were effective at distinguishing participants from East and southern Africa, by country (Kenya, Uganda, Tanzania, South Africa and Botswana) and by self-reported ethnicities reported in >2% of participants (S10 Fig). Using the first five PCs, we then assigned all participants to one of nine ancestry clusters based on model-based clustering, which has previously been shown to reduce population stratification bias. Subsequently, we used the Random Forests algorithm to identify 357 SNPs that were most predictive of geographic region (East Africa versus southern Africa) and the nine ancestry clusters. These 357 SNPs were collectively able to differentiate the ancestry clusters but were much less important for predicting ancestry than self-reported ethnicity and geographic region (S11 Fig). The final Goldengate SNP chip included the 357 ancestry SNPs along with 27 SNPs that maximized the probability that all participants had a different genotype at one or more loci in order to ensure that DNA samples came from unique individuals. The genotyping chip was used on DNA from 65 potential case-control pairs. After eliminating samples that failed QC (1 case failed cryptic relatedness, and longitudinal samples from 2 controls suggested potential sample heterogeneity), verifying matching on ancestry cluster and identifying controls with highest cumulativeHIV-1 exposure scores over all visits, 50 case-control pairs were selected for complete genome sequencing for the Discovery analysis (24 male cases with matched controls, and 26 female cases with matched controls). When characterized by the two strongest components of the exposure score, namely mean plasma HIV-1 RNA and proportion of follow-up visits where no-condom use was reported, the Discovery stage controls were verified as sampled from the highest exposure strata (Table 1). Cases also had a median time to HIV-1 seroconversion of 11. 8 months, while the median duration of follow-up without seroconversion for controls was 22. 8 months. WGS was performed by Complete Genomics, Inc (CGI), using published methodology [68], and with samples blinded as to case and control assignment. SNV calling was performed by CGI using proprietary software cgatools version 1. 5. 0 build 31 (dev) and human reference genome NCBI Build 37. Overall, CGI sequence quality and aggregate descriptive characteristics were similar for case and control genomes (S2 Table). Genome annotations were from CGI based on the National Center for Biotechnology Information (NCBI) refSeq database. Annotations used for this rare SNV analysis included non-synonymous protein coding sequence, 3’- and 5’-untranslated regions, and splice donor/acceptor sites with exon/intron boundaries identified through the University of California-Santa Cruz (UCSC) refFlat database. [69] To validate WGS sequence data we identified 8 SNVs in the WGS data that lacked an rsID in dbSNP and were identified in either of the two genes targeted for replication analysis (CD101 and UBE2V1). SNVs not previously reported to dbSNP have a higher probability of being false positives than the bulk of SNVs in dbSNP, so these provide “tough” test cases for the CGI calls for validation of true positive rare variants. These SNVs were Sanger sequenced from the individual genomes originally used to generate the WGS. All 8 were successfully validated (S3 Table). Prior to analysis, each SNV was classified as “functional” or not, based on CGI annotation indicating that the SNV either (1) altered protein coding, (2) was at a splice-site or (3) in a UTR sequence region. Only these genic categories of functional polymorphisms were included in this analysis. Case-control comparison of SNVs aggregated by genic region (“variant burden”) was accomplished through logistic regression analyses with burden scoring derived specifically for rare variant analysis by Morris and Zeggini with the p-value based on the likelihood ratio test (their RVT1) [24] The Wald test p-value was found to be unreliable—much too conservative—for genic regions with high imbalance between cases and controls, which are exactly the situations we seek to find. Although cases and controls were matched by ethnic origin, the first 3 principal components (PC) in the logistic regression of each genic region were used to control for potential residual population stratification. To focus on less common or rare variants and higher effect sizes, the RVT1 ignores variants with MAF above a specified cutoff (>0. 125 in our analysis). The MAF cut-off of 0. 125, which is somewhat higher than typically considered rare, was used to allow for random variation in the observed MAF for rare SNVs in the population, as well as to allow for enrichment of rare SNVs in the sample due to the selection of extremes. Because the study participants are sub-Saharan Africans, external estimates of MAFs for the observed SNVs were not available. This testing prioritization strategy resulted in 18,354 genic regions tested (each region with an RVT1 p-value produced for the null hypothesis of no difference in variant burden within the region between cases and controls) and included a total of 284,632 functional variants in these tests. Principal components created by LD-pruning variants with MAF > 0. 03 were used to adjust the RVT1 tests. Genic regions with the lowest p-values were considered for replication. These also have the lowest False Discovery Rates [70] (FDRs). Based on FDR analysis, the two most significant regions had a probability of 0. 83 that at least one was a true positive (not a false discovery). These two regions were moved forward to the MIPs Replication stage based on this high probability of having a true positive result. We planned and implemented a formal frequentist p-value-based criterion for declaring significance of the replication of CD101 and UBE2V1 genic region variants contributing to the aggregate scores in the Discovery stage. We broke the variants into smaller aggregates for the replication analysis (four aggregates for CD101, two for UBE2V1) in an a priori attempt to further enrich at least some of these aggregates with a higher proportion of true positive variants. This reflected the fact that Discovery stage test regions include both variants that are significantly associated with outcome (i. e. , signal) and those not significantly associated with outcome (i. e. , noise). If sub-groups of the originally aggregated variants can be enriched for truly associated variants having a common direction of effect, the statistical power for replication of such an enriched sub-group can be increased considerably relative to the larger, noisier aggregate [19,28] despite the need for a multiple-test correction due to increasing the number of tested aggregated variant groups. Based on statistical reasoning, variants with large effect sizes and the same direction of effect are most likely to be driving the Discovery finding. To formulate sub-groups prior to the Replication analysis, “by-variant” association results were tabulated and examined for size and direction of effect. Confidence intervals are wide for these by-variant results due to the low MAF/low-power issue, but they are not completely uninformative. To complement these results, we also compared by-variant results to genome and protein maps of the regions (Fig 1; S4 & S5 Tables). Sixteen of 24 predicted functional variants in CD101 had by-variant estimates with increased risk (HRs>1), which is the direction of effect for the aggregate RVT1 CD101 result that we seek to replicate. Fourteen of these were designated as PRVs based on significance of by-variant tests (S3A Fig). These 14 CD101 PRVs were divided into Ig-like (N = 5 PRVs), cytoplasmic (N = 5), 3’-UTR (N = 2) and splice site domain (N = 2) subgroups. Similarly, 11 functional variants in UBE2V1 identified as PRVs (S3B Fig) divided into 5’-UTR (N = 6) and 3’-UTR (N = 5) subgroups. Although a fairly common strategy in replication of rare variant results is to include all rare variants within a region in the replication test whether or not these variants were among those identified in the Discovery stage. For example, such a strategy would be used to study low-density lipoprotein receptor (LDLR) variants [27,28] (OMIM #606945), given that the structure and function of LDLR is well-known, providing reasonably high confidence that novel missense or truncation RVs in particular sub-regions of the LDLR gene will have a deleterious association. In this study, we elected not to include in replication tests any variants seen only in the Replication stage and not in the Discovery stage for four reasons: (1) identification of individual variant contributions is important at this point to identify variation in any specific structural protein components affecting the outcome (and point toward mechanism); (2) identification of individual variant contributions for variants that are only modestly rare makes it feasible to use genotype tests for presence of individual risk variants in future studies or risk assessments; (3) addition of new variants in the tests at the Replication stage obscures the contribution of variants from the Discovery stage alone; and (4) we did not expect a large contribution from new/novel/rare variants in the Replication stage because of the sample size limitations at this point, and therefore were not concerned by potential loss of statistical power via this choice. In addition to the primary replication tests that include subsets of variants identified in the Discovery stage (four tests planned for CD101 and two for UBE2V1), we performed exploratory testing on all additional variants found in the Replication stage to take the most advantage of these findings while limiting the formal replication tests as above. All exploratory statistical test results are reported by nominal p-value and FDR to provide some quantification of evidence for association, but no formal multiple test correction is given for these findings given their often highly correlated nature (making a Bonferroni correction misleading) and due to the post hoc character in some cases. We measured sexual exposure behavior defining Protected-sex Index (PI) as the proportion of study visits for which only abstinence or 100% condom use was reported. Simulation studies showed that statistical power was highest when the Replication cohort included individuals with PI ≤0. 6 while individuals with higher levels of protected sex were excluded. We set the Replication cohort to be the 262 baseline uninfected individuals each having PI ≤ 0. 6 who remained after the Discovery stage. Due to the economy of sequencing multiple individuals using multiple inversion probe sequencing (MIPs) [33] and the intention to demonstrate proof of principle in our power calculations, we successfully sequenced an additional 986 individuals (thirteen 96-well plates). This set of 1248 individuals was balanced on sex and ethnicity relative to seroconverters and non-seroconverters at a 1: 6 ratio, a typical sampling strategy. These additional 986 individuals had higher exposure scores than the average participant in the serodiscordant couple cohorts based on high plasma HIV-1 RNA levels of their infected study partner, but many reported low levels of or no unprotected sex. Therefore, this group of additional participants was not considered by itself. There were 19 individuals failing the MIPs procedure due to low DNA concentration, resulting in 1229 individuals with 261 in the Replication sample and 968 in the auxiliary sample. MIPs can be used for targeted sequence capture, followed by massively parallel sequencing of captured products. This strategy is efficient and cost-effective for sequencing multiple candidate genes in modest to large sample sets [33]. MIPs were designed to target all coding exons plus 10 additional base pairs of intron/exon flanking sequence for CD101 (RefSeq NM_004258), and UBE2V1 (RefSeq NM_021988). In total, 74 MIPs (Integrated DNA Technologies, Coralville, IA) were designed and pooled in equimolar ratios and a test library of 24 African control samples from Centre d’Etude du Polymorphisme Humain (CEPH) and 23 Caucasian in-house control samples was produced and evaluated. For MIPs that failed to produce the minimum required average sequence depth of 60X in the test library, their concentration in the final pool was increased according to their level of under-performance in order to increase coverage above the average threshold. The optimized MIP pool (S12 Table) was then used to generate libraries and targeted sequence from the ethnically-matched Replication Sample. The MIP pool was phosphorylated with T4 Polynucleotide kinase and T4 Ligase Buffer (New England Biolabs, Ipswich, MA). The reaction was held at 37°C for 45 minutes and then de-activated at 65°C for 20 minutes. For the majority of samples (N = 1096), 100ng of genomic DNA per sample was used in a capture reaction with a 200: 1 MIP to gDNA ratio. For samples with less than 100ng of DNA available (N = 228), a range of 40ng to 100ng was used in the capture. During a 24-hour reaction at 60°C, HemoKlen Taq (New England Biolabs, Ipswich, MA) was used to capture the target regions and Ampligase (Epicentre, Madison, WI) was used to circularize constructs. E. coli Exonuclease I and Exonuclease III (New England Biolabs, Ipswich, MA) were used to enzymatically clean capture reactions. The reactions were held at 37°C for 30 minutes and then de-activated at 95°C for 2 minutes. Captured products were amplified using iProof HF Master Mix (Bio-Rad, Hercules, CA) during 20–22 cycles of PCR. This PCR was performed using a generic forward primer and reverse primers containing a generic portion along with a unique 8bp molecular tag used to barcode the captured products for each sample. Between 40 and 96 barcoded samples were pooled into each library and the libraries were purified using Agencourt AMPure XP (Beckman Coulter, Indianapolis, IN) magnetic beads in the ratio 0. 9: 1 beads to samples. Libraries were eluted into EB buffer (Qiagen, Valencia, CA) and quantified via the Broad Range Quant-iT dsDNA Assay Kit (Life Tech, Waltham, MA) using a SpectraMax Gemini XPS Fluorometer (Molecular Devices, Sunnyvale, CA). Libraries were then pooled and sequenced on an Illumina MiSeq using 300 cycle paired end (v2) reagents (Illumina, San Diego, CA). Custom oligonucleotides complementary to sequences in the MIP constructs were used. Each run contained between 136 and 192 samples. Libraries were diluted and denatured according to Illumina’s standard procedure with final loading concentrations ranging from 8 to 10pM. Individual fastq files were generated by the MiSeq Reporter Software (v2. 5). The resulting fastq files were aligned to Hg19 with the Burrows-Wheeler Aligner (BWA v0. 5. 9-r16), and a multi-sample variant call file was generated using the Genome Analysis Tool Kit (GATK v2. 4-9-g532efad). Annotated was performed using Variant Effect Predictor v82 (http: //www. ensembl. org/info/docs/tools/vep). The average sequencing depth per sample across all sites was 198X. A total of 1229 individuals’ sequences passed quality control. Given that a subset of these samples (N = 138) had less than 100ng of total DNA and therefore could have yielded high variant missingness in these samples but with good data quality in the remainder of samples, we set the threshold of 12% missingness for excluding a variant from the Replication stage tests below (compared to a conventional GWAS threshold of 5% missingness). A total of 314 variants passed quality control including a Hardy Weinberg Equilibrium test cut-off of p>3. 18x10-5. Testing for association between genomic variation and HIV-1 acquisition risk in the Replication stage was performed using straightforward, standard censored data methods, since follow-up time in the three cohorts was variable with most observations being censored. This provides greater statistical power than using a subset of individuals who have reached a minimum follow-up time or the outcome to create a case-control sample, as it allows all follow-up time to contribute to the estimated risk associations. Kaplan-Meier plots and Cox proportional hazards models were used. We performed the primary replication tests by scoring the aggregate risk variable as 1, if an individual carried a minor allele for any of the variants in the aggregate, and 0 otherwise. This variable was tested for association with time-to-seroconversion via the Cox model. Risk variables were tested jointly in the Cox model, as well, and results were checked for possible confounding by sex, cohort, age, country/ethnic group (five countries with Kenya broken into three major ethnic groups plus others, a surrogate for ancestry, S4 Fig). Principal components that were available on the subset of Replication individuals from the GWAS also were checked as potential confounders/adjustment variables on this subset. Because we and others [71,72] have found bacterial vaginosis (BV) to be associated with increased risk of sexually acquired HIV-1 infection, BV is a potential confounder if associated with either CD101 or UBE2V1. BV assessment (as we have previously defined it [73]) was available on a subset of Replication individuals and was tested as a potential confounder on this subset. Cytokine levels at baseline entry to the cohorts were measured in a nested case-control subset of individuals from a previous study [37], wherein laboratory methods are described in detail. We assessed for differences in these cytokine levels among CD101 variant carriers and non-carriers. This analysis includes all individuals who had both baseline cytokine measurements performed and were among the individuals included in the Replication stage MIPS sequencing or Discovery stage WGS in order to determine their CD101 variant genotypes. Three CD101 missense variants with increased-risk in the Discovery stage data also had high increased-risk in the Replication stage data and were all in the Ig-like aggregate that replicated overall (chr1: 117554421, chr1: 117560058, chr1: 117568500). We separated individuals with cytokine measurements into individuals with any minor allele for any of these three variants and individuals without any of the three alternate variants. Individuals in the second group might have other CD101 variants, but it is impossible to separate individuals with certainty into those with and without any CD101 causal variants; but this grouping will result in a conservative estimate of the difference between groups, as some individuals with causal CD101 variants might be included in the group “without causal” variants and diminish the difference between groups. (That is, any misclassification here will not produce false positive results but attenuate differences instead.) We then compared cytokine measurement distributions between these two groups. Since more than half of the cytokines had median values at the limit of detection in the overall group, we compared the two groups for a difference in the 75th percentile rather than the median. This test assesses for a difference in the cytokine distributions primarily using information from those values that are above the limit of detection. (Five cytokines measured in the original report [37] had 75th percentile values that were at the lower limit of detection and were not included in this analysis). Specifically, logistic regression was employed with the outcome being an indicator of whether a cytokine measurement was above the overall 75th percentile for both groups combined. The independent variable was the CD101 carrier status indicator and the regression was adjusted for the panel (batch) for each measurement. Several batches had significantly different means for some cytokines, and adjusting for batch generally increased the significance of the difference between gene variant groups. Mean differences in cytokine levels also were compared between groups to ensure that the logistic regression results were consistent with differences in means, though the latter will be affected by the large number of ties at the lower limit of detection.
Antiretroviral drugs for pre-exposure prophylaxis (PrEP) or treatment significantly reduce risk of HIV-1 acquisition and transmission, but face challenges of increasing access, maintaining high adherence, and selecting viral resistance. Improved understanding of the molecular determinants of HIV-1 sexual transmission could provide new public health HIV-1 prevention interventions. Factors proven to impact sexual HIV-1 transmission risk include epidemiologic exposure (e. g. , level of virus in the transmitting partner and frequency of unprotected sex), presence of genital inflammation, and host genetic variants common in the population. Rare or intermediate frequency genetic variants are an increasingly recognized reservoir of complex human disease-causing factors, but are not well studied in HIV-1 infection. However, the low frequency of these variants reduces statistical power to detect disease associations. Aggregating variants in a common biological domain (e. g. , a gene) can increase power for identifying variants with a common direction of effect. We report comparison of whole genome sequences from HIV-1 exposure extremes—highly-HIV-exposed individuals who remained HIV-uninfected and lower-exposed individuals who became HIV-infected. We discover and replicate associations between HIV-1 risk and aggregate variation in two genes, CD101 and UBE2V1 that increase directly with the level of HIV-1 exposure. These genes may modulate host inflammation thereby identifying molecular mechanisms linking genital inflammation to HIV-1 infection, possibly leading to novel candidate host-directed HIV-1 prevention interventions.
blood cells
genome-wide association studies
innate immune system
medicine and health sciences
immune physiology
cytokines
pathology and laboratory medicine
body fluids
immune cells
pathogens
immunology
microbiology
alleles
retroviruses
viruses
immunodeficiency viruses
ethnicities
developmental biology
rna viruses
genome analysis
molecular development
white blood cells
genomics
animal cells
medical microbiology
hiv
microbial pathogens
t cells
hiv-1
viral replication
genetic loci
blood plasma
immune system
people and places
blood
anatomy
cell biology
virology
viral pathogens
physiology
genetics
biology and life sciences
population groupings
cellular types
computational biology
lentivirus
regulatory t cells
organisms
human genetics
Whole genome sequencing of extreme phenotypes identifies variants in CD101 and UBE2V1 associated with increased risk of sexually acquired HIV-1