Transcriptome analysis of primary monocytes from HIV-positive patients with differential responses to antiretroviral therapy

Background Despite the significant contributions of monocytes to HIV persistence, the HIV-monocyte interaction remains elusive. For patients on antiretroviral therapy, previous studies observed a virological suppression rate of >70% and suggested complete viral suppression as the primary goal. Although some studies have reported genetic dysregulations associated with HIV disease progression, research on ex vivo-derived monocytic transcriptomes from HIV+ patients with differential responses to therapy is limited. This study investigated the monocytic transcriptome distinctions between patients with sustained virus suppression and those with virological failure during highly active antiretroviral therapy (HAART). Methods Genome-wide transcriptomes of primary monocytes from five HIV+ patients on HAART who sustainably controlled HIV to below detection level (BDL), five HIV+ patients on HAART who consecutively experienced viremia, and four healthy HIV sero-negative controls were analyzed using Illumina microarray. Pairwise comparisons were performed to identify differentially expressed genes followed by quantitative PCR validation. Gene set enrichment analysis was used to check the consistency of our dataset with previous studies, as well as to detect the global dysregulations of the biological pathways in monocytes between viremic patients and BDLs. Results Pairwise comparisons including viremic patients versus controls, BDL versus controls, and viremic patients versus BDLs identified 473, 76, and 59 differentially expressed genes (fold change > 2 and FDR < 0.05), respectively. The reliability of our dataset was confirmed by gene set enrichment analysis showing that 6 out of 10 published gene lists were significantly enriched (FDR < 0.01) in at least one of the three pairwise comparisons. In the comparison of viremic patients versus BDLs, gene set enrichment analysis revealed that the pathways characterizing the primary functions of monocytes including antigen processing and presentation, FcγR mediated phagocytosis, and chemokine signaling were significantly up-regulated in viremic patients. Conclusions This study revealed the first transcriptome distinctions in monocytes between viremic patients and BDLs on HAART. Our results reflected the outcome balanced between the subversion of the monocyte transcriptome by HIV and the compensatory effect adapted by host cells. The up-regulation of antigen presentation pathway in viremic patients particularly highlighted the role of the interface between innate and adaptive immunity in HIV disease progression.


Introduction
Monocytes, a key cell type in innate and adaptive immunity, have a propensity to differentiate into macrophages or dendritic cells [1,2]. This differentiation ability, along with the activities of antigen presentation, migration, chemotaxis, and phagocytosis [3], enables them to play crucial roles in HIV pathogenesis. Though less permissive to HIV infection than T cells and macrophages [4,5], monocytes could be infected by HIV [6,7] and infectious virus can be isolated from circulating monocytes in untreated patients and highly active antiretroviral therapy (HAART) responders [4,8]. By harboring and trafficking HIV into various tissue compartments through differentiating into tissue macrophages or dendritic cells, monocytes serve as important viral reservoirs [9,10]. During therapy, monocytes can maintain HIV replication throughout HAART as antiretroviral drugs may not block viral replication in monocytes as efficiently as they do in CD4+ T cells [10]. Moreover, undifferentiated monocytic precursor cells (such as CD34+ progenitor cells) infected by HIV may pass the virus to progeny monocytes and keep on renewing the viral pool in peripheral blood monocytes [11,12].
Despite the significant contributions of monocytes to HIV persistence, the underlying pathogenic mechanism is not fully understood. To better understand HIV pathogenesis at the genomic level, genome-wide transcriptomic studies of monocytes and monocyte-derived macrophages (MDM) have been carried out. For example, studies using monocytes/MDM infected by HIV in vitro revealed the key areas of monocyte dysfunctions related to inflammation [13], cytokine networks [14], cell cycle [15], cytoskeleton [16], and signaling pathways [17,18]. Other studies using ex vivo-derived monocytes identified an anti-apoptosis gene signature in viremic patients [19], a mixed phenotype with both increased and decreased pro-inflammatory features in patients with high viral load [20,21], and a novel candidate gene NAMPT correlating with the viral load in therapy-naïve patients [22]. Recently, by comparing the monocyte transcriptomes from HIV+ progressors and therapy naïve non-progressors, we have shown the systematic alteration of the interrelated pathways such as Toll-like receptor (TLR) signaling and cytokine-cytokine receptor interaction in viremic patients [23]. Although these studies have provided large datasets to facilitate our understanding, current knowledge on the dysregulations of monocytic transcriptome during HIV disease progression remains far from complete. In particular, none of the previous studies has looked into the global dysregulations of the biological pathways in monocytes from patients with sustained virus suppression versus those with virological failure during HAART, as we previously did on T cell subsets [24]. Regarding the virological suppression rate by HAART, the study on the patients receiving HAART for 12 months in Nigeria has observed a virological suppression rate of 76.7% versus a virological failure rate of 23.4% [25], whereas another study on the UK cohort has reported that 73.5% of the patients initiating HAART achieved complete virological suppression within 6 months [26].
In order to get a better insight into the dysregulations of monocytic transcriptomes from HIV+ patients with differential responses to antiretroviral therapy, this study analyzed transcriptomes of primary circulating monocytes from HIV+ patients on HAART who sustainably controlled HIV to below detection level (BDL), HIV+ patients on HAART who consecutively experienced viremia (VIR), and 4 sero-negative controls (CTR; Table 1) using Illumina HumanHT-12 v3 Expression BeadChip. Our analysis at the gene set level has shown that in the comparison of VIR versus BDL, the pathways characterizing the primary functions of monocytes including antigen processing and presentation, FcγR mediated phagocytosis, and chemokine signaling were significantly up-regulated in the VIR group. These results reflected the outcome balanced between the subversion of monocyte transcriptome by HIV and the compensatory effect adapted by host cells during disease progression following therapy. In particular, the up-regulation of antigen presentation pathway in the VIR group highlighted the role of the interface between innate and adaptive immunity in HIV disease progression.

Cluster analysis and identification of differentially expressed genes
Genome-wide transcriptomes of primary monocytes from 5 HIV+ patients on HAART who consecutively experienced viremia, 5 HIV+ patients on HAART who sustainably controlled HIV to below detection level, and 4 healthy HIV sero-negative controls (5 VIR, 5 BDL, and 4 CTR; Table 1) were analyzed using Illumina microarray. The hierarchical clustering analysis revealed that the VIR group formed an independent cluster from the BDL group and these HIV+ groups further combined into a distinct cluster from the CTR group ( Figure 1). Pairwise comparisons between the three groups were carried out and differentially expressed genes (DEGs) with FDR < 0.05 and fold change > 2 were identified for each comparison. For the comparison of VIR versus CTR, 473 DEGs were identified (324 up-regulated and 149 down-regulated; Additional file 1). For the comparison of BDL versus CTR, 76 DEGs were found (45 up-regulated and 31 downregulated; Additional file 1). When the VIR group was compared to the BDL group, 59 DEGs were detected (48 up-regulated and 11 down-regulated; Table 2). These 59 DEGs were uploaded to DAVID (Database for Annotation, Visualization, and Integrated Discovery) for the detection of DEGs overlapping with the genes in HIV interaction database [23]. Fourteen DEGs were present in HIV interaction database, covering 24% of our list (Table 2), which gave the initial confirmation of the reliability of our dataset at the discrete gene level (Additional file 2).

qPCR validation and the comparison to previous microarray studies
To further confirm the DEGs from microarray analysis, mRNA expression levels of the selected DEGs were measured by quantitative PCR (qPCR; Table 3). The DEGs were selected based on the coverage of different levels and directions of fold change, different ranges of FDR values, and/or biological significance. The cohort for qPCR validation consisted of 10 viremic patients, 10 BDLs, and 9 healthy controls (including all the original samples used in the microarray). The fold changes for each pairwise comparison evaluated by qPCR were fully consistent with the results obtained from microarray, which confirmed the reliability of our microarray data.
We then compared our dataset with published DEG lists derived from the studies on monocyte/MDM transcriptomes modulated by HIV since 2002 [15,16,19,21,22,[27][28][29][30]. The GSEA showed that in our transcriptome dataset, 6 out of the 10 published gene lists were significantly enriched (FDR < 0.01, the most stringent cut off ) in at least one of the 3 pairwise comparisons, whereas the remaining 4 gene lists reached the relaxed significance level (FDR < 0.25) in all of the 3 pairwise comparisons (Table 4). Nevertheless, the highly significant enrichment of the majority of these gene lists (FDR < 0.01) demonstrated general consistency of our data with previous studies.

Pathway analysis by GSEA
As we aimed to identify monocyte dysfunction in relation to HIV disease progression, our subsequent pathway analysis focused on the comparison between VIR   and BDL groups. GSEA was attempted using the KEGG pathways including 186 gene sets. Unlike DEG approach, GSEA uses the whole gene expression dataset to identify enriched pathways and is better able to detect small coordinated changes in gene expression in the context of gene set [31,32].
In the comparison of VIR versus BDL, 26 pathways were significantly up-regulated in the VIR group (FDR < 0.05/stringent cut off; Additional file 3), whereas no pathway down-regulated in the VIR group passed this threshold. These 26 significantly up-regulated pathways can be grouped into 3 functional categories including immune-related pathways (n = 10), disease (n = 10), and metabolism pathways (n = 6; Table 5; Additional file 3).
We chose to focus on the immune-related pathways as they have the most direct relevance to the immune dysfunction of monocytes.
Out of the ten immune-related pathways (Table 5), two were involved in cell differentiation and development (hematopoietic cell lineage and neurotrophin signaling), two were associated with transendothelial migration (cell adhesion molecules and adherens junctions), and the remaining six covered all the major aspects of innate immunity, including chemokine signaling, IgA production, complement cascade, phagocytosis, lysosome, and antigen presentation. All the core enrichment genes contributing to the up-regulation of the immune-related pathways are listed in Table 5. Since the pathways of   antigen presentation, phagocytosis, and chemokine signaling characterize the primary functions of monocytes, we further inspected their core enrichment genes.

Up-regulation of antigen presentation pathway in VIR versus BDL
The antigen processing and presentation pathway (HSA04612) was significantly up-regulated in the VIR group compared to the BDL group (FDR = 0.016). The enrichment plot and the heat map of the core enrichment genes for antigen presentation pathway are displayed as representatives to illustrate the GSEA output in Figure 2. Figure 2B shows not only the coordinated up-regulation of these core enrichment genes in the VIR group as a combination from all the viremic patients, but also the variations in gene expression between subjects within each group. For example, while VIR2 had lower expression of the core enrichment genes than VIR1, VIR2 still exhibited higher expression than the majority of the patients from the BDL group. Nineteen out of 88 gene members of this pathway were core enrichment genes (Table 5

Promoter motif analysis by GSEA
Promoter motif analysis by GSEA used gene sets that contained gene members sharing the same transcription factor binding site defined in the TRANSFAC database [33]. This analysis can identify the coordinated changes of the genes under the control of a certain transcriptional regulator. The GSEA revealed that 37 and 42 gene sets were significantly up-and down-regulated in the VIR group (FDR < 0.05) compared to the BDL group (Additional file 4). While the majority of the significantly altered gene sets contained regulatory motifs matching for the annotated transcription factors, 12 gene sets contained motifs which did not match any known transcription factor. The most significantly up-and downregulated gene sets in the VIR group (both FDR = 0.03) had gene members sharing the regulatory motif for transcription factors NFκB and MYC, respectively.

Discussion
Our study has provided the first snapshot into the transcriptome distinctions of the primary monocytes between HIV+ patients on HAART who consecutively experienced viremia and HIV+ patients on HAART who sustainably controlled HIV to below detection level ( Table 1). The main objective of the study was to identify genomic signatures associated with HIV disease progression. The relevance of the identified DEGs was initially confirmed by querying HIV interaction database using DAVID (Table 2; Additional file 2). The comparison between our dataset . Bottom, plot of the ranked list of all genes. Y axis, value of the ranking metric; X axis, the rank for all genes. Genes whose expression levels are most closely associated with the VIR or BDL group get the highest metric scores with positive or negative sign, and are located at the left or right edge of the list. Middle, the location of genes from the antigen processing and presentation pathway within the ranked list. Top, the running enrichment score for the gene set as the analysis walks along the ranked list. The score at the peak of the plot is the enrichment score (ES) for this gene set and those genes appear before or at the peak are defined as core enrichment genes for this gene set. B. Heat map of the core enrichment genes corresponding to A. The genes that contribute most to the ES, i. e., genes that appear in the ranked list before or at the peak point of ES, are defined as core enrichment genes. Rows, genes; columns, samples.

Range of colors (red to blue) shows the range of expression values (high to low).
and previous microarray studies on monocyte/MDM transcriptomes further confirmed the reliability of our data ( Table 4). The previously identified gene sets, which were highly significant (FDR < 0.01) in our dataset, reflected the altered biological functions including cytokine networks, cell cycle, signaling pathways, metabolism, immune responses, and transcriptional regulation. These biological themes were also represented by the gene sets which only reached the relaxed significance threshold (FDR < 0.25). The failure of these gene sets to achieve higher statistical cut-offs could possibly be attributed to certain biological and/or technical variations across multiple microarray studies. For example, differences in the cohort such as age, gender, and ethnicity could contribute to the biological variations, whereas microarray inter-platform differences could contribute to the technical variations. Compared to DEG approach which focused on discrete genes, GSEA enabled a more comprehensive detection of genes contributing to the enrichment of the pathways correlated with disease progression. Therefore, the subsequent discussion focused on the pathways significantly associated with progressive phase of HIV disease revealed by GSEA. Since the immune-related pathways have the most direct relevance to HIV disease and this aspect has been under constant investigations, our subsequent discussions shall center on the immune-related pathways including antigen presentation, phagocytosis, and chemokine signaling. At the transcriptome level, we observed an overall upregulation of antigen processing and presentation pathway in the VIR group compared to the BDL group (FDR = 0.016; Figure 3; Table 5), which was manifested by the upregulation of both endogenous (MHC-I) and exogenous (MHC-II) signaling branches. Five molecules from MHC-II pathway were detected as core enrichment genes including IFI30, CD74, HLA-DOA, HLA-DRB1, and HLA-DRB5 (Table 5). The coordinated up-regulation of HLA-DR was consistent with the previous studies reporting higher expression levels of surface HLA-DR on monocytes from HIV+ individuals compared to sero-negative controls [34][35][36]. As for the MHC-I pathway, the majority of the core enrichment genes were reported to have interactions with HIV, such as HSP90AB1 [37] and PA28 (PSME1) [38]. Moreover, 3 MHC-I molecules including HLA-C, HLA-E, and HLA-F were also detected as core enrichment genes with the coordinated up-regulation in the viremic patients. The importance of HLA-C was highlighted in a very recent study showing that increased HLA-C expression was correlated with increased likelihood of cytotoxic T lymphocyte (CTL) responses and frequency of viral escape mutation [39]. In addition, the substantially increased expression of MHC-I including HLA-C has been reported in monocytes in HIV progressors compared to HAART-suppressed patients, which was suggested to contribute to the generation of the dysfunctional naive CD8-low T cells that emerged during disease progression [40]. Taken together, our observation that the up-regulation of HLA-C expression on monocytes was associated with plasma viremia could possibly reflect a compensatory effect imposed by impaired CTL responses. However, this compensatory effect may be subverted by HIV, resulting in further CD8+ T cell dysfunction as proposed by Favre et al. [40].
As antigen presentation pathway is right at the interface between innate and adaptive immunity, its significant upregulation in the viremic patients provides direct evidence that adaptive immune response is perturbed through the components of innate immunity during HIV disease progression. This interference is systematic throughout the process of antigen presentation as demonstrated by the core enrichment genes covering not only the aforementioned MHC molecules, but also the genes associated with antigen digestion, loading and transportation (Table 5). In fact, the dysregulation of antigen presentation pathway is not the only pathway altered at the interface. Other innate immunity pathways crucial for the development of adaptive immunity have also been associated with HIV disease progression by this and previous studies [23,24,41]. For example, this study detected that complement cascade regulating both B and T cell responses [42] was significantly up-regulated in the viremic patients (Table 5). Consistently, complement pathway was also found to be up-regulated in the viremic patients versus BDLs by our previous study on primary CD4+ T cells [24]. In addition, a recent study has reported that monocytes and complement system contributed to the tuberculosis-associated immune reconstitution inflammatory syndrome in HIV-TB co-infected patients [41]. Altogether, these studies have demonstrated the importance of complement components in HIV disease progression. Recently, another key component of innate immunity, Toll-like receptor signaling (TLR) pathway, has been found to be significantly down-regulated in monocytes from viremic patients versus long-term non-progressors [23]. The different direction of the changes may reflect the various aspects of HIV-host interactions that contribute to disease progression, such as HIV persistence and impairment of T cell functions. Despite this difference, all of the aforementioned studies point towards the adaptive immunity being perturbed at the interface where innate and adaptive immunity interact during HIV disease progression.
Consistent with the previous reports on phagocytosis dysfunction in monocytes upon HIV infection [43][44][45], we observed the significant up-regulation of FcγR-mediated phagocytosis pathway (HSA04666) in the VIR group versus the BDL group (FDR = 0.029; Figure 4). The expression of FCGR1A (FcγRI/CD64), the sole high affinity receptor for monomeric IgG, was coordinately increased along with other core enrichment genes in the viremic patients. This observation was consistent with the previous study reporting the increased expression of FcγRI on monocytes from acute HIV infection patients compared to chronically-infected individuals and healthy controls [46]. In addition, the previous study suggested that monocytes might expand and/or up-regulate the expression of this high affinity FcγR in response to the burst of viral replication. This host response may also account for the increased FCGR1A expression in the viremic patients on HAART, as detected in this study.
The perturbation of phagocytosis pathway was further manifested by the up-regulation of the genes encoding for kinases in early signaling events such as SRC kinases (HCK and LYN) and SYK, as well as the genes encoding for proteins involved in oxidative burst (SPHK1/NADPH oxidase) and cytoskeleton rearrangement ( Table 5). The majority of the core enrichment genes, which spread along the signaling branches of PI3K -AKT, PIP5K -VASP/ WASP/WAVE -ARP2/3, and VAV -RAC -PAK1 -CFL1, have also been implicated previously in HIV pathogenesis. A recent study using gene knockdown and pharmacological inhibition of SRC and AKT identified them as mediators of HIV-induced inhibition of autophagocytosis in bystander macrophages/monocytic cells [47]. The elevated AKT activity was also detected in HIVinfected macrophages, and PI3K/AKT inhibitors were suggested as a novel therapy for interfering with the establishment of HIV reservoirs [48]. In the signaling branch pointing towards ARP2/3, WAVE2 was involved in the activation of the actin polymerization nucleator ARP2/3 complex, which played a role in the migration of the viral core components toward the cell nucleus and the efficient infection of HIV in cell lines [49,50]. In addition, HIV Nef was found to activate the signaling branch of VAV -RAC -PAK (p21), which markedly enhanced the NADPH oxidase response [51]. Nef also mediated CFL1 phosphorylation, which was crucial for maintaining actin homeostasis [52]. Taken together, the up-regulated expression of the genes involved in various kinase cascades and the gene encoding for FcγRI may reflect the HIV interference with the host genetic machinery on the one hand, and the compensatory processes adapted by the host on the other hand during virological failure.
Chemokine signaling pathway (HSA04062) was also significantly up-regulated in the VIR group versus the BDL group (FDR = 0.038), which was manifested by the systematic up-regulation of the genes encoding for chemokine ligand 2 (CCL2/MCP-1) and chemokine receptors including CCR1, CCR2, CCR3, CCR7, CCR9, CX3CR1, CXCR5, and XCR1 (Table 5; Figure 5). The up-regulation of MCP-1 in viremic patients was also reported by Pulliam et al. [20], which was suggested to be involved in enhancing HIV production and spread [53]. In addition, it was shown that the sequential activation of SRC, MAPKs, and PI3K -AKT -NFκB pathways resulted in the increased expression of MCP-1 [54]. Given the fact that the majority of these genes overlapped with the core enrichment genes we detected, it is thus plausible to hypothesize that the up-regulation of MAPK and NFκB braches signaling through Gα subunit may be related to the coordinated up-regulation of MCP-1, which could contribute to HIV spread. This potential link could be one of the mechanisms underlying the association of the up-regulation of chemokine signaling pathway with virological failure during HAART.
In the comparison of VIR versus BDL, the promoter motif analysis yielded a list of significantly up-and downregulated gene sets with members of each gene set containing the same binding site for a certain transcription factor (Additional file 4). The transcription factors implicated by these significant gene sets included some of the well-known regulators, such as SP (SP1 and SP3), NFκB (RELA binding motif), AP1, AP2, and CREB1, which have been demonstrated to play crucial roles in HIV transcription regulation [17]. Interestingly, the most significantly up-regulated gene set in the VIR group (FDR = 0.03) contained the binding motif for RELA, which was also identified as the core enrichment gene in chemokine signaling pathway in the VIR group. This overlap not only demonstrated the consistency between pathway and promoter motif analysis, but also implicated the complexity of the transcriptional regulation underlying HIV-monocyte interaction. In addition to the aforementioned transcription factors, some of the significant gene sets contained the binding motifs without matching any known transcription factor, which may suggest potential novel regulators that warrant future investigations.
A few limitations of this study should be noted. First, this study used a cross-sectional design, which could not provide dynamic findings from a longitudinal perspective. Secondly, this study used a relatively small sample size appropriate for the pilot investigation and future studies using larger sample size are thus warranted to further confirm the results. Finally, although the GSEA identified a panel of significantly altered gene sets with high relevance to HIV disease progression during therapy, the predictive nature of GSEA as the common limitation of statistical tools should be noted and interpretations should be made with caution. To overcome this limitation, future biological experiments should be conducted to directly confirm and further explore these findings.

Conclusions
This study has revealed the first transcriptome distinctions in monocytes between HIV+ patients on HAART who consecutively experienced plasma viremia and HIV+ patients on HAART who sustainably controlled plasma viremia to below detection level. Compared to the BDL group, the pathways characterizing the primary functions of monocytes including antigen processing and presentation, FcγR mediated phagocytosis, and chemokine signaling, were significantly up-regulated in the VIR group. Our results reflected the outcome balanced between the subversion of monocyte transcriptome by HIV and the compensatory effects adapted by host cells. Furthermore, the altered pathways of antigen presentation and complement cascade highlighted that HIV manipulated adaptive immune response via innate immunity components at the interface of innate-adaptive immunity during disease progression. These data offered new comparative insights into the perturbed genetic networks of ex vivo-derived monocytes subverted by HIV during disease progression on therapy. In-depth functional studies on the regulation of these pathways and the corresponding core enrichment genes along with proteomic analysis may further confirm our findings and provide detailed molecular mechanisms underlying HIV-monocyte interaction.

Patient profiles and collection protocol
Five HIV+ patients on HAART who sustainably controlled HIV to below detection level (viral load < 40 copies of HIV RNA/ml), five HIV+ patients on HAART who consecutively experienced viremia, and four healthy HIV sero-negative controls were studied (Table 1). Patients with viral load < 1,000 were chosen to represent patients on HAART with virological failure because the previous study has demonstrated that the transcriptome profiling of these patients was more homogenous [23]. It has also been shown that monocyte transcriptomes from patients with viral load < 1,000 and those with higher viral load exhibited similar dysregulated pathways during disease progression. Patients in the VIR and BDL groups were on HAART for > 52 weeks. The BDL and VIR groups had a broad range of CD4+ T cell counts (100-796 cells/μl) as our objective was to find unbiased differences of gene expression profiling of monocytes between patients on HAART with sustained virus suppression and virological failure regardless of the degree of T cell decline. This grouping criterion by viral load has been successfully used by the previous studies [20,23,24]. These patients received two NRTIs (zidovudine, lamivudine, stavudine, emtricitabine, tenofovir) in association with one or two protease inhibitors (darunavir, ritonavir, indinavir, saquinavir, atazanavir). Ten patients were from the HIV clinic at Westmead Hospital and the four healthy controls were from the Australian Red Cross Blood Service in Sydney. This study was approved by the Sydney West Area Health Services Research Ethics Committee, and all blood samples were collected after individual informed written consent.

Purification of CD14+ monocytes and RNA isolation
A single blood sample (10-20 ml in EDTA) was obtained from each patient. After separation of plasma, primary PBMCs were isolated immediately after obtaining blood samples by Ficoll-gradient centrifugation and purified. This aspect was strictly followed in our experiments because of previously described lower RNA yields and possible changes in gene expression profiles upon storage of blood [55]. CD14+ monocytes were then obtained by positive isolation with antibody-conjugated magnetic beads according to the manufacturer's instructions (Miltenyi Biotech, Germany) with a purity > 98.6% as verified by flow cytometry. Binding of antibody to CD14 does not trigger signal transduction and a previous study has clearly demonstrated that CD14 positive selection does not alter cellular transcriptome by comparing gene expression profiles in parallel using either positive or negative selection [55]. Total RNA was isolated from purified cells using RNeasy Mini kit (Qiagen Pty Ltd., Clifton Hill, Victoria, Australia) with an integrated step of on-column DNase treatment.
cRNA preparation, microarray hybridization and scanning RNA quality was checked by Agilent Bioanalyzer and RNA Integrity Scores were higher than 9 for all the samples (Table 1). cRNA amplification and labeling with biotin were performed using Illumina TotalPrep RNA amplification kit (Ambion, Inc., Austin, USA) with 250 ng total RNA as input material. cRNA yields were quantified with Agilent Bioanalyzer and 750 ng cRNAs were hybridized to Illumina HumanHT-12 v3 Expression BeadChips (Illumina, Inc., San Diego, USA). Each chip contains 12 arrays and each array contains >48,000 gene transcripts, of which, 46,000 are derived from human genes in the National Center for Biotechnology Information (NCBI) Reference Sequence (RefSeq) and UniGene databases. All reagents and equipment used for hybridization were purchased from Illumina, Inc. According to the manufacturer's protocol, cRNAs were hybridized to arrays for 16 hours at 58°C before being washed and stained with streptavidin-Cy3. Then the beadchips were centrifuged to dry and scanned on the Illumina BeadArray Reader confocal scanner. To minimize the batch effect, the microarray chips were all processed at the single site using the same platform with the identical setting of the parameters by the same experimenter. The microarray dataset has been submitted to GEO (Accession Number GSE52900).

Analysis of differentially expressed genes
The quality of the entire data set was assessed by box plot and density plot of bead intensities, density plot of coefficient of variance, pairwise MAplot, pairwise plot with microarray correlation, cluster dendrogram, and nonmetric multidimensional scaling using R/Bioconductor and the lumi package [56]. Based on the quality assessment, all 14 samples were deemed suitable for further analysis. Data normalization was performed using log2 transform and a robust spline normalization (RSN) implemented in the lumi package for R/Bioconductor [56,57]. Cluster analysis of gene expression profiling was carried out using dist and hclust functions from R stats package. Euclidean distance and complete linkage were used for distance metric and linkage criterion, respectively. To reduce false positives of differentially expressed genes, genes below detectable limit (based on a detection p value cutoff 0.01) were removed from the dataset. A linear model fit in conjunction with an empirical Bayes statistics was used to identify candidate DEGs [58]. P values were corrected for multiple testing using FDR adjustment implemented in lumi package. Pairwise comparisons for the 3 groups were carried out and candidate DEGs with fold change >2 and FDR <0.05 were identified for each of the comparisons. The list of the DEGs derived from the comparison of VIR versus BDL was uploaded to DAVID for the detection of the DEGs showing overlap with the genes in HIV interaction database at NCBI [59,60].
Gene set enrichment analysis GSEA was used for the comparison of our dataset with the published DEG lists from the previous studies (in vivo and ex vivo since 2002) [15,16,19,21,22,[27][28][29][30], the investigation of global dysregulations of the biological pathways, and the promoter motif analysis. For the comparison with previous studies, 10 DEG lists were used from the studies on monocyte/MDM transcriptomes modulated by HIV (Table 4). For the pathway investigation, the gene sets were from MsigDB [32], catalog C2 functional sets, subcatalog KEGG pathways, which included 186 gene sets from pathway databases. For the promoter motif analysis, C3 motif gene sets which contained gene members sharing the same transcription factor binding site were used. This collection included 615 gene sets and each of them was annotated by a TRANSFAC record [33,61].
Instead of focusing on discrete DEGs, we analyzed the entire transcriptome data with GSEA to identify genes coordinately regulated in predefined gene sets from various biological pathways [32]. For each group comparison, GSEA was performed using the normalized data of entire 48,803 transcripts (GSEA version 2.07, Broad Institute http://www.broad.mit.edu/gsea). First, a ranked list was obtained by ranking all genes according to the correlation between their expression and the group distinction using the metric signal to noise ratio. Then the association between a given gene set and the group was measured by the non-parametric running sum statistic termed the enrichment score (ES), which was calculated by walking down the ranked list (increasing ES when encountering a gene in the given gene set and decreasing ES when encountering a gene not in the gene set). To estimate the statistical significance of the ES, a nominal p value was calculated by permuting the genes 1,000 times. To adjust for multiple hypothesis testing, the maximum ES was normalized to account for the gene set size (NES) and the false discovery rate (FDR) corresponding to each NES was also calculated. Along with the pathway enrichment results, the details report for each significant pathway was simultaneously generated. This report listed the details of each gene member in columns, one of which indicated whether this gene member was "core enrichment gene" or not. The core enrichment genes account for the enrichment signal of the pathway and the inspection of them can reveal a biologically important subset within the pathway [32].

Real-time quantitative PCR
Ten genes and 14 pairs of group comparison were selected for validation based on the coverage of different levels and directions of fold change, different ranges of FDR values, and/or biological significance. Purified total cellular RNA was used for reverse transcription with oligo d (T) and Superscript III followed by RNase H treatment (Invitrogen Life Technologies). The cDNA was then subjected to qPCR in a 96-well format in triplicate reactions with defined primers and SYBR Green (Invitrogen Life Technologies). The qPCR reactions were carried out for the extended cohort consisting of 10 viremic patients, 10 BDLs, and 9 healthy controls (including all the original samples used in the microarray) using Mx3005P™ QPCR System (Stratagene). The mean expressions of the tested genes in each group were obtained and the housekeeping gene GAPDH was used as an internal control and the normalizer for all data. The fold change was calculated by the relative quantitation method 2 -(ddCt) . Primer sequences for each transcript are available from the authors upon request.