A total of 494,656 SNPs passing quality control were tested for association with levels of HIV-1 replication in monocyte-derived macrophages (MDM) using linear regression. With data from 191 healthy blood donors whose MDM ranked in the bottom quartile with lowest (n = 95 donors) or top quartile with highest (n = 96 donors) Gag p24 production 14 days post infection with HIV-1 ^[39] (Table 1), we found strongest associations for SNPs in the genes PDE8A (rs2304418, p = 2.4×10⁻⁶ and rs12909130, p = 8.3×10⁻⁶), UBR7 (rs2905, p = 7.0×10⁻⁶), MOAP1 (rs1046099, p = 9.9×10⁻⁶ and rs1270629, p = 1.09×10⁻⁵), DYRK1A (rs12483205, p = 2.2×10⁻⁵) and SPOCK3 (rs17519417, p = 2.5×10⁻⁵) (Figure 1, Table S1). The two SNPs in PDE8A were found to be in high linkage disequilibrium (LD; r² = 0.97), whereas only a moderate degree of LD was found between rs1046099 and rs1270629 in MOAP1 (r² = 0.54) (Table 2). Table 2 shows all other SNPs associated with HIV-1 replication in MDM (cutoff p value = 5×10⁻⁵; n = 16), and includes information about LD, location, number of donors homozygous for the minor allele (MIN) and empirical p values after permutation testing. Genotyping results for none of these 16 SNPs violated Hardy-Weinberg equilibrium.

The empirical p values for linear regression using 10⁷ permutations of the genotype for each SNP in Table 2 were in close agreement with the asymptotic p values (Table 2). However, none of these SNPs remained statistically significant after a conservative correction for multiple testing (Bonferroni threshold of p<1×10⁻⁷; Figure S1). We next investigated the association between the SNPs in the five most promising genes (rs2304418 in PDE8A, rs2905 in UBR7, rs1046099 and rs1270629 in MOAP1, rs12483205 in DYRK1A and rs17519417 in SPOCK3) in a second group of blood donors (replication cohort). MDM from an independent group of 32 healthy blood donors were infected with HIV-1 and Gag p24 levels were measured 14 days after infection. One donor was found to be homozygous for the 32 base pair deletion in CCR5, rendering the cells completely resistant to infection with CCR5-using HIV-1, and results obtained with macrophages from this donor were excluded from further analysis. The associations for the SNPs in PDE8A, UBR7, MOAP1 and SPOCK3 with Gag p24 production by MDM could not be replicated in this small group of 31 donors (data not shown). However, in this replication cohort we again found a strong association between SNP rs12483205 in DYRK1A and in vitro HIV-1 replication in macrophages (p = 0.0034) (Figure 2, Table S2), also after correction for cell number and normalization (p = 0.0081; n = 28, information on cell number was missing for 3 donors) (data not shown, Table S2). This association remained significant after correction for multiple testing (Bonferroni corrected p = 0.020 (0.0034×6 SNPs for which we tried to replicate the association)) and when calculating the empirical p value for linear regression using 10⁵ permutations of the genotypes (p = 0.004) (Table S1).

The DYRK1A SNP rs12483205 genotypes were then also determined for the 202 donors with MDM that had intermediate Gag p24 levels in the supernatant ^[39]. Adding results from this group of donors did not change the strength of the association between the rs12483205 genotype and HIV-1 replication in MDM (p = 2.09×10⁻⁵; n = 393) (Figure S2A, Table S3). As expected, combined analysis of both the initial group of donors (n = 393) and the replication cohort (n = 28) revealed an increase in the strength of the association for rs12483205 in DYRK1A (p = 4.84×10⁻⁶) (Figure S2B, Table S4).

The 32 base pair deletion in CCR5 strongly affects replication of HIV-1, also in macrophages in vitro^[39]. MDM from donors that were heterozygous for this 32 base pair deletion (n = 33) had significantly lower HIV-1 replication than donors without the deletion (n = 158) (p = 0.0036, two-sample t-test) (Table S5). There were significantly more CCR5 wt/Δ32 heterozygous donors present in the rs12483205 heterozygous (HZ) and rs12483205 MIN groups (p = 0.02, Fisher Exact test). After adjusting for the CCR5 wt/Δ32 heterozygous genotype by using it as a covariate in multivariate analysis, the strength of the association between the DYRK1A rs12483205 genotype and HIV-1 replication in MDM decreased (from p = 2.2×10⁻⁵ to p = 1.2×10⁻⁴), but persisted (Table S5). Moreover, the association of rs12483205 in DYRK1A with HIV-1 replication in MDM in the replication cohort of 31 donors, was independent of the CCR5 Δ32 genotype (p = 0.0022, and p = 0.007 when corrected for cell number and normalized for date of isolation; Table S2), suggesting only a spurious relationship between the CCR5 locus on chromosome 3 and SNP rs12483205 in DYRK1A located on chromosome 21. Indeed, this was supported by results from the multivariate analysis on the total group of 393 donors (p = 3.8×10⁻⁵, Table S3) and on the combined group of 421 (393+28 donors from the replication cohort) individuals (p = 9.7×10⁻⁶, Table S4).

Next, we investigated the association between HIV-1 load or disease course in vivo and the six most promising SNPs identified in our GWAS (rs12483205 in DYRK1A, rs2304418 in PDE8A, rs2905 in UBR7, rs1046099 and rs1270629 in MOAP1 and rs17519417 in SPOCK3). Results from the Amsterdam Cohort Study (ACS) GWAS on HIV/AIDS disease progression (Van Manen et al., manuscript in preparation) and three previously published HIV-1 cohort studies with GWAS data were used: the CHAVI cohort ^[44], the GRIV ^[45], ^[46] and the MACS ^[47]. To avoid possible overlap between the replication cohorts, MACS samples were removed from the CHAVI cohort and the analysis was performed only using data from the European subset (EURO-CHAVI).

Illumina Human Hap BeadChips were used for genotyping DNA samples in the EURO-CHAVI cohort (550 K/1 M), the GRIV (300 K) and the ACS (300 K), whereas in the MACS they used Affymetrix Human 500 K GeneChips for the first stage discovery analysis ^[47]. Using different genotyping platforms results in genotype data for a different set of SNPs, which occasionally makes it necessary to estimate the unobserved genotypes, a process called SNP imputation ^[48]. When the imputation score (a measure for the reliability of the imputation) was below 0.8, imputation was considered unreliable and data was therefore not used. Results for SNPs rs2905 (UBR7), rs1046099 (MOAP1), rs1270629 (MOAP1) and rs12483205 (DYRK1A) in the MACS, and rs1046099 (MOAP1) and rs17519417 (SPOCK3) in both the GRIV and the ACS, were based on SNP imputation. The imputation score for rs12483205 in DYRK1A was 0.94, and <0.8 for rs2905, rs1046099 and rs1270629. For rs1046099 and rs17519417 the imputation scores were 0.88 and 1 respectively, in the GRIV and in the ACS.

Table 3 shows the results for association testing between disease progression or viral load and the identified SNPs in DYRK1, PDE8A, UBR7, MOAP1 and SPOCK3. No significant associations or trends were found for the SNPs in PDE8A, UBR7 or MOAP1. However, SNPs rs12483205 in DYRK1A and rs17519417 in SPOCK3 were significantly associated with in vivo endpoints in at least one of the four independent cohorts analyzed. While no association was observed in the GRIV and the ACS, we found an association between rs12483205 in DYRK1A and disease progression in the MACS (p = 0.0048, Table 3) and in the EURO-CHAVI cohort (p = 0.035, Table 3). In both cohort studies the minor allele of rs12483205 was associated with slower disease progression, which was consistent with the decreased in vitro replication of HIV-1 in MDM for the minor allele of rs12483205.

Furthermore, there was an association between SNP genotype and progression to AIDS for rs17519417 in SPOCK3 in the MACS ^[47] (p = 0.012, Table 3), with the C allele associated with slower progression. This was also in agreement with the association with lower in vitro HIV-1 replication in MDM in our study.

The DYRK1A mRNA has multiple splice forms. According to the NCBI Entrez Gene database ^[49] the alternative splicing can yield four different DYRK1A transcript variants, referred to as transcripts 1, 2, 3 and 5 (Figure 3). In addition to differences in the number or size of exons, there is also variation in the length of both the 5′ and 3′ untranslated region (UTR). Transcript variant 1 encodes the longest isoform (763 amino acids), which is a fraction larger than isoform 2 (754 amino acids). Isoforms 3 and 5 lack a C-terminal part of the full length protein (179 and 234 amino acids respectively), and consequently miss among others a poly-His domain. SNP rs12483205 in DYRK1A lies ∼900 base pairs downstream of the most distant 5′ UTR fragment of transcript variant 3 (Figure 3). PCR analysis with primers able to discriminate between all four known transcript variants (Figure S3) confirmed the presence of transcript variants 1, 2 and 3 in MDM, and in the glioma cell line U87 that was used as a positive control for the PCR (Figure 4). Transcript variant 5 was present in U87 cells, but was not detected in MDM (Figure 4). The DYRK1A gene region is devoid of known SNPs that are in high LD with rs12483205; currently available data on LD in this region does not reveal any SNP with an r²>0.5 in the Caucasian population ^[50]–^[52].

This study has been conducted in accordance with the ethical principles set out in the declaration of Helsinki, written informed consent was obtained from all participants and was approved by the Medical Ethics Committee of the Academic Medical Center in Amsterdam and the Ethics Advisory Body of the Sanquin Blood Supply Foundation, The Netherlands.

We previously determined the ability of HIV-1 to replicate in monocyte-derived macrophages (MDM) from 429 different healthy seronegative blood donors ^[39]. In brief, Gag p24 levels were measured in MDM culture supernatant 14 days post infection by an in-house enzyme-linked immunosorbent assay. To correct for differences in the number of viable MDM present at day 14 post infection, p24 levels were expressed per 10,000 cells. Since monocyte isolations were performed in four time frames and by two operators, p24 levels were normalized by dividing through the median per period and operator. These normalized p24 levels were subsequently used as a measure for in vitro HIV-1 replication in MDM.

After excluding donors that were homozygous for the 32 base pair deletion in CCR5 (n = 5) or not from European decent (n = 30), we selected 192 individuals whose MDM gave the highest (n = 96) or the lowest (n = 96) p24 production in vitro, for SNP genotyping; thus representing two groups of donors with MDM that had a more extreme phenotype (Table 1). Inclusion of donors with extreme phenotypes is known to increase power in genetic association studies ^[93]. These two groups with a more extreme phenotype were separated by a group with intermediate HIV-1 replication in MDM which consisted of 202 donors. There was no difference in age or distribution of males and females between the group with MDM that had low HIV-1 replication and the group with MDM that had high Gag p24 production (p = 0.95, Pearson chi-square test and p = 0.37, two sample t-test, respectively; Table 1). One donor from the group with low in vitro HIV replication in MDM was excluded for further analysis because the corresponding DNA samples did not pass the quality control (SNP call rate <0.98; see paragraph “Quality control of SNP genotyping data” below).

Cells (monocytes and lymphocytes) used to replicate the genome-wide association study (GWAS) findings were derived from blood samples collected from an independent group of 32 healthy blood donors (replication cohort).

Peripheral blood lymphocytes (PBL) from each donor were used for DNA isolation using the QIAamp “DNA blood mini kit” (Qiagen, Valencia, CA, USA). For the two groups of donors whose MDM had highest (n = 96 donors) or lowest (n = 96 donors) in vitro HIV Gag p24 production, we used the Illumina Infinium Human Hap 610-Quad BeadChip (Illumina, San Diego, CA, USA) ^[94] for SNP genotyping. ABI TaqMan® SNP genotyping assays were used to genotype DNA from the replication cohort for rs12483205 (C_31609775_10), rs12909130 (C_1342209_10) as a proxy for rs2304418, rs2905 (C_3236245_20), rs1046099 (C_7585751_10) and rs17519417 (C_33238869_10) (Applied Biosystems, Carlsbad, CA, USA). All TaqMan® assays were run on a LightCycler® 480 system (Roche, Basel, Switzerland) using the following amplification cycles: 10 min 95°C; 50 cycles of 15 sec 95°C, 1 min 60°C.

One donor DNA sample did not pass the minimal genotyping call rate threshold of 98% and this sample was excluded for further SNP association analysis. The Illumina 610Q SNP bead chip contained 620,901 markers to detect common genetic variation, including 21,890 markers for copy number variation (CNV). After excluding a total of 126,245 SNPs (minor allele frequency <0.05, or SNP call frequency <0.98, or SNPs for which there were no donors homozygous for the minor allele) and copy number variation markers, 494,656 SNPs were left to be tested for association between genotype and in vitro HIV-1 replication MDM. Violation of Hardy-Weinberg equilibrium was assessed post-analysis for all SNPs with p<5×10⁻⁵. Gender calls in BeadStudio (version 3.1.3.0) were compared with self-reported gender and were found to match for all donors (n = 191).

Donors who reported that either one or more parents or grand parents were born outside of Europe had been excluded from further analysis (as described in the paragraph “Study population” above). The genetic homogeneity of the genotyped donors whose MDM gave lowest (n = 96) or highest (n = 96) Gag p24 production in vitro was confirmed by both Structure ^[95] and Eigenstrat ^[96] analysis. As expected, the Q-Q plot displaying the normality of the p value distribution did not show deviations from what is expected under the null hypothesis (λ = 1.0024; λ = 1.0 is expected with a null hypothesis; Figure S4), indicating little effect of population stratification.

We set a lower detection limit (background signal) for the normalized p24 values based on the highest results that we found for a donor that was homozygous for the 32 base pair deletion in CCR5. Normalized p24 values were log10 transformed to allow for parametric testing.

Association was tested assuming an additive relationship between the genetic variants (wild-type, heterozygous and homozygous for the minor allele) and phenotype, using linear regression. Permutation testing was done to empirically determine the corresponding p value based on the observed linear regression F-statistic relative to the distribution of the F-statistic estimated from 10⁷ permutations of the genotypes. All R code and data used for the calculations can be found as supplementary online material (Tables S1, S2, S3, S4 and S5 and text file “R Scripts S1”).

Linkage disequilibrium (LD) between SNPs was calculated using data obtained from genotyping DNA from our studied population (n = 191) with the Illumina SNP chip, using Haploview software (version 4.2) ^[50]. SNP genotype imputations for our study population were performed using Impute software ^[48] and the Caucasian HapMap population (release 21) as the reference panel.

Seven day old MDM and U87 cells were used for RNA isolation and subsequent cDNA synthesis using the Roche “high pure RNA isolation” and “transcriptor first strand cDNA synthesis” kit (Roche, Basel, Switzerland). To uniquely detect the presence of transcript variants 1, 2, 3 and 5, we used primer pairs [gene: 5′-TGATATTGTCATGTTACAGAGGCGG-3′] (forward) and [gene: 5′-CTGACGCACCTGGGGACTG-3′] (reverse), [gene: 5′-TGTCTCTGAGGTTCTTTTCCAGTG-3′] (forward) and [gene: 5′-AGCACCCTCTCAATTCCCAATGCC-3′] (reverse), [gene: 5′-GCTCGCACGTGGTTCATTTGCT-3′] (forward) and [gene: 5′-TCCTTAGACAGGAACGTCATGAACCT-3′] (reverse), and [gene: 5′-CAGGAGGACCTGGTGGGCGA-3′] (forward) and [gene: 5′-TGCTGACGCACCTGAGCTTG-3′] (reverse) respectively (Figure S3). For the detection of transcript variants 1, 2, 3 and 5 the following amplification cycles were used: 2 min 95°C; 35 cycles of 30 sec 95°C, 30 sec 62°C (64°C for transcript 1), 30 sec 72°C (90 sec for transcript 1); 10 min 72°C.