Differential Adipose Tissue Gene Expression Profiles in Abacavir Treated Patients That May Contribute to the Understanding of Cardiovascular Risk: A Microarray Study

Objective To compare changes in gene expression by microarray from subcutaneous adipose tissue from HIV treatment naïve patients treated with efavirenz based regimens containing abacavir (ABC), tenofovir (TDF) or zidovidine (AZT). Design Subcutaneous fat biopsies were obtained before, at 6- and 18–24-months after treatment, and from HIV negative controls. Groups were age, ethnicity, weight, biochemical profile, and pre-treatment CD4 count matched. Microarray data was generated using the Agilent Whole Human Genome Microarray. Identification of differentially expressed genes and genomic response pathways was performed using limma and gene set enrichment analysis. Results There were significant divergences between ABC and the other two groups 6 months after treatment in genes controlling cell adhesion and environmental information processing, with some convergence at 18–24 months. Compared to controls the ABC group, but not AZT or TDF showed enrichment of genes controlling adherence junction, at 6 months and 18–24 months (adjusted p<0.05) and focal adhesions and tight junction at 6 months (p<0.5). Genes controlling leukocyte transendothelial migration (p<0.05) and ECM-receptor interactions (p = 0.04) were over-expressed in ABC compared to TDF and AZT at 6 months but not at 18–24 months. Enrichment of pathways and individual genes controlling cell adhesion and environmental information processing were specifically dysregulated in the ABC group in comparison with other treatments. There was little difference between AZT and TDF. Conclusion After initiating treatment, there is divergence in the expression of genes controlling cell adhesion and environmental information processing between ABC and both TDF and AZT in subcutaneous adipose tissue. If similar changes are also taking place in other tissues including the coronary vasculature they may contribute to the increased risk of cardiovascular events reported in patients recently started on abacavir-containing regimens.


Introduction
Antiviral therapies for HIV infection have been associated with abnormalities in serum lipids, [1] an increased incidence of type 2 diabetes mellitus, [2] and increased cardiovascular disease. [3][4][5][6][7] However the metabolic and adipose tissue effects of antiretroviral regimens are not the same with all regimens. [8] Regimens containing tenofovir disoproxil fumarate (TDF) have demonstrated less effect on serum lipids and subcutaneous adipose tissue than those containing zidovidine (AZT). [9] Moreover treatment with abacavir-containing regimens has been associated with increased risk of cardiovascular disease in some studies [10][11][12][13][14] but not others [15][16][17][18] We have previously reported changes in the expression of genes controlling lipid metabolism and mitochondrial respiratory chain in biopsies taken from patients 6 and 18-24 months after starting antiretroviral regimens containing different nucleoside analogue regimens combined with efavirenz. [19,20] Here we report results from a microarray study performed on subcutaneous adipose tissue biopsies from the same study.

Patients and Methods
Details of the patients have been previously reported and the design of the study is summarised in Fig. 1. [19,20] Briefly, the study was undertaken in two parts. In the first part iliac crest fat biopsies were performed on 32 HIV-positive patients naïve to antiretroviral therapy who were randomised to receive either zidovidine (AZT)/lamivudine (3TC n = 15) in fixed dose preparation or tenofovir (TDF)/emitricitabine (n = 17), again in fixed dose preparation. Both groups also took efavirenz (EFV). All patients underwent extensive biochemical screening, tests of body fat distribution (whole body DEXA and abdominal CT scans) and iliac crest biopsy under local anaesthesia, as previously described. [19] Tests were performed before starting treatment, at six months and again at 18-24 months after initiating therapy. 15 HIV negative controls underwent similar investigations.
In the second non-randomised part of the study patients who were on their initial antiretroviral regimen containing abacavir (ABC) and 3TC in fixed-drug combination plus EFV for 6 months (n = 8) or 18-24 months (n = 11) were also tested in a similar way. Patients and controls from both parts of the study were matched for age, ethnicity, gender, weight, and BMI. None had ever had an AIDS defining disease. All three HIV groups were also matched with regards to pre-treatment CD4 count, viral load [19,20] and family history of diabetes and cardiovascular disease. The median age (range) for the controls was 37 , for AZT 33 , for TDF 35.5 (23-53) and for ABC 37 . The median (IRQ) CD4 count before treatment were AZT 187 (153-268), TDF 234 (163-253) and ABC 215 (75-306). Pre-treatment viral loads log 10 (IRQ) were 5.1 (4.4-5.7), 4.6 (4.2-5.3) and 4.5 (4.4-4.9) respectively for AZT, TDF and ABC. All but two patients (both on AZT) had undetectable viral load when tested at Competing Interests: Funds obtained from "Gilead Sciences", "Pharmacia Upjohn" and "Dupont Pharmaceuticals" were provided as unrestricted educational grants. The above companies have no say in the design of the study, played no role in its conduct or the writing of the study. None of the companies have any ownership of the results, nor have any of the authors accepted consultancies, have products in development, patents or any other financial connections with the above companies. At the time of the microarray analysis one of the authors (KP) was employed by Beckman Coulter Genomics where he performed the initial statistics. Beckman Coulter Genomics was paid to perform the microarray and the bio-statistics. Since then KP has moved from the company and now works for ILS Genomics. He has continued to perform the biostatistics in his private time and entirely as a co-contributor to the study. Neither company has any claims to the data, control over the publication or ownership of any of the material. This does not alter the authors' adherence to PLOS ONE policies on sharing data and materials. 6 month following treatment (209 and 109,638 copies/ml) but all had undetectable viral loads at 18-24 m. Both admitted to poor compliance. None developed any clinical evidence of lipodystrophy during the period of study and post treatment weight, BMI, DEXA and abdominal CT results were not significantly different between groups at both time points after treatment. Serum lipids, fasting glucose and other parameters were similar between groups at each time point. [19,20] Samples for microarray were available from 103 biopsy specimens (Fig. 1). There was no significant difference in weight, BMI, central and peripheral fat, and visceral and subcutaneous fat, or blood biochemistry in the patients undergoing microarray study at any time point.

Microarray methods
RNA was extracted for microarray and array data generated using a Qiagen RNeasy Lipid Tissue Mini Kit. Total RNA was converted into labeled cDNA using the WT-Ovation Pico RNA Amplification System and samples were hybridized to an Agilent Whole Human Genome Microarray 4×44K (AMADID 014850, Agilent Technologies). The Agilent arrays were scanned with the Agilent Technologies DNA Microarray Scanner G2505C with default settings. Scanned images were extracted and initial quality control performed in Agilent Feature Extraction (AFE) software version 10.7.3.1. Additional quality assessment was performed using the R arrayQualityMetrics package. [21] As no major experimental problems could be detected, all microarrays were included in subsequent analysis steps. The dataset was read into R and processed using the Agi4x44PreProcess package with options to use the AFE Processed Signal as foreground signal and BG Used signal as background adjusted signal. [22] The data were then normalized by the quantile method using the limma function normalize Between Arrays. [23] Probeset filtering was performed using the Agi4x44PreProcess filter method using AFE provided flags to identify features with quantification errors of the signal, reducing the 45,015 features on the Agilent Whole Human Genome Microarray 4×44K array by 34%, or 15,339 features. Residual technical batch effects were corrected using the ComBat method implemented in the SVA R package. [24] Since patients were not randomised to the ABC treatment all samples were treated as unrelated. Gene expression differences between treatment arms were assessed using the linear modeling features implemented in the limma package in R [25]. Multiple testing across contrasts was performed using the limma global method with Benjamini-Hochberg approach for control of false discovery rate (FDR). Gene set enrichment analysis of GO terms and KEGG pathways was performed using the conditional hypergeometric test implemented in the GOstats R package [26]. Statistical significance of gene set enrichment was performed using the Benjamini-Hochberg correction method implemented in the EMA R package [27]. Microarray data have been submitted to the Gene Expression Omnibus (GEO) and are available under NCBI accession number GSE 62117

Ethics Statement
The study was approved by the South Birmingham Research Ethics Committee (Reference 05/Q2707/40) and the Medicines and Healthcare Products Regulatory Agency-MHRA (Eudract number 2005/004021-26). All patients signed an informed consent to participate in the study.

Global Gene Expression
Multi-dimensional scaling (MDS) on all 103 samples was performed using profile level gene expression data from individual samples to see if subjects group according to treatment, time points or both. Distances were plotted using the limma plotMDS function with default parameters. Along the first dimension there was a clear separation of samples by samples time whereas in the second dimension the greatest distance was between the ABC treatment groups (6 m and 18-24 m) versus the other treatment groups including the untreated patients and the controls (Fig. 2).
Groups of differentially expressed genes were next generated using the R/Bioconductor limma library. Statistical criteria for identification of differentially expressed genes were FDR adjusted PÃ < 0.01 and fold change greater than 2.0. The cross sectional comparisons Table 1 showed the difference between groups was particularly striking at 6 months with ABC showing 4,185 differentially expressed genes (DEGs) compared to controls while the number of DETs for AZT and TDF were 895 and 381 respectively. Gene expression between the three groups at 18-24 months appear to have converged. At 18-24 m there were only 5 DETs between AZT and TDF, 88 between AZT and ABC and 251 between TDF and ABC (Table 1).

Gene enrichment analysis a. Comparison with controls
To gain mechanistic insight statistical list DEGs for each treatment contrast were tested for association to Gene Ontology (GO) terms and KEGG terms using the conditional Figure 2. The multi-dimensional scaling plot (MDS) shows that by and large ABC 6 and ABC 18-24 samples are separated from other samples in the study. Distances in the plot represent root-mean-square deviation (Euclidean distance) for the top 500 genes that best distinguish those samples. hypergeometric test algorithm implemented in R package GOStats (Falcon and Gentleman, 2007). Statistical significance of gene set enrichment was performed using the Benjamini-Hochberg correction method. Selected KEGG pathways showing over-representation in the three treatment groups compared to the controls are shown in the ( Table 2). Adherence junction was significantly enriched compared to controls in the ABC group at both time points but not in the other two treatment groups. Six months after treatment 24 genes were up-regulated and 8 genes were down regulated compared to controls in the ABC group (<0.01). At 18-24 months after treatment the corresponding figures were 22 upregulated and 7 down regulated (p<0.01) Individual adherence junction genes significantly up-regulated or down-regulated in the ABC treatment arm compared to controls are tabulated in Table 3.
Focal adhesions and tight junctions were also significantly enriched compared to the controls in the ABC treated group at 6 months but not in the other treatment groups. (Table 2) GO terms significantly enriched with ABC but not the other treatment groups at 6 months after treatment include cell adhesion ( b. Comparison between treatment arms at each time point We next performed enrichment analysis using statistical list of DEG from contrast of the three different treatment regimens at 6 and 18-24 months (excluding the HIV negative samples) in a cross sectional analysis. A Venn analysis was performed for each time point and Volcano plots were created filtering by fold change >1.5 and Benjamini-Hochberg corrected p-value <0.05 at 6 months and 18-24 months (Fig. 3). Six months after treatment enrichment of gene pathways involved with WNT signalling and leukocyte transendothelial migration were seen in ABC compared to patients receiving TDF and neuroactive ligand receptor interaction, leukocyte transendothelial migration and cytokine-cytokine receptor interactions compared to patients receiving AZT (Table 4).
KEGG pathways significantly over expressed when ABC was compared with the combined other two treatment groups at six months (intersection in Venn analysis) included ECM Table 3. Expression of genes controlling adherence junction in ABC treated group in comparison with controls at 6 months and 18-24 months after treatment.
At 18-24 months there were insufficient genes for a KEGG pathways or GO term gene enrichment analysis.
c. Comparison with HIV treatment naïve Enrichment analysis was next performed using statistical list of DEG in the three treatment groups and HIV infected treatment naïve subjects. The AZT and TDF pre-treatment samples were combined (n = 31) to form the HIV treatment naive group. These were then contrasted with the three treatment groups at 6 m and 18-26 months using limma R/Bioconductor limma library. A Venn analysis was performed for each time point and the unique gene sets expressed were tested for over representation analysis of KEGG pathways and GO-terms using the standard Hypergeometric test.
Significant number of genes were only available for a comparions between treatment groups six months after treatment, but not at the later time point. There was enrichment of the MAPK signalling pathway for ABC group compared to HIV treatment naïve patients (p<0.05) which was not seen with the other two groups. None of the pathways relating to environmental processing or signal transduction was enriched in the TDF and AZT treatment groups.

d. Individual gene analysis
In view of the findings showing enrichment of pathways relating to endothelial function we examined a number of genes relating to endothelial adhesion molecules some of which that had previously been reported to be affected by HIV and its treatment as well as some genes involved with tight junction, adherence junction, leukocyte transendothelial migration, and cytokine-cytokine receptor interaction. As can be seen (Table 5) there was over expression of some genes with ABC treatment in comparison with AZT and TDF. Differences between the three groups disappeared in all but one of the genes examined at 18-24 months. Moreover, when we performed enrichment analysis against the GO term 'biological adhesions' (GO:22610) compared to ABC six months after treatment, of the 39 genes expressed, 21 genes in the TDF group and 12 genes in the AZT group were either significantly up-or down regulated (p<0.05). There was no significant difference when AZT and TDF were compared at 6 months. However at 18 m the three regimens had converged such that significant enrichment was only seen in two genes for TDF and one gene for AZT compared to ABC, with the two latter treatment groups again showing no difference when compared against one another. A similar difference was observed when the three treatment groups were compared with controls with 29 genes in GO term 'biological adhesions' registered as significant at 6 months for ABC against none for the other two groups while the number of genes being significantly   Table 6. Compared to controls 48 and 17 genes controlling GO terms 'cell junction assembly' and 'apical junction' respectively were up-regulated and 17 and 4 genes were down-regulated (p<0.01). Similarly, for GO term 'cell adhesion', 172 genes were up-regulated and 17 were down-regulated in ABC-treated patients for 6 months compared to controls. A complete list of Kegg pathways that are enriched for all three treatment groups at different times after initiating treatment compared to controls is shown on Table 7.

e. Validation of microarray
In order to validate the microarray study we looked at genes that we had previously shown to be significantly up or down-regulated in the adipose tissue by real time PCR on the same subcutaneous fat samples. [19,20] Thus we found that compared to the controls hexose 6phosphate dehydrogenase (H6PD) was up-regulated (fold change 4.1, p<0.0001) in the AZT treated group as was also found with PCR. Expression of 11-beta hydroxysteroid dehydrogenase-1 (11βHSD-1) was down-regulated compared pre-treatment expression for both AZT and TDF (p = 0.001 and p<0.01 respectively) in keeping with real time PCR. Similarly, uncoupling protein 2 (UCP-2) was up-regulated for AZT treatment group at 18-24 months (p = 0.009) compared to pre-treatment values and UCP-3 down-regulated for both AZT and TDF (p = 0.007 and p = 0.004 respectively), as previously reported for PCR. Cytochrome BB (CYBB) and cytochrome B5 were down-regulated (p = 0.03 and p = 0.009 respectively) after 18-24m AZT treatment. Moreover CYBB was significantly up-regulated in the TDF 18-24m group compared to AZT treatment after 18-24 m (p = 0.001). These findings confirm our previous report that treatment with AZT, compared to TDF appears to impair mitochondrial oxygen transport genes. [19,20] Biological validation also comes from the observation that KEGG pathway enrichment was observed for fatty acid metabolism (TDF, p = 0.04), glycolysis and gluconeogenesis (TDF, p = 0.04), and steroid hormone biosynthesis (AZT, p = 0.03) at six months. At 18-24 months after treatment PPAR signalling pathway (OR 5.4, p = 0.023) and biosynthesis of fatty acids (OR 9.2, p = 0.026) were enriched in all three treatment groups, an observation in keeping with return to health reported in many studies, and the observed increase in limb and central fat in the first 24 weeks of antiviral therapy. [28] The consistency of the direction of changes in individual genes in the cross-sectional comparisons across the treatment groups at both time intervals (Table 6) strengthens further the validity of the study.

Discussion
In this study we generated microarray data from mRNA prepared from subcutaneous fat samples from HIV negative controls, treatment naïve HIV infected individuals, and after 6 months and 18-24 months into their initial antiretroviral treatment. The regimens contained AZT, TDF or ABC in combination tablets with lamivudine or emtricitabine. All three regimens also contained efavirenz. The patients were matched for age, ethnicity, sex and pre-treatment CD4 count and viral loads. None had experienced any AIDS defining conditions. We have shown that at all time points patients receiving ABC displayed a markedly different pattern of gene expression within subcutaneous adipose tissue, with a large percent of the overall variance accounted for by ABC compared to the other two treatment groups. The differences were particularly marked at 6 months after treatment initiation but were still evident at 18-24 months. After 6 months of treatment, 4185 genes showed significant differentially expressed sequences (DET) in patients receiving ABC relative to the controls while the number of DETs in the AZT and the TDF group were 895 and 381 respectively. After 18-24 months of treatment the differences between the three groups had narrowed. (Table 1) We next examined the specific pathways over expressed in the three treatment groups compared to controls and to treatment naïve HIV-infected patients. All groups showed an enrichment of PPAR signalling and biosynthesis of fatty acids in agreement with previous studies of the effect of antiretroviral therapy [19,20,29,30] and a return to health. [8,28] The ABC group, however showed enrichment of a number of pathways at both time points that was not observed with the other two groups. The pathways of interest are involved with cell processes and cell communication (adherence junction, focal adhesions, tight junctions) and environmental information processing (WNT signalling, ECM receptor interaction, leukocyte transendothelial migration, neuroactive ligand interactions, MAPK signalling). Some of these may have important roles in cardiovascular disease. [31][32][33][34][35] With the greatly increase survival observed with the use of combination anti-retroviral therapy long term complications of such therapies have assumed increasing importance. HIV itself [36] and antiviral treatments are associated with an increased incidence of cardiovascular disease [10] which is not entirely explained by traditional risk factors including adjustment for lipids [11]. In large cohort studies the increased risk of cardiovascular disease was associated with the use of protease inhibitors [10] and current or recent use of ABC and less consistently, didanosine. [10][11][12][13][14]18] Pooled analysis of two cohort studies [11,14] demonstrated a relative risk of 1.91 (95% CI: 1.5, 2.4) for use of ABC within previous six months and risk of myocardial infarction [37]. Other cohort studies [18,38] or randomised control trials (RCT) [16,39] have, on the other hand, not confirmed the association between ABC and cardiovascular events. Patients recruited in RCT, however, tend to be younger and followed for a shorter period of time than cohort studies. Hence despite their limitations results from large cohort studies cannot be entirely ignored [37]. Current International AIDS Society -USA Panel guidelines [40] advise avoidance of ABC use in patients at high risk of cardiovascular disease.
Possible mechanisms for the cardiovascular effects of ABC include an association with inflammatory markers, [41][42][43] abnormalities in endothelial cell function, [32,44] and platelet dysfunction. [45][46][47][48] These results have not been confirmed by other studies. [49][50][51][52] Thus the potential role of ABC in cardiovascular disease and the possible mechanisms through which this can be effected remain to be established.
The observation that the increased risk of cardiovascular events was only seen with recent or current use of ABC, [11,14,18,41] particularly in those in the higher Framingham risk group [11] suggests an additive effect of ABC to an already pre-existing underlying pathological mechanisms such as plaque instability, endothelial function, or thrombosis. [53] Our study provides another possible mechanism that requires further research. Tight junctions (TJ) focal adhesions and adherence junctions (AJ) were overrepresented in the KEGG pathway and 'biological adhesions' in GO terms 6 months after ABC treatment. In particular TJ and AJ are interlinked and important for maintaining tissue architecture [33,54,55] and have been implicated in cardiovascular disease. The gene VE-cadherin at AJ up-regulates TJ gene claudin-5 by limiting the nuclear accumulation of beta-catenin, which represses the claudin-5 promoter. [54] We showed down regulation of claudin-5 (CLDN5) and up-regulation of catenins (CTNNA1, CTNNB1) compared to controls and to AZT and TDF at 6 months (Tables 3 and  5) while there was no observable change in the expression of VE-cadherin (data not shown). These changes may result in increased endothelial permeability in subcutaneous adipose tissue in patients treated with ABC. [54] Were such time related changes in vascular permeability to be replicated in other tissues, and in particular coronary vasculature, an additional possible mechanism for the observed increase in the incidence of cardiovascular disease in patients starting abacavir treatment with underlying cardiovascular risk would be suggested. Further studies could provide insight into this suggested mechanism.
Pathways involved with environmental information processing (ECM-receptor interactions, MAPK signalling, neuroactive ligand interaction, Wnt signalling pathway) tended to be either enriched earlier in ABC ( Table 2) or were significantly different at 6 months (Tables 4 and 5) and GO term 'biological adhesions' 'cell adhesion' and 'cell matrix adhesion', 'cell-cell signaling' and 'signaling' were overrepresented in the ABC group compared to the other two groups.
Cell adhesion molecules (CAMs) [31,33] and inflammatory markers are also believed to play an important role in atherosclerosis. However, reports of abnormalities in circulatory CAMs and inflammatory markers in endothelial function after ABC treatment are conflicting. [42,42,44,51,52,[56][57][58] as are the nature of the relationship between circulating inflammatory markers and the expression of adhesion molecules at the endothelial surface. [34,59,60] De Pablo et al. reported that ABC significantly induces leucocyte accumulation through activation of MAC-1 which reacts with its endothelial ligand ICAM-1 both in vitro and in vivo [56,61]. The demonstration of enrichment of genes involved with leukocyte transendothelial migration (Table 4) and up-regulation of ICAM 1, ICAM-4 and ICAM-5 compared to controls and, in the case of the latter two compared to the other two groups, after six-month ABC treatment (Table 5) are compatible with these findings.
Our study does have several limitations. Patients receiving AZT and TDF were randomised while those on ABC were only matched as to gender and ethnicity. While there was no statistical difference in weight (pre-treatment and at test), BMI, CD4 count, smoking history, and family history of diabetes and cardiovascular disease between the three groups, an unintended bias cannot be excluded. Moreover the number of subjects in each treatment group was necessarily small, and not all subcutaneous fat samples were available for study, further limiting the conclusions. However we adjusted for multiple comparisons by using the Benjamini-Hochberg False Discovery Rate (FDR) as well as Bonferroni multiple testing corrections. Finally as discussed above any extrapolation for gene expression from samples obtained from one tissue must only be guardedly extended to the whole subject, not least because we were unable to comment on post-transcriptional and regulation or post-translational modifications of gene expression.
Our study has highlighted novel mechanisms by which ABC may interfere with endothelial adhesion process and endothelial function in subcutaneous adipose tissue, particularly early in the treatment. In particular abnormalities in tight junction and adherence junction brought about by ABC use, if observed in other tissues, may be a contributing factor to cardiovascular events in persons already with some underlying risk for atherosclerosis [62], and may require further study.