Combined ChIP-Seq and transcriptome analysis identifies AP-1/JunD as a primary regulator of oxidative stress and IL-1β synthesis in macrophages

Jund regulates macrophage gene expression that controls primarily oxidative stress and IL-1β synthesis

A description of the macrophage function related to different levels of JunD expression in different inbred rat is summarised in Table 1. In order to identify genes under the transcriptional control of JunD in primary macrophages, expression levels of Jund were first silenced by RNA interference in the WKY BMDMs that over-express JunD (Figure 1) [12]. Following confirmation of knockdown by qRT-PCR and Western blot (Figure 2A and C), the samples were subjected to microarray analysis (Additional file 1: Figure S1A). Because AP1/JunD regulates lipopolysaccharide (LPS)-TLR4 mediated cytokine secretion in primary human macrophages [12, 14], microarray analysis was also performed in samples treated with Jund siRNA (48h) and stimulated with LPS for 8 h (Figure 2B and Additional file 1: Figure S1B). Genome-wide analysis of the BMDM transcriptome identified 1672 differentially expressed genes between unstimulated BMDMs transfected with Jund siRNA compared to scrambled control and 1,476 differentially expressed genes following 8 hours of LPS (100ng/ml) stimulation (FDR< 5%). JunD acted both as an activator and as a repressor of gene expression in BMDMs. Transfection of WKY BMDMs with Jund siRNA, thereby lowering Jund expression, resulted in the reduced expression of 868 genes out of 1672 (~ 50%) demonstrating that Jund had an activatory role in gene transcription. Alongside this, Jund knockdown also increased the expression of 804 genes (~ 50%) demonstrating that Jund could also have a repressive effect on transcription. After eight hours of LPS stimulation, 638 genes had reduced expression in the Jund siRNA knockdown group compared to controls and 838 genes demonstrated higher expression following Jund knockdown. Validation by qRT-PCR of a set of differentially expressed genes following Jund siRNA knock-down that encompassed a range of fold change differences between the two siRNA groups confirmed the microarray findings (Figure 3A and B, Additional file 2: Table S1) in 21 out of the 23 (91%) genes selected. The directional change (activation or repression) of the microarray data was confirmed in all 23 genes.

Strain	JunD levels in macrophages*	Macrophage activation**
WKY	+++	+++
WKY.LCrgn2	+	+
Lewis	+	+

* Jund expression levels and JunD/AP-1 protein binding measured by qRT-PCR and TransAm assay [12].

** Macrophage activation assessed by Fc receptor mediated oxidative burst [12].

JunD expression levels determine the extent of the JunD cistromes

Gene ontology analysis of the JunD-bound genes in the different strains and states identified marked similarities between the enriched terms in basal WKY BMDMs (Additional file 2: Table S3) and WKY.LCrgn2 BMDMs (Additional file 2: Table S4) as well as LPS stimulated WKY.LCrgn2 BMDMs (Additional file 2: Table S5) covering multiple autoimmune disease pathways and core macrophage functions such as antigen processing. This suggested that regardless of expression level, JunD binds to genes involved in core macrophage processes. Following LPS stimulation in WKY BMDMs enrichment was seen for immune processes linked with the response of the macrophage to stimulation such as intracellular signalling cascades and the MAPK signalling pathway (Additional file 2: Table S6) highlighting the role played by JunD in the LPS response.

Integration of microarray and ChIP-Seq datasets identifies primary JunD targets in macrophages

In order to identify the set of primary JunD targets, we integrated all the genomic datasets and identified the genes that correlated closely with the expression pattern of Jund during the LPS timecourse that were also differentially expressed following siRNA knockdown of Jund and showed a ChIP-Seq peak in WKY BMDMs. This identified two major networks (basal and LPS) with 24 genes in the basal state (Figure 5, upper panel) and 36 genes after LPS stimulation (Figure 5, lower panel) as primary JunD targets through which JunD mediates its effect on macrophage activation. The primary JunD targets correlate with Jund expression levels in two independent datasets (microarray analysis on WKY and WKY.LCrgn2 BMDM LPS time course and microarray on Jund siRNA in WKY BMDMs) and have a JunD ChIP-Seq peak in WKY BMDMs. These are genes where the expression is under direct control of JunD/AP1 binding, correlating with cellular JunD levels. The transcription factor Runx1 was identified as primary Jund targets together with several genes involved in oxidative stress such as Trpv4, Vav2, Ifi30, Nqo2, and P2ry2, Ctnnb1 and Bcl2l1 Additional file 1: Figure S5).

JunD expression levels determine active IL-1β secretion in primary macrophages and nephritic glomeruli

Since the previous role of JunD is regulating macrophage oxidative burst is known [12] and our data shows that JunD regulates genes involved in IL-1β synthesis (Il1b, Nlrp3, Klrb1a), we investigated whether JunD regulates active IL-1β secretion following Nlrp3-inflammasome activation in primary rat macrophages. Western Blot analysis of mature IL-1β in BMDMs primed with LPS and activated by ATP demonstrated a reduction in mature IL-1β between WKY and WKY.LCrgn2 BMDMs (Figure 6A). Importantly, there was a marked increase in mature IL-1β production in LEW.WCrgn2 BMDMs compared to LEW confirming that the Jund congenic interval was able to alter mature IL-1β production (Figure 6A). These results were confirmed by ELISA for IL-1β in LPS primed and ATP stimulated BMDMs which showed significant differences between the production of IL-1β by WKY BMDMs compared to WKY.LCrgn2 and LEW BMDMs (Figure 6B). The production of IL-1β in nephritic glomeruli from WKY, WKY.LCrgn2, LEW and LEW.WCrgn2 rats was examined (Figure 6C). This showed a significant reduction in IL-1β production in all the other strains compared to WKY demonstrating a role for the Crgn2 congenic interval in IL-1β production.

Animals

WKY (WKY/NCrl) and LEW (LEW/Crl) rats were purchased from Charles River (Margate, UK). Single congenic rats were generated by introgressing the Crgn2 QTL from chromosome 16 from a LEW donor onto a WKY recipient background and vice versa as previously described [12]. All procedures were performed in accordance with the United Kingdom Animals (Scientific Procedures) Act.

BMDM culture

BMDMs were prepared from the femurs of parental and congenic strains using previously described methods [12]. Femurs from adult (8-10 weeks) rats were isolated and flushed with Hanks buffer (Life Technologies). Total bone marrow derived cells were plated and cultured for 5 days in Dulbecco’s modified Eagle’s medium (Life Technologies) containing 25 mM Hepes (Sigma), 25% L929 conditioned medium, 25% decomplemented fetal bovine serum (Biosera), penicillin (100 U/ml, Invitrogen), streptomycin (100 μg/ml, Invitrogen) and _L-glutamine (2 mM Invitrogen). The cells were characterised as macrophages by ED-1 staining. Basal macrophages were left unstimulated whilst stimulated cells were stimulated with 100 ng/ml lipopolysaccharide (Sigma). Following stimulation, cells for gene expression analysis were homogenized in TRIzol (Invitrogen) and stored at -80°C.

siRNA inhibition of Jund expression

siRNA knockdown was carried out as previously described [12]. Briefly, on day 5 of culture, WKY BMDMs were replated in six-well plates (1x10⁶ cells per well) in DMEM (Invitrogen) overnight and transfected for 48 hours with siGENOME SMARTpool for Jund (100 nM, Dharmacon) or siGENOME non-targeting siRNA pool as the scrambled control siRNA using Dharmafect 1 (1:50, Dharmacon) as a transfection reagent in OPTIMEM medium (Invitrogen). The siRNA sequences used in the siGENOME SMARTpool for Jund are listed in Additional file 2: Table S7. siRNA knockdown was confirmed with quantitative PCR (detailed below) and Western blotting (Figure 1).

RNA extraction and microarray preparation

Total RNA was extracted using the TRIzol method and purified using RNeasy Plus spin columns (Qiagen). 100ng of RNA was amplified, labelled and hybridised to Rat Gene 1.0 ST arrays (Affymetrix, Santa Clara, CA, USA) using the Ambion WT Expression Kit (Life Technologies) as per manufacturer’s instructions. For timecourse expression analysis, four BMDM preparations from four biological replicates were used for each timepoint and condition. For siRNA expression analysis, four BMDM preparations from at least two biological replicates were used for each timepoint and condition. The microarray data is available in MIAME-compliant (minimum information about a microarray experiment) format at the Array Express database (http://www.ebi.ac.uk/arrayexpress) under accession code E-MEXP-3469.

Microarray data analysis

CEL intensity files were produced using GeneChip Operating Software version 1.4 (Affymetrix) and quality tested using the Affymetrix Expression Console v1.1.2. All 32 files in the timecourse data set and 16 files in the siRNA dataset were suitable for further analysis. Probe-level data was normalised using robust multichip average (RMA) [48, 49]. A custom definition file was created using up-to-date probe information [50] and filtered to exclude probes containing the 2,520,602 single nucleotide polymorphisms present between the WKY and LEW genomes (Santosh Atanur, MRC Clinical sciences centre, personal communication). The moderated T test with 40,000 permutations implemented in Statistical Analysis of Microarrays (SAM) version 3.0 was used to identify differentially expressed genes at an FDR threshold of 5% and timecourse analysis was performed using EDGE with 40,000 permutations and a 5% FDR threshold [51]. Hierarchical clustering analysis was performed using MultiExperiment Viewer (MeV) v4.8 [52, 53] with the Euclidean distance measure. Gene ontology analysis was carried out using the functional annotation tools within DAVID, the Database for Annotation, Visualisation and Integrated Discovery v6.7 [54, 55].

Quantitative PCR

All qPCRs were performed with an ABI 7900 Sequence Detection System (Applied Biosystems, Warrington, UK). A two-step protocol was used as previously described [56] beginning with cDNA synthesis with iScript select (Bio-Rad) followed by PCR using SYBR Green Jumpstart Taq Ready Mix (Sigma). A total of 10ng of cDNA per sample was used. All samples were amplified using a set of 4 biological replicates with three technical replicates used per sample in the PCR. Sequence detection software (SDS) version (Applied Biosystems) was used to obtain the Ct values. Results were analysed using the comparative Ct method and each sample was normalised to the reference gene Hprt, to account for any cDNA loading differences. The primer sequences used for the qRT-PCR validation of microarray data are listed in Additional file 2: Table S8.

Chromatin immunoprecipitation (ChIP)

BMDMs were left in the basal condition or stimulated for 2 hours with 100ng/ml LPS. Cells were fixed for 10 minutes with 1% formaldehyde and ChIP performed using ChIP-IT Express as per manufacturer’s instructions with some modifications. Sonication was carried out using Covaris S2 (Woburn, Massachusetts, USA) in a volume of 300μl sonication buffer with the following settings; 20% duty cycle, intensity 8, 200 cycles/burst, cycle length 30 s for 28 to 30 cycles dependant on cell count. The ChIP lysate was immunoprecipitated using 2 μg of JunD antibody (Santa Cruz sc74-X) or negative IgG control (sc-2026) overnight. Cross links in the immunoprecipitated chromatin and control input chromatin were reversed by heating the samples at 65°C for 5 hours followed by proteinase K digestion for 1 hour. Samples were purified using Qiagen MinElute columns as per manufacturer’s instructions prior to downstream analysis.

High throughput sequencing

Single read library preparation and high throughput single read sequencing for 36 cycles was carried out on an Illumina Genome Analyser IIx according to the manufacturer’s protocols (Illumina) with some modifications. Each immunoprecipitated sample for library preparation was the product of five separate technical replicates of immunoprecipitation pooled together and purified using Qiagen MinElute columns and quantified using the hsDNA Qubit assay (Invitrogen, Life Technologies, California, USA). A 1:20 dilution of Illumina adaptors was used at the adaptor ligation step to avoid adaptor dimer formation and a 1:2 primer dilution used to prevent dimerization during the PCR amplification stage. The samples were quantified using the hsDNA Qubit assay (Invitrogen, Life Technologies, California, USA) and the size range analysed using a HS DNA chip on a 2100 Bioanalyser (Agilent Technologies, West Lothian, UK) prior to submission for qPCR analysis and cluster generation and sequencing by the CSC/IC Genome Core facility.

qPCR of ChIP enriched DNA

Immunoprecipitated DNA fragments were analysed by real-time PCR. Primers used are listed in Additional file 2: Table S9. All samples were amplified using a set of 3 biological replicates with three technical replicates used per sample. After an initial denaturation step of 94°C for 2mins, the samples were cycled 40 times at 94°C for 15s, 60°C for 1 min and 72°C for 1 min with data collection performed during the 72°C elongation step. Sequence detection software (SDS) version 2.3 (Applied Biosystems) was used to obtain the Ct values and each sample was analysed with reference to Ct values for matched control ‘Input’ non-immunoprecipitated chromatin. A standard curve of 1:5 dilutions of Input DNA was used to calculate the % Input level of the transcription factor binding at the investigated locus. The standard curve was constructed from the Input DNA sample for the appropriate strain and condition and analysed within the SDS software v2.3. For each test gene the % Input levels were then determined using %total = 2^ΔCt × (% of input sample used) where Δ = Ct (Input) - Ct (sample). Fold change over IgG was expressed using 2^-ΔCt where ΔCt= ΔCt_Jund – ΔCt_IgG.

ChIP-Seq data analysis

Sequencing of the ChIP-Seq libraries was carried out on the high throughput Illumina Genome Analyzer II. Initial data processing was performed using Illumina Real Time Analysis (RTA) v1.6.32 software (equivalent to Illumina Consensus Assessment of Sequence and Variation, CASAVA 1.6) with default filter and quality settings. Quality filtered reads were then realigned to the reference rat genome (RGSC3.4) using the Burrows Wheeler Alignment tool v0.5.9 (BWA) [57]. Read ends were trimmed if Phred-scaled base quality scores dropped below 20. Reads that uniquely mapped to the reference genome were used to detect areas of enrichment with BayesPeak v1.1.3 [22, 23] using a posterior probability threshold of 0.9. A stringent posterior probability threshold of 0.9 was used to filter all bins passing the threshold to form the final contiguous peak regions to produce a more accurate reflection of true peak calls. An over fitting diagnostic was performed using λ₁ < 0.7 and score < -2.25 to filter out regions which showed no enrichment but had a high enough background for the algorithm to call peaks in. Peak regions were annotated using the gene intervals annotator (GIN) implemented in CARPET using a gene priority approach to give the associated transcript I.D., the associated gene feature and the distance of the peak to the nearest transcriptional start site [58]. HOMER was used to predict motif occurrence within peaks [32, 59] with default settings for a maximum motif length of 12 base pairs. The outputs of the peak calling algorithms were visualised in the Integrative Genomics Viewer [60] using custom WIG files generated from the output data generated by BayesPeak. The latter were generated by extending each mapped read by 200 bp and then by using 10 bp bins the overlapping tag count was generated.

Integrated data analysis and identification of primary JunD targets

JunD gene expression patterns in WKY and WKY.LCrgn2 BMDMs over the LPS time course were used for Spearman correlation analysis with the rest of the transcripts on the microarrays (P < 0.001 cut off). In order to integrate the three different datasets, the list of significantly correlated set of transcripts was filtered and annotated according to two criteria: significant differential expression following Jund siRNA knockdown (both basal and after LPS stimulation) and the presence of the JunD ChIP-Seq peaks in WKY BMDMs (both basal and after LPS stimulation). Two separate networks (basal and LPS) were built with Cytoscape version 2.8.3 showing primary JunD targets.

Isolation and culture of rat nephritic glomeruli and detection of IL-1β by Western Blot and ELISA

Glomeruli were isolated as previously described [13]. After 48 hours of incubation at 37°C, nephritic glomeruli supernatants were collected and stored at -20°C for IL-1β sandwich ELISA analysis. For IL-1β detection in BMDMs by Western Blotting, cells were primed with LPS (1 μg/ml, 3hours) and stimulated with ATP (5 mM) for 30 minutes. BMDMs were then scraped and both cells and supernatant were collected, filtered using Amicon ultra centrifugal filters for protein purification and concentration (10 kDa cut-off, Millipore, UK) and the concentrated samples were diluted with 5× sample buffer containing 200 mM Tris-HCl, 6% SDS, 2mM EDTA, 4% 2- Mercaptoethanol, 10% glycerol and boiled for 10 minutes. The samples were then resolved by SDS polyacrylamide gel electrophoresis (PAGE) and transferred to an Immobilon-P Transfer Membrane (Millipore). Rabbit polyclonal anti IL-1β (New England BioLabs, UK) was used to detect the mature IL-1β. To assess secreted IL-1β in BMDMs and nephritic glomeruli, cell supernatants were subjected to sandwich ELISA and secreted IL-1β amounts were determined using a standard curve with rat recombinant IL-1β according to manufacturer’s instructions (R and D Systems).

For the confirmation of siRNA knockdown of Jund, WKY BMDMs were plated into six-well plates at a density of 1 × 10⁶ cells per well and treated with Jund specific siRNA or a scrambled oligonucleotide for 48 hours before total protein was extracted for Western Blot analysis using the above technique. For JunD detection a specific rabbit polyclonal anti-Jund antibody from Santa Cruz Biotechnology (USA) was used. This blot is representative of 4 different experiments performed with 4 biological replicates.