Transcriptome Analysis of Genes Regulated by Cholesterol Loading in Two Strains of Mouse Macrophages Associates Lysosome Pathway and ER Stress Response with Atherosclerosis Susceptibility

Cholesterol loaded macrophages in the arterial intima are the earliest histological evidence of atherosclerosis. Studies of mouse models of atherosclerosis have shown that the strain background can have a significant effect on lesion development. We have previously shown that DBA/2 ApoE−/− mice have aortic root lesions 10-fold larger than AKR ApoE−/−mice. The current study analyzes the response to cholesterol loading of macrophages from these two strains. Macrophages from the atherosclerosis susceptible DBA/2 strain had significantly higher levels of total and esterified cholesterol compared to atherosclerosis resistant AKR macrophages, while free cholesterol levels were higher in AKR cells. Gene expression profiles were obtained and data were analyzed for strain, cholesterol loading, and strain-cholesterol loading interaction effects by a fitted linear model. Pathway and transcriptional motif enrichment were identified by gene set enrichment analysis. In addition to observed strain differences in basal gene expression, we identified many transcripts whose expression was significantly altered in response to cholesterol loading, including P2ry13 and P2ry14, Trib3, Hyal1, Vegfa, Ccr5, Ly6a, and Ifit3. Eight pathways were significantly enriched in transcripts regulated by cholesterol loading, among which the lysosome and cytokine-cytokine receptor interaction pathways had the highest number of significantly regulated transcripts. Of the differentially regulated transcripts with a strain-cholesterol loading interaction effect, we identified three genes known to participate in the endoplasmic reticulum (ER) stress response, Ddit3, Trib3 and Atf4. These three transcripts were highly up-regulated by cholesterol in AKR and either down-regulated or unchanged in loaded DBA/2 macrophages, thus associating a robust ER stress response with atherosclerosis resistance. We identified significant transcripts with strain, loading, or strain-loading interaction effect that reside within previously described quantitative trait loci as atherosclerosis modifier candidate genes. In conclusion, we characterized several strain and cholesterol induced differences that may lead to new insights into cellular cholesterol metabolism and atherosclerosis.


Introduction
Atherosclerosis is a complex and progressive pathology of arteries that can be initiated by the accumulation and entrapment of lipoproteins in the extracellular matrix of the sub-endothelial intima layer [1,2].Monocytes migrate across the endothelial monolayer into the intimal layer where they differentiate into macrophages and take up modified low density lipoproteins by scavenger receptors [3,4].Once inside the cells modified LDL is degraded in lysosomes, depositing free cholesterol in the lysosomes that is transported to the endoplasmic reticulum for conversion into cholesterol esters and storage in lipid droplets [5,6].Cholesterol esters inside the lipid droplets can be hydrolyzed back into free cholesterol; and, recently, this hydrolysis has been shown to be mediated via autophagy delivering cholesterol esters to lysosomal acid lipase [7].Free cholesterol can be exported out of the cells through passive diffusion, and/or protein mediated efflux via ATP-binding cassette transporters ABCA1, ABCG1 or scavenger receptor SR-BI [8][9][10][11][12].Macrophages cannot limit their uptake of cholesterol as the expression of scavenger receptors is not subject to feedback regulation by cellular cholesterol content [13,14].An imbalance between cholesterol influx and efflux leads to excessive accumulation of cholesterol in macrophages transforming them into cholesterol-engorged foam cells characteristic of fatty streaks [15][16][17].
The genetic background strain has been shown to effect atherosclerotic lesion area in hyperlipidemic ApoE2/2 mice [18].We have previously shown that chow-diet fed DBA/2 ApoE 2/2 mice have aortic root lesion areas that are 10-fold larger than those in AKR ApoE 2/2 mice at 16 weeks of age [19].To gain insight on the effects of cholesterol loading on macrophage gene expression and its relation to atherosclerosis susceptibility, we incubated bone-marrow derived macrophages from atherosclerosis resistant AKR ApoE 2/2 and atherosclerosis susceptible DBA/2 ApoE 2/2 mice with acetylated LDL (AcLDL).We found that DBA/2 cells accumulate more total and esterified cholesterol, but less free cholesterol than the AKR macrophages.We investigated the effect of cholesterol on the macrophages transcriptome and identified transcripts that were differentially regulated by strain or cholesterol loading, and those transcripts in which there was a strain-loading interaction effect.Gene set enrichment analysis showed that some of the identified transcripts are involved in the lysosome pathway, and we found several transcripts that play a role in the endoplasmic reticulum stress response.We hypothesized that some of these differentially regulated transcripts could be responsible for the observed difference in atherosclerosis susceptibility between the AKR and DBA/2 strains.We identified differentially expressed transcripts as possible atherosclerosis modifier candidate genes residing within atherosclerosis QTLs previously characterized in a strain intercross between these two strains [19].

Mice
Age matched AKR and DBA/2 ApoE 2/2 mice were maintained on chow diet and sacrificed at 16 weeks of age.Plasma total cholesterol levels at time of sacrifice were 372.8647.78mg/ml and 912.76104.3mg/ml for AKR ApoE2/2 and DBA/2 ApoE2/2 mice, respectively.Femurs were collected to isolate and culture bone marrow monocytes.
Ethics Statement.All experiments were performed in accordance with protocols approved by the Cleveland Clinic Institutional Animal Care and Use Committee.

Total, free, and esterified cholesterol quantification
Bone marrow derived macrophages (BMM) from AKR and DBA/2 ApoE 2/2 mice were cultured in p100 cell culture dishes for 13 days with 15% L-cell conditioned media, a source of macrophage colony stimulating factor (M-CSF).We performed control studies to compare BMM cholesterol loading by 24 h incubation with 50 mg/ml AcLDL, copper oxidized LDL, or native LDL; and, we only observed robust cholesterol loading with AcLDL at this dose, justifying its use in our subsequent studies.The cells were loaded for 48 h with 50 mg/ml AcLDL to allow foam cell formation with unloaded cells used as controls.Lipids were extracted and total, free and esterified cholesterol levels were quantified and normalized to cell protein as previously described [20].Quadruplicate samples were assayed in triplicate and significance was determined by ANOVA with Newman-Keuls posttest using GraphPad Prism software.

Loading of macrophages with acetylated LDL for transcriptome profiling
Differentiated BMM in p100 dishes from AKR and DBA/2 ApoE 2/2 mice were incubated for 24 hours in the presence or absence of 50 (experiment 2) or 100 mg/ml (experiment 1) AcLDL in quadruplicate.Different AcLDL doses were used due to altered efficiency of cholesterol loading of independent AcLDL preparations.After the incubation, the media and cholesterol were aspirated off, cells were washed twice, scrapped in 1 ml ice-cold PBS and transferred to microcentrifuge tubes.Cell pellets were obtained by centrifugation in a microcentrifuge at 5,000 rpm for 2 minutes at RT and kept on ice for subsequent use.

Isolation of total RNA from BMM cell pellets
Total RNA was isolated from cell pellets using the RNeasy Mini Kit (Qiagen, Valencia, CA), following the manufacturer's protocol.On-column digestion with RNase-free DNase (Qiagen, Valencia, CA) was performed to remove genomic DNA.DNase was removed in the subsequent washing steps.RNA integrity was tested by overnight incubation of 200-500 ng of total RNA at 37uC and observation of the 18S and 28S ribosomal RNA bands on a 1.2% agarose/ethidium bromide gel.

Hybridization and detection of gene transcripts
An aliquot of total RNA for each sample (, 2 mg) was submitted to the Genomic Core at Lerner Research Institute.Total RNA was used as template to synthesize single stranded cDNA following the Illumina protocol.Single stranded cDNA is then converted into double stranded cDNA, purified and in vitro transcribed (in the presence of biotinylated UTP and CTP) to produce biotin-labeled cRNA.Purified labeled cRNA from experiment 1 and 2 samples were hybridized respectively to MouseRef-8_v1 and Ref-8_v2 microarray chips (Illumina) for 16 hours at 58uC, according to the manufacturer's instructions.Eight samples were profiled on one chip with either 24,613 or 25,697 transcript probes per sample.To reduce the chip-to-chip variability two control and two cholesterol-loaded samples from each strain were put on each microarray chip.After hybridization, the arrays were washed several times and stained with streptavidin-Cy3.After the staining, the arrays were washed and scanned on the Illumina BeadArray Reader using Illumina BeadScan image data acquisition software.

Microarray data analysis
Expression data were quantile normalized and log2 transformed using the R-package lumi [21].1,342 unique probes that hybridized to transcript sequences containing a strain-specific polymorphism or an indel (insertion/deletion) were excluded from further gene expression-related analysis.Sequence variations between AKR and DBA/2 (SNPs and insertions-deletions) were obtained from sequencing data generated by the Mouse Genomes Project: http://www.sanger.ac.uk/resources/mouse/genomes/.Probes were also removed if no sample had a detection p-value less than 0.01.9,316 and 11,919 probes remained for further analysis in experiments 1 and 2, respectively.To determine strain effects on the transcriptome, only the results of unloaded bone marrow macrophages were analyzed.Log2 fold change, unadjusted p-values, and false discovery rate adjusted p-values were determined separately for the two independent datasets with the limma package in R [22,23] using a fitted linear model with the following equation: Y i = b 1 Strain+b 2 Loading+b 3 Strain-Loading, with the strain and loading bs as additive covariates and the strainloading interaction b as an interactive covariate.
Gene set enrichment analysis was performed for all expressed transcripts to identify possible pathways altered by strain, loading and strain-loading interactions.The romer function in the limma package was used to test for gene set enrichment analysis [24,25] with the Molecular Signature Database from the Broad Institute as gene sets [26].We present and discuss further only the pathways that were significantly enriched using the permutation test p-value threshold of 1.0E-04.The MIAME compliant microarray dataset from this study is available through the Gene Expression Omnibus server, accession number GSE38736.

Real-Time quantitative PCR (qPCR)
To validate the findings of gene transcript measurements by microarray, we assessed the expression of selected gene transcripts by using qPCR approach.The expression levels of tribbles homolog 3 (Trib3), activating transcription factor 4 (Atf4), ceroidlipofuscinosis neuronal 5 (Cln5), and ATP-binding cassette subfamily G member 1 (Abcg1) were quantified relative to the endogenous control gene, ubiquitin C (Ubc) using pre-designed TaqMan gene expression assays (Applied Biosystems).Mean fold changes for each sample were calculated as previously described [27].

Western blot
Unloaded and cholesterol loaded BMM from AKR and DBA/2 were lysed with 0.5 mL of M-PER mammalian protein extraction reagent (Pierce Biotechnology).Cell lysates (45 mg protein/lane) were loaded, separated on a 4-15% gradient polyacrylamide gel and transferred to PVDF membranes by semi-dry electroblotting.After blocking with rapid blocking buffer (Amresco, Cat # M325) for 1 hour, the membrane was incubated overnight with rabbit antibody to TRIB3 (Proteintech, Cat # 13300) at 4uC.The membrane was washed and further incubated with HRP goat antirabbit IgG (Amresco, Cat # N791) secondary antibody for 1 hour.The membrane was then exposed to an enhanced chemiluminescent system, and bands were visualized by exposure to X-ray film.After stripping, GAPDH protein was visualized as loading control using goat antibody to GAPDH (Abcam, Cat # ab9483), followed by HRP rabbit anti-goat IgG as described above.TRIB3 band density was quantified using Total Lab TL120 software (Nonlinear Dynamics) and normalized to the respective GAPDH band density.

AKR and DBA/2 macrophages respond differently to cholesterol loading
Bone marrow derived macrophages of the two strains were incubated with 50 mg/ml AcLDL for 48 hours.We observed that macrophages from the atherosclerosis susceptible DBA/2 strain had significantly higher levels of total cholesterol (p-value,0.0001) and cholesterol esters (p-value ,0.0001) compared to atherosclerosis resistant AKR macrophages (Figure 1A and 1B).In contrast, AcLDL-loaded DBA/2 macrophages had significantly lower levels of free cholesterol (p-value,0.05, Figure 1C).We have recently found that the strain difference in cholesterol ester metabolism is due to reduced cholesterol ester turnover in DBA/2 cells, which we attribute to decreased autophagic flux; and, autophagy is the primary mechanism responsible for the degradation of cholesterol esters in lipid droplets [28].The current study focuses on the transcriptome differences between macrophages from the DBA/2 and AKR strains and the changes in response to cholesterol loading in these macrophages.Our objective was to identify possible genes and pathways that may play a role in the observed strain-specific differences in response to cholesterol loading and atherosclerosis.

Hierarchical clustering
Bone marrow derived macrophages from mice of both strains were incubated with or without AcLDL for 1 day (n = 4 per group, total of 8 groups) in two independent experiments.Total RNA was applied to Illumina expression arrays and analyzed as described in Materials and Methods.Hierarchical clustering analysis was performed on all 32 samples (Figure S1).The two independent experiments were assayed on different versions of the Illumina array, using cells from different sexes with different batches of AcLDL, and were treated and analyzed in different years.Based on the observed clustering profile of these two groups, we decided not to pool the samples, and thus we analyzed the expression data from the two experiments separately for strain, loading and strainloading interaction.Furthermore, we used the two studies to identify the overlapping group of genes and pathways that consistently respond to cholesterol loading in a strain-shared or strain-specific fashion.

Strain differences on BMM transcriptome
Experiment 1 samples.3,059 transcripts were identified as differentially expressed between AKR and DBA/2 unloaded macrophages at a stringent false discovery rate (FDR) adjusted pvalue,0.01 in experiment 1 samples.Table 1 shows the top 10 most significant differentially expressed transcripts, while the entire list is shown in Table S1.229 differentially expressed transcripts were found with a 2-or higher fold change in expression between the two strains, of which 137 transcripts were expressed higher in DBA/2 and 92 transcripts were expressed higher in AKR macrophages.
Experiment 2 samples.A similar analysis was performed on the experiment 2 samples and 1,703 transcripts were found to be differentially expressed between AKR and DBA/2 macrophages (FDR adjusted p-value,0.01;Table 1

top ten and Table S2 for all results
).Of the 418 differentially expressed transcripts with a 2-or higher fold change, 220 transcripts were expressed higher in DBA/2 and 198 were expressed higher in AKR macrophages.
Combined results and discussion.522 differentially expressed transcripts overlapped among the 3,059 and 1,703 strain All data are mean6SD, N = 4 per group using the average of triplicate assays per sample.P-values were calculated by ANOVA with Newman-Keuls posttest, showing only the strain differences after loading.There were no significant strain differences in unloaded cells.doi:10.1371/journal.pone.0065003.g001significant differences identified in the experiment 1 and 2, respectively (significance threshold set at 0.01 FDR, Table S3), making these transcript changes the ones that we are most confident of.Among the top ten differentially expressed transcripts, two and three were in the histocompatibility gene family in experiment 1 and 2, respectively; however, different members of genes in this family were observed in the two experiments.The only gene to make it to the top ten in both experiments was Prcp, encoding prolylcarboxypeptidase (also known as angiotensinase C), expressed significantly higher in DBA/2 macrophages.Transcripts encoding for transmembrane glycoprotein nmb (Gpnmb) and napsin A aspartic peptidase (Napsa) were expressed significantly higher in AKR macrophages.
Gene set enrichment analysis.To identify the common biological pathways most relevant to the genes that differ in expression between AKR and DBA/2 BMM, transcriptomes from both experiments were subjected to Gene Set Enrichment Analysis (GSEA) using KEGG pathways, and we report here only the pathways identified as significantly enriched in both.Strain effects on geneset enrichment were found for the hematopoietic cell lineage, chemokine signaling, toll like receptor signaling and aldosterone regulated sodium reabsorption pathways (permutation test p-value,0.0001,Table S4).
In conclusion, significant basal gene expression differences were observed between the two strains in both experiments, which need to be considered in the following analysis of gene expression changes in the AKR and DBA/2 macrophages in response to cholesterol loading.

Cholesterol loading effect on BMM transcriptome
To identify the differentially regulated transcript in response to cholesterol loading, the expression data were fitted in a linear model with strain as an additive variable and strain-loading interaction as an interactive variable.This model identified transcripts whose expression was either up-regulated or downregulated in one or both strains upon cholesterol loading.
Experiment 1 samples.3,758 transcripts were identified as differentially expressed in response to cholesterol loading in AKR and DBA/2 macrophages at an FDR adjusted p-value,0.01(Table 2

top ten and Table S5 for all results
).There were 261 differentially expressed transcripts with a 2-fold or higher change in expression upon cholesterol loading in one or both strains, of which 127 transcripts were up-regulated and 134 down-regulated.
Experiment 2 samples.A similar number of significantly differentially expressed transcripts were found in macrophages in response to cholesterol loading (3,308 transcripts; FDR adjusted pvalue,0.01;Table 2, Table S6).Of the 567 differentially expressed transcripts with a 2-fold or higher change, 236 transcripts were up-regulated upon cholesterol loading versus 331 that were down-regulated in one or both strains.
Combined results and discussion.The cholesterol loading effect on gene expression was largely reproducible in both experiments, despite microarray platform differences.There were 2,475 cholesterol regulated transcripts with identical probes on the two array platforms, and of these, 1140 transcripts were significantly regulated by cholesterol loading in both experiments (Table S7).The cholesterol-loading induced fold changes of these transcripts were also well correlated between the two experiments (R 2 = 0.68, Figure S2).Transcripts that were significantly up-regulated in response to cholesterol loading in both experiments include: tribbles homolog 3 (Trib3), hyaluronoglucosaminidase 1 (Hyal1), vascular endothelial growth factor A (Vegfa), etc. Transcripts that were significantly down-regulated in both experiments include: chemokine (C-C motif) receptor 5 (Ccr5), lymphocyte antigen 6 complex, locus A (Ly6a), interferon-induced protein with tetratricopeptide repeats 3 (Ifit3), etc.The related purinergic receptors P2ry13 and P2ry14 were also down regulated in both experiments, with P2ry13 identified as the most significantly altered transcript in response to cholesterol loading in experiment 1, while P2ry14 was identified among the top 10 most significant transcripts in experiment 2. P2ry14 encodes for a G-protein coupled receptor expressed in a subpopulation of bone-marrow hematopoietic stem cells [29].P2ry13 has been shown to play a role in hepatic uptake of holo-HDL particles, particularly in SR-BI deficient liver cells; and, P2ry13 knockout mice are reported to have decreased reverse cholesterol transport [30].However, whether P2ry13 deficiency has an effect on plasma HDL levels in mice is controversial [30,31].The most significantly altered transcript in response to cholesterol loading in experiment 2 was chemokine (C-C motif) receptor 5 (Ccr5), which was down-regulated by this treatment.Emerging evidence from both human and mouse studies supports important role(s) played by the Ccr5 receptor and its ligand Ccl5 in atherogenesis, a detailed description of which is provided in a recent review by Jones et al. [32].
Gene set enrichment analysis.Eight pathways were significantly enriched in transcripts regulated by cholesterol loading in both experiments: lysosome, cytokine-cytokine receptor interaction, primary bile acid biosynthesis, allograft rejection, aminoacyl tRNA biosynthesis, autoimmune thyroid disease, hematopoietic cell lineage, and type I diabetes mellitus (permutation test p-value,0.0001,Table S4).We looked more carefully at the lysosome pathway because it had the highest number of genes involved and of the recent discovery that macrophage cholesterol ester is mobilized for efflux via autophagy, with the cholesterol ester hydrolyzed by lysosomal acid lipase [7].The number of significant cholesterol regulated transcripts in this pathway was 43 and 45 in experiment 1 and 2, respectively (Table 3), with 25 in common in both experiments.GSEA allows both the few strongly differentially expressed transcripts and the many weakly differentially expressed transcripts to factor into the enrichment analysis.Thus, most of the lysosomal pathway genes in Table 3 are only modestly regulated by cholesterol loading.Nevertheless, many minor changes in gene expression in the same pathway may add up to a large overall effect in that pathway, in this instance lysosome function.
Most of the lysosome pathway genes regulated by cholesterol showed a larger fold change in macrophages from the atherosclerosis resistant AKR than the susceptible DBA/2 strain, and a systematic analysis of the strain-cholesterol loading interaction is provided below.These experimental-validated regulated transcripts include the following lysosomal acid hydrolases: 1) proteases represented by cathepsins (CtsC, CtsE and CtsZ); 2) the peptidase napsin (Napsa); 3) glycosidases (Hyal1 and Gla); 4) palmitoyl-protein thioesterase (Ppt1); 5) phosphatase (Acp2); and 6) ceramidase (Asah1).Point mutations in the ASAH1 gene leads to the lysosomal storage disorder Farber disease; and, Asah1+/2 mice develop an advanced lipid storage disease in many organs, most prominently in liver [33].The conserved set of cholesterol regulated genes also  includes: 1) several major and minor lysosomal membrane proteins (Lamp2, Cd68 and Cln5); 2) adapter-related protein complex subunits beta, mu and sigma (Ap3d1, Ap3m2 and Ap1s1) that are involved in the transport between lysosomes and Golgi; 3) sortilin 1 (Sort1); and 4) mucolipin1 (Mcoln1).Of these nonhydrolytic lysosomal pathway genes, two stood out as potentially relevant to atherosclerosis Mcoln1and Sort1.Mcoln1 encodes for a protein that co-localizes with endocytosed material that accumulates in lysosomes, and it plays a role in the exit of lipids from the lysosome and the trafficking of MHCII to the plasma membrane [34].In addition, Mcoln1-deficient neurons have defective autophagic flux leading to increased levels of LC3-II and P62 [35].Sort1 is a genome wide association study (GWAS) hit for LDL-cholesterol and coronary artery disease, and has been shown to play a role in hepatic apoB lipoprotein secretion and LDL uptake [36].
SREBP and LXR motifs in cholesterol regulated transcripts.Since GSEA for sequence motifs was not very informative for the cholesterol loading datasets (data not shown); we examined motifs for two well known sterol regulated transcription factors.The classical example of cholesterol regulation of gene expression involves the down regulation of genes containing the sterol responsive element (SRE).This regulation is mediated by sterol control of sterol regulatory element binding protein (SREBP) processing in the ER and Golgi, such that high sterols repress its processing and low sterols permit its processing into a positively acting transcription factor [37].Thus, we searched for the V$SREBP1_02 motif (KATCACCCCAC target sequence) motif within 2 kb of the transcription start site among all expressed transcripts in the two experiments.We identified 24 expressed transcripts associated with the V$SREBP1_02 motif in the first experiment.11 transcripts (46%) were significantly regulated by cholesterol loading with 6 up regulated and 5 down regulated (Table S8).In the second experiment, we identified 35 expressed transcripts associated with the V$SREBP1_02 motif, 12 of which (34%) were significantly regulated by cholesterol loading with 4 up regulated and 8 down regulated.Seven replicated transcripts with this motif met our criteria for significant regulation by cholesterol loading in both experiments, two up and five down regulated (Table S8); however, the overlap between the two sets of   cholesterol regulated SERPB1-associated transcripts was not significantly different compared to the overlap of all cholesterol regulated genes (one-tailed Fisher's exact test, p = 0.15).These seven replicated transcripts are not directly involved in cholesterol biosynthesis.A similar finding was reported in mouse macrophages deficient for SREBP1a, in which the most highly regulated transcripts did not include those coding for classical cholesterol biosynthesis genes [38].Overall, we were surprised that several of the SREBP motif-containing genes were induced upon loading, a state where SREPB is expected to be low and the SRE containing genes are repressed.However, there are now well known examples of SREBP acting as a transcriptional repressor, such as SREBP repression of IRS2 transcription in liver [39].The oxysterol activated transcription factor LXR heterodimerizes with RXR and binds to genes harboring LXR responsive elements, often leading to sterol mediated up-regulated gene expression, as demonstrated for Abca1, Abcg1, and apoE [40].We searched for the LXR DR4 motif TGACCGNNAGTRACCC within 2 kb of the start site of expressed transcripts.We identified 32 expressed transcripts associated with the LXR DR4 motif in the first experiment.17 transcripts (53%) were significantly regulated by cholesterol loading with 6 up regulated and 11 down regulated (Table S9).In the second experiment, we identified 31 expressed transcripts associated with the LXR DR4 motif, of which 11 transcripts (55%) were significantly regulated by cholesterol loading with 7 up regulated and 4 down regulated.Six transcripts with this motif met our criteria for significant regulation by cholesterol loading in both experiments, two up and four down regulated (Table S9); however, the overlap between the two sets of cholesterol regulated LXR-associated transcripts was not significantly different compared to the overlap of all cholesterol regulated genes (one-tailed Fisher's exact test, p = 0.45).The two replicated up-regulated genes were Abcg1 and Stac2.Abcg1 is a wellknown LXR target gene, but we did not identify Abca1 another well-known target, and apoE is not expressed in our apoE-deficient macrophages.The expression of Abca1 on our microarrays was very low, possibly due to poor probe design, so we quantified its expression in the experiment 1 samples by qPCR, and determined that it is indeed up-regulated by cholesterol loading in both strains (Figure 2).In a combined ChIP and gene expression study in human THP1 macrophages, LXR agonist treatment was found to lead to widespread up and down regulation of genes adjacent to confirmed LXR binding sites [41], confirming that LXR may act as either a transcriptional activator or repressor.

Cholesterol loading-strain interaction effect on BMM transcriptome
A fitted linear model using strain and loading as additive covariates was used to identify the transcripts with a significant cholesterol loading-strain interaction effect.This interactive effect identifies transcripts that have different directions or degrees of cholesterol regulation between the two strains, for example, a transcript that is up regulated in AKR and down regulated in DBA/2, or a transcript that is highly up regulated in AKR but only moderately up regulated in DBA/2 and vice-versa.
Experiment 1 samples.1,929 probes were identified with a significant loading-strain interaction effect at an FDR adjusted p-  value,0.01,with several transcripts independently identified by multiple probes (Table S10).The top 10 most significant transcripts for strain-loading interaction effect were: Slamf9, Chac1, Trib3, Vwf, Gja1, Tgfbi, Dok2, Gadd45a, Gdf15 and Dner.The response of three of these transcripts to cholesterol loading in both strains is shown in Figure 3, with Slamf9 highly down-regulated by loading in AKR (25.4 fold) and not significantly changed in DBA/2; Trib3 highly up-regulated in AKR (6.1 fold) and moderately up-regulated in DBA/2 (1.2 fold); and Dner downregulated in AKR (22.5 fold) and up-regulated in DBA/2 (1.4 fold).
Experiment 2 samples.965 probes were identified with a significant loading-strain interaction effect at an FDR adjusted pvalue,0.01,with several transcripts independently identified by multiple probes (Table S11).The top 10 most significant probes for strain-loading interaction effect were: Ddit3, Aqp9, Trib3, Nurp1, Cox6a2, Asns, Ccr5 (represented by three independent probes), and Slc1a4.The response of three of these transcripts to cholesterol loading in both strains is shown in Figure 4 Combined results and discussion.There were 213 probes with highly significant strain-loading interactions that were conserved in both experiments (presented in Table S12).Ddit3, identified as the transcript with the most significant loading-strain interaction effect in the second experiment samples, was also observed with a significant strain loading effect in the first experiment.This gene encodes the DNA-damage inducible transcript 3 proteins, more widely known as CHOP (C/EBP homologous protein).CHOP is a transcription factor whose expression is up-regulated by endoplasmic reticulum (ER) stress response and the unfolded protein response (UPR), which under prolonged stress can trigger apoptosis [42,43].Two other genes known to participate in the ER stress pathway, Trib3 and Atf4 [44][45][46], were also among the list of genes with strain-loading interactions in both experiments with similar regulation as seen for Ddit3, highly up regulated by cholesterol loading in AKR, and either down regulated or unchanged in loaded DBA/2 cells.Tabas and colleagues have previously shown that free cholesterol loading of macrophages induces ER stress including CHOP production [47,48].Thus, our finding of increased expression of ER stress related genes in AKR, the strain with increased free cholesterol loading (Figure 1), fits well with the effects of free cholesterol on ER stress described by Tabas [47,48].
ER stress, CHOP, and Trib3 have previously been implicated in atherosclerosis.In ApoE2/2 mice fed a western-type diet to induce advanced lesions, CHOP deficiency led to 35% smaller aortic lesions and 50% less plaque necrosis when compared to controls, despite similar levels of plasma lipoproteins (49).Similar results were obtained on the Ldlr-deficient background [49].Adenoviral knockdown of Trib3 in ApoE-Ldlr double knockout also led to smaller aortic lesions [50].Our finding of increased ER stress response gene expression in cholesterol loaded macrophages from the atherosclerosis resistant AKR strain appears to differ from the above findings where decreased ER stress response was associated with smaller lesions.There are several potential reasons that might explain this discrepancy.One explanation is that the atherosclerosis studies comparing AKR and DBA/2 strains were performed in chow fed mice at an early time point prior to lesion necrosis, and Tabas has postulated that lesion macrophage apoptosis in early stage atherosclerosis is associated with increased phagocytic clearance of apoptotic cells and diminished lesion progression.Another explanation may be that ER stress can also trigger autophagy, which can be protective against lesion progression [51].Autophagy also plays a role in cholesterol ester catabolism by delivering foam cell lipid droplets to the lysosome where lysosomal acid lipase cleaves cholesterol esters into free cholesterol, which is the substrate for cholesterol efflux [7].

Validation of data by quantitative Real-Time PCR (qPCR)
To confirm the microarray data we performed quantitative realtime PCR for four selected transcripts: three highly regulated ones, Trib3, Atf4 and Cln5, along with and LXR regulated transcript Abcg1.The expression levels for each transcript in the unloaded and cholesterol loaded samples from experiment 1 were quantified relative to the ubiquitin C (Ubc) control, which was found to be the least variable transcript in this study among the list of endogenous qPCR controls available from Applied Biosystems (,3% coefficient of variation among the 16 samples in experiment 1).For each of these four transcripts, the expression levels quantified by qPCR were found to be consistent with the microarray data.Linear regression analysis revealed a highly significant positive correlation for each tested transcript (R 2 values ranged from 0.66 to 0.97, Figure 5).

Western Blot Analysis
We compared Trib3 mRNA levels (microarray) and protein levels (Western blot) in unloaded and cholesterol loaded AKR and DBA/2 BMM (Figure 4C and 4D).A large increase in Trib3 upon cholesterol loading in AKR cells was observed at both the mRNA and protein levels.In contrast, upon cholesterol loading in DBA/2 cells there was no appreciable change in Trib3 mRNA or protein.Thus, the strain-specific response to cholesterol was reproducibly observed at both the mRNA and protein levels.However, the basal levels of Trib3 mRNA in the AKR and DBA/2 cells were similar, while Trib3 protein levels were much higher in DBA/2 cells.Differences observed between mRNA and protein levels are not uncommon in complex biological samples [52].Previous studies comparing mRNA and protein abundance among different mouse strains have found limited correlations (mean R = 0.27), such that the mean association between protein and mRNA levels is only 7.3% comparing various strains [53].Furthermore, this study showed that clinical traits among the various mouse strains were more strongly correlated with transcript levels than protein levels [53].

Identification of Candidate Genes for Ath22, Ath26 and Ath28 QTLs
Several studies have integrated genetics and genomics to identify plausible candidate genes for complex diseases.Previous studies in our lab have identified atherosclerosis QTLs from an AKRxDBA/2 ApoE2/2 intercross [19], of which Ath28 QTL on chromosome 2, Ath22 on chromosome 15 and Ath26 on chromosome 17 were recently validated in a second independent F 2 cohort [54].We searched the list of genes differentially  expressed by strain, cholesterol loading, or the cholesterol loadingstrain interaction, whose differential expression was conserved in both experiments.Table 4 shows this list of differentially regulated atherosclerosis modifier candidate genes that reside within these three QTL intervals each defined as 1-LOD score drop from the QTL peak.More candidates were identified within the Ath26 QTL interval on chromosome 17, as this is the largest and most genedense QTL interval.S4 Significantly represented biological pathways identified in experiment 1 and 2 datasets for strain, loading and strain-loading interaction effect by Gene Set Enrichment Analysis.

(XLS)
Table S5 Genes with a significant cholesterol loading effect in experiment 1 samples.The fold change after cholesterol loading is represented as the beta-coefficient for the fitted linear model.Positive log2 fold change, transcript expressed higher in response to cholesterol loading; negative log2 fold change, transcript expressed lower upon cholesterol loading.(XLS) Table S6 Genes with a significant cholesterol loading effect in experiment 2 samples.(XLS) Table S7 Regulated transcripts with a significant cholesterol loading effect that overlap between experiment 1 and 2 datasets.(XLS) Table S8 Significant transcripts for cholesterol loading effect, with promoter regions 2 Kb to +2 Kb from the start site of transcription containing the motif KAT-CACCCCAC, which matches annotation for SREBF1: sterol regulatory element binding transcription factor 1.

(XLS)
Table S9 Significant transcripts for cholesterol loading effect, with promoter regions 2 Kb to +2 Kb from the start site of transcription containing the motif TGACCGNNAGTRACCC, which matches annotation for NR1H3: nuclear receptor subfamily 1, group H, member 3 (DR4).

(XLS)
Table S10 Genes with a significant strain-cholesterol loading interaction effect in experiment 1. Strength of loading interaction effect on fitted model (with direction relative to expected DBA/2 response to cholesterol loading).(XLS) Table S11 Genes with a significant strain-cholesterol loading interaction effect in experiment 2.

(XLS)
Table S12 Regulated transcripts with a significant cholesterol loading-strain interaction effect that overlap between experiment 1 and 2 datasets.(XLS)

Figure 1 .
Figure 1.Total (A), esterified (B), and free (C) cholesterol mass in unloaded and AcLDL loaded AKR and DBA/2 ApoE2/2 macrophages.All data are mean6SD, N = 4 per group using the average of triplicate assays per sample.P-values were calculated by ANOVA with Newman-Keuls posttest, showing only the strain differences after loading.There were no significant strain differences in unloaded cells.doi:10.1371/journal.pone.0065003.g001 genes between the two experiments are shown in bold.#Calculated as log2 of AKR loaded/AKR unloaded average expression, positive numbers higher in loaded, negative numbers higher in loaded.*Calculatedas log2 of DBA loaded/DBA unloaded average expression doi:10.1371/journal.pone.0065003.t003

Figure 2 .
Figure 2. Gene expression and validation of microarray data by quantitative real-time PCR in experiment 1 and 2 macrophages.Abca1 (A) and Abcg1 (B) expression was calculated relative to ubiquitin C (Ubc) gene expression, an endogenous control whose expression remained unchanged under the conditions of experiment.Values are expressed as mean6SD (N = 4).Different numbers above bars show p,0.001 (A) or p,0.05 (B) by Newman-Keuls ANOVA posttest, while similar numbers above bars show no significant differences.doi:10.1371/journal.pone.0065003.g002

Figure 5 .
Figure 5. Validation of microarray expression data.Linear regression analysis of microarray expression and qPCR data for Trib3(A), Abcg1 (B), Atf4 (C) and Cln5 (D) performed in experiment 1 unloaded and loaded AKR and DBA/2 macrophages.Microarray data were not log2 transformed for this analysis.doi:10.1371/journal.pone.0065003.g005 strain effect b Transcripts significant for cholesterol loading effect c Transcripts significant for cholesterol loading-strain interaction effect # Ath28 QTL confidence interval on chromosome 2 extends from 160 to 181Mb $ Ath22 QTL confidence interval on chromosome 15 extends from 3 to 33Mb & Ath26 QTL confidence interval on chromosome 17 extends from 12 to 65Mb doi:10.1371/journal.pone.0065003.t004

Figure
Figure S1Hierarchical clustering analysis of 32 samples included in the study.Four replicates each for control and cholesterol loaded macrophages from AKR and DBA/2 strain, both independent experiments were included (loaded; unloaded samples).(PDF)FigureS2Conservation of cholesterol induced changes in macrophage gene expression in two independent experiments.Linear regression analysis of log2 fold changes of the 1,140 overlapping transcripts between experiment 1and experiment 2 dataset that are significantly regulated by cholesterol loading in one or both strains.P-value of linear regression ,0.0001 (PDF)

Table 1 .
Top 10differentially expressed transcripts between AKR and DBA/2 unloaded macrophages ranked by p-value.

Table 2 .
Top 10 significantly up-regulated or down-regulated transcripts in response to cholesterol loading ranked by p-value.
# The fold change after cholesterol loading is represented as the b-coefficient for the fitted linear model.Positive log2 fold change, transcript expressed higher in response to cholesterol loading; negative log2 fold change, transcript expressed lower upon cholesterol loading.*FDR adjusted p-value based upon permutation.doi:10.1371/journal.pone.0065003.t002

Table 3 .
Significantly regulated transcripts upon cholesterol loading involved in lysosome pathway, ranked by fold change (loaded/unloaded).

Table 4 .
Differentially expressed transcripts conserved in both experiments that reside within Ath28, Ath22 and Ath26 QTLs.

Table S1 Differentially
Expressed Transcripts between AKR and DBA in Unloaded BMM (experiment 1).Table S2 Differentially Expressed Transcripts between AKR and DBA in Unloaded BMM (experiment 2).Table S3 Differentially Expressed Transcripts between AKR and DBA/2 that Overlap between Experiment 1 and 2 Datasets in Unloaded BMM.(XLS)