Association of HADHA expression with the risk of breast cancer: targeted subset analysis and meta-analysis of microarray data

The primary dataset

The microarray dataset used in the present study is available for public use on the Gene Expression Omnibus website http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS3716 of the National Institutes of Health, USA. Details of the study subjects on whom these microarray studies were conducted have been described previously [5]. Briefly, the dataset comprises microarray data collected through Affymetrix Human Genome U133A platform that measures expression of 22,283 genes. The data were collected using histologically normal epithelium from four sets of subjects--those who underwent reduction mammoplasty (n = 18), those who underwent preventive mastectomy (n = 6), estrogen receptor positive (ER+) breast cancer patients (n = 9) and estrogen receptor negative (ER-) breast cancer patients (n = 9). The data were available in normalized format.

Targeted subset analysis

Replication of the results: meta-analyses

We also aimed to ensure that the results obtained from one microarray dataset were robust and could be replicated in other datasets. We queried the Oncomine database and retrieved microarray data from other relevant studies. We studied the association of gene expression with three outcomes--risk of breast cancer, risk of metastasis and risk of recurrence. We then combined these datasets meta-analytically using the random effects model of DerSimonian and Laird [13, 14]. For these analyses the effect size was measured and expressed as the standardized mean difference (SMD) and its 95% confidence intervals. The Oncomine website reports the results as means, medians, quartiles and minimum and maximum values. Since the random-effects model assumes normal distribution of the effect measures, we first estimated the mean and standard error for each group (for example, for subjects with breast cancer; subjects with a metastatic event or subjects with recurring breast cancer) using the method described by Hozo et al. [15] We then estimated the SMDs. To depict the potential variability in the HADHA expression based on the probes used by individual studies, we conducted the meta-analyses separately for each combination of the study and the probe used. Each comparison represented a specific combination of the included study and the reporter used in the study. The between-study heterogeneity in this meta-analysis was examined using the I² statistic. Since expression data on all individual subjects was available for the outcome of risk of breast cancer, we also conducted individual patient data (IPD) meta-analysis [16]. For this we used the clustered unconditional logistic regression analyses [16, 17] with disease status as a dichotomous dependent variable, comparison-specific z-scores as the predictor variable and comparison indicator as the clustering variable. Comparison-specific z-scores were estimated as the relative deviates (mean expression/standard deviation of expression) within each comparison group.

Other statistical analysis

To quantify and test differential gene expression, we used two-tailed Student's t tests for unpaired samples. The clinical and statistical significance of the findings were presented as volcano plots. To account for multiple testing, we estimated the false discovery rates (q values) using the QVALITY software program [18]. Discriminant utility of each gene was assessed using non-parametric receiver operating characteristic (ROC) curve analysis. To group subjects based on their HADHA expression, we used a k-means clustering approach. All statistical analyses were conducted using Stata 10.0 software package (Stata Corp, College Station, Texas). We aimed for a type I error rate of 0.05 and a false discovery rate of 0.15.

Differential expression analyses

Using the shortlisted set of 47 genes shown in Table 1, we first determined if these genes were differentially expressed in subjects with cancer (n = 18) and those without (n = 24). The volcano plot (Figure 1A) showed that seven of the 47 genes were significantly differentially expressed between these study groups. These genes included five over-expressed genes (AQP7, PLTP, PCK2, GCDH and ARSA) and two under-expressed genes (ACSL5 and HADHA). Of these, HADHA was the most significant statistically. To account for the possible covariance among these gene expression values we conducted stepwise multivariate analyses using unconditional logistic regression and observed that only two genes--HADHA and ARSA--were retained in the final model (Figure 1B). This model explained 35% of inter-individual variability in breast cancer susceptibility with a predictive accuracy of 86.8%. Interestingly, when HADHA expression was removed from this model the ARSA lost its statistical significance but removal of ARSA did not affect the statistical significance of HADHA. This indicates that HADHA gene expression was the most important statistical predictor of altered risk of breast cancer.

Does ER status influence the expression of HADHA?

To examine if this association could be influenced by the ER status, we conducted three sets of analyses. First, we studied whether HADHA expression was different based on the ER status. We found that the mean HADHA was not significantly differentially expressed by ER status (mean HADHA expression in subjects with ER + breast cancer = 6.00; in subjects with ER-breast cancer = 5.90; p = 0.462). Second, we adjusted the standard error estimates for the ER status using clustered logistic regression and observed that the statistical significance for the HADHA gene expression further increased (p = 0.0001) while that of the ARSA gene decreased (p = 0.082) indicating that the association of HADHA was unlikely to have been influenced by the ER status. Third, we constructed volcano plots and conducted stepwise logistic regression analyses by comparing the ER + and ER-subjects separately with subjects without cancer as the reference group. We observed (Figure 1C-F) that HADHA gene expression was the only consistent predictor across ER status but more so in the ER-subjects. Indeed, the q value for the HADHA gene was 0.15 for the cancer versus no cancer comparison, 0.13 for the ER-versus no cancer comparison but 0.88 for the ER + versus no cancer comparison. Two other genes (UCP3 and DGAT1) were retained in the final model of stepwise regression analyses when ER-subjects were compared to the no cancer group however this association was not observed when ER + subjects were compared to the same reference group.

Graded risk of breast cancer based on HADHA expression

We next considered whether the association of HADHA gene expression with risk of breast cancer exhibited a threshold effect or whether it was a graded dose-response. For this, we used two approaches. First, we normalized the gene expression in the no cancer group to 100%. We found (Figure 1G) that the HADHA expression had fallen to 73% (95% CI 64%-83%) in subjects with cancer; with a higher expression in ER + subjects (76% of the no cancer group, 95% CI 61%-91%) than in ER-subjects (70% of the no cancer group, 95% CI 55%-85%). Second, the k-means clusters (which explained 95.9% of the variability in HADHA expression) clearly demonstrated a dose-response association (Figure 1H) such that more severe down-regulation of HADHA was associated with a greater risk of being in the breast cancer group.

Meta-analyses of the differential expression of HADHA

We first observed that the mean expression levels for HADHA probes (expressed as log transformed values) were widely different across the six studies (Zhao et al. [24]:-0.33, Radvanyi et al. [21]: 3.08, Richardson et al. [22]: 5.29, Karnoub et al. [20]: 2.97, Turashvili et al. [23]: 3.92 and Finak et al. [19]:-2.65). We therefore transformed these values into comparison-specific z-scores (mean expression for a comparison/standard deviation of expression for that comparison). Upon this z-transformation, all the studies had a mean z-score of 0 and a standard deviation of 1. We conducted meta-analyses on the HADHA expression z-scores. Using the DerSimonian and Laird model, we observed (Figure 2) that the summary SMD (filled diamond in Figure 2) was-0.48 (95% CI-0.84--0.11). Considering the statistical properties of SMD it is possible to transform this into probability [25]. This transformation indicated that there was an average 18.4% reduction in expression of HADHA (95% CI 4.5%-30.0%) in cases of breast cancer as compared to normal subjects. Interestingly, this significant reduction in the expression of HADHA was observed in spite of the high degree of heterogeneity (I² 64.6, p < 0.001, pie-chart in Figure 2) between the comparisons due to different cancer subtypes, reporters used in various studies and other study characteristics.

We observed that the invasive ductal carcinoma (p = 0.046) and unspecified invasive breast carcinoma (p = 0.005) showed a significant under-expression of HADHA gene but lobular carcinoma (p = 0.781), invasive lobular carcinoma (p = 0.780) or invasive mixed carcinoma (p = 0.717) did not show a significant alteration of HADHA gene expression. Alternatively, we conducted the IPD meta-analysis using logistic regression analyses. We found that the odds ratio for breast cancer was 0.74 (95% CI 0.60-0.92) after clustered analyses. Thus, there was a 26% reduction in the risk of breast cancer per unit increase in z-scores. These values show a striking resemblance with the findings observed in the Graham et al. dataset and demonstrate the replicability of our findings.

We also investigated if HADHA expression was associated with an altered risk of metastasis and recurrence. For risk of a metastatic event we found nine studies [26–33] representing 430 metastatic events and 1104 metastasis-free cancers (Figure 3). Subjects who developed a metastatic event during follow-up had a reduced HADHA expression (summary effect size-0.65, 95% CI-1.47-0.16) but this was not statistically significant (p = 0.117). Also, there was a very high degree of between-study heterogeneity (I² 98.8%). Similarly, for the outcome of the risk of recurrence (Figure 4), we found that there were 10 studies [19, 27, 29, 30, 32, 34–38] representing 484 recurring and 890 non-recurring breast cancers. Meta-analysis demonstrated that although there was a consistent decrease in average HADHA expression in patients with a recurring form of breast cancer (summary effect size-0.60, 95% CI-1.44-0.24), the finding was neither statistically significant (p = 0.160) nor homogeneous (I² = 98.7%) across studies.

Biological plausibility

Biological significance of our novel observations should be considered in the light of the following facts. First, the HADHA gene (chromosomal location 2p23) codes for the four alpha chains in the 8-meric mitochondrial trifunctional protein (TFP) [40]. This enzyme performs three cardinal functions in the β-oxidation of long chain fatty acids by catalyzing the activities of the 2-enoyl-CoA hydratase (ECH), L-3-hydroxyacyl-CoA dehydrogenase (HACD) and 3-ketoacyl-CoA thiolase (KACT). Of these three, the first two enzymes (ECH and HACD) are specifically catalyzed by the alpha chains of TFP. Severe deficiency (< 50% of normal activity) of TFP is known to be associated with life-threatening manifestation of the long chain 3-hydroxyacyl-CoA dehydrogenase deficiency [41]. However, the effects of a milder deficiency of TFP (for example, activity between 50%-80% of the normal) are currently unknown. Our results indicate that breast cancer patients had 18-30% decreased expression of HADHA gene. We therefore hypothesize that there may be a compromised metabolism of long chain fatty acids in breast cancer due to a relative deficiency of the alpha chains of TFP. In this context, it is noteworthy that a recent large genome-wide association study [42] found a strong association of breast cancer with a polymorphism in the gene encoding enoyl CoA hydratase domain containing 1 (ECHDC1), which also partakes in the integrity of the TFP.

Second, the efficacy of β-oxidation of n-3 and n-6 long chain fatty acids can be tissue- and location- specific. For example, in rat livers it has been shown that the n-3/n-6 ratio influences peroxisomal but not mitochondrial β-oxidation [43]. In contrast, mitochondrial β-oxidation of long chain fatty acid has been implicated in breast cancer pathogenesis [42]. We also could not demonstrate a significant association of the genes involved in the PPAR-γ pathway reinforcing the possibility that mitochondrial rather than peroxisomal β-oxidation of long chain fatty acids may be more critical in breast carcinogenesis. Third, HADHA occupies an important position in the network of genes that have been implicated in autophagy and apoptosis [44]. Finally, triangulation of the following facts lends additional credence to our observations: i) intact epithelium of mammary glands has the ability to act as stem cells for carcinogenesis [2]; ii) n-3 long chain fatty acids have the ability to target such stem cells [45]; and iii) HADHA is involved in the mitochondrial β-oxidation of long chain fatty acids. Together these observations from published literature strongly support the biological plausibility of our finding that HADHA is differentially expressed in subjects with and without breast cancer.

Limitations

Our study has all the limitations implicit in any microarray association study and meta-analyses. In addition, there are three more limitations. First, although there is a strong circumstantial evidence that favors an inference of HADHA expression-breast cancer association, it must be realized that robust functional studies are required before this association can be conclusively claimed. Our study does not have a component of functional assays that can help put these results in a biological perspective. Second, due to limitations imposed by the microarray platform used in the primary study, we could not evaluate the potential association of a large number of additional lipid and fat metabolism related genes with the risk of breast cancer. Inclusion of those genes may not only affect the q values associated with HADHA but also may provide a more comprehensive understanding of the role of fatty acids in breast cancer. Thirdly, although consistent, the observed differential expression of HADHA with cancer progression (as reflected by risk of metastasis and recurrence) is statistically inconclusive.