The genetic basis for differential macrophage response to infectious and inflammatory stimuli

The differences in response to infectious and inflammatory stimuli in AJ and B6 mice have a genetic component [19, 38, 43]. Therefore, having established differences in parasite growth and activation phenotypes in AJ and B6 BMDM, we sought to establish the genetic basis for the differences in macrophage activation and toxoplasmacidal activity between AJ and B6 using the AXB/BXA mice. We stimulated BMDM obtained from 26 age-matched female AXB/BXA RI mice with IFNG+TNF, IL-4, CpG, or LPS and measured the amount of NO, urea, IL-12, IL-10, and CCL22 produced. We also infected IFNG+TNF stimulated BMDM with Toxoplasma that express luciferase and measured relative parasite growth by measuring luciferase activity. These phenotypes exhibited a continuous distribution in the RI mice, which is characteristic of quantitative traits (S1A–S1F Fig).

Due to the stable and unique combination of blocks of parental alleles in their genomes, the AXB/BXA RI mice are particularly suited for quantitative trait locus (QTL) mapping [55]. Therefore, we used the AXB/BXA genetic map (containing 934 informative genetic markers [56]) in a genome-wide scan in R/qtl [57] to identify the genomic regions that modulate differences in AJ and B6 BMDM phenotypes. To correct for multiple testing on the 934 genetic markers, we performed 1000 permutation tests on the individuals relative to their phenotypes [58] to obtain adjusted p-values for each QTL (Table 1). Parasite growth in the IFNG+TNF-stimulated BMDM mapped to chromosome 3 (147.7 Mb), proximal to the Guanylate binding protein (Gbp) locus (142.6 Mb) that is known to modulate the murine macrophage response to Toxoplasma [59–62]. Even though AJ and B6 BMDM display distinct polarization states following IFNG+TNF or IL-4 stimulation, except for IL-12 and CCL22, we did not observe statistically significant QTL peaks for any of the activation phenotypes. Instead we observed a suggestive QTL for NO on chromosome 12 (Table 1) and a second QTL for NO on chromosome 4 (S2A Fig). The additive QTL for NO also mapped to chromosome 4 (S2B Fig). Similar to observations in the parental BMDM, the AJ allele at the chromosome 12 QTL was associated with high amounts of NO (S2C Fig). However, the AJ allele at the chromosome 4 QTL was associated with low amounts of NO (S2D Fig). Although mapping in cis to Arg1 on chromosome 10, the QTL for the amount of urea did not reach statistical significance. Nevertheless, the allele effect at the suggestive urea QTL was consistent with the parental allele effect on urea.

Finally, we selected the top 2 QTL for each phenotype (where there was marginal difference between the LOD scores for the second and third largest QTL, we picked both) and grouped the BMDM based on the genotypes at each QTL. We then estimated the QTL inheritance pattern by comparing the values for the corresponding phenotype amongst the genotypes using a one-way ANOVA with Tukey’s-Post-test (StatPlus, AnalystSoft Inc) [63] (Table 2).

Divergent genomic loci modulate macrophage transcriptional response to infectious and inflammatory stimuli

Discrete transcriptional programs modulate the response of immune cells, including macrophages, to stimuli such as pathogens and immune factors [24, 31, 33]. As such transcriptional profiles can be used to gain insight into the intricate and incipient molecular networks that modulate complex traits, such as macrophage activation [24, 34, 64, 65]. Consequently, we investigated whether specific transcriptional programs modulate the differential activation of AJ and B6 BMDM. To do this, we performed high throughput RNA-sequencing (RNA-seq) on BMDM obtained from the same 26 age- and sex-matched AXB/BXA mice described above and their progenitors (28 samples in total), before (resting, controls) and after infection with Toxoplasma or stimulation with IFNG+TNF, or CpG. For each sample we generated at least 100 million paired-end reads, except for the IFNG+TNF-stimulated BMDM samples that were sequenced on a single end. Although all the samples were sequenced once, due to the unique recombination of the parental alleles and homozygosity at each of the informative genetic markers (934), there are at least 4 replicates (the minimum number of mice having the same allele at each marker) for each marker. Thereafter, we processed the RNA-seq data as previously described [25]. Briefly, we aligned the RNA-seq reads to the mouse reference genome (NCBI build GRcm38, downloaded from Illumina iGenomes; https://support.illumina.com/sequencing/sequencing_software/igenome.html) using TopHat [66]. To avoid read alignment bias due to sequence polymorphisms between AJ and B6 genomes, we made a synthetic reference genome in which all the polymorphic nucleotides between AJ and B6 were converted to a neutral nucleotide [67]. However, and consistent with a recent report [68], allele bias did not significantly affect read alignment to the genome. On average, about 70% of reads in each sample uniquely mapped to the synthetic genome, which was about 1% less than the number of reads uniquely aligned to the iGenome. Henceforth, unless otherwise stated, all RNA-seq data presented herein were processed using the synthetic genome. Transcript abundance was estimated using Cufflinks [69] and reported as fragment per kilobase exon per million reads (FPKM).

Next, using the FPKM values from each of the 26 RI BMDM and the corresponding AXB/BXA genetic map, we mapped the genomic loci that modulate gene expression (expression QTL, eQTL) [32] using R/qtl [57]. As previously described [25], we performed 1000 permutations to correct for multiple testing across the 934 informative genetic markers in the AXB/BXA cross. Next, we used a false discovery rate (FDR) ≤ 10%, calculated in the qvalue package [70], to correct for multiple testing on the transcripts and to nominate significant eQTL. Finally, to allow for meaningful comparisons, we only included in the downstream analyses, for each stimulation condition, eQTL for transcripts with an average FPKM ≥ 5 across the 26 RI BMDM. The linkage analyses and the subsequent filtering steps were performed separately in resting, IFNG+TNF-stimulated, CpG-stimulated, and Toxoplasma-infected BMDM and identified, 131, 367, 688 and 1008 significant eQTL, respectively, dispersed throughout the genome (Fig 2A and 2B and S1A–S1D Dataset). Thus, consistent with previous studies [24, 71], the genetic background influences macrophage gene expression profile. Importantly the different stimuli induced transcriptional programs that were modulated by distinct eQTL hotspots, which can be exploited to unravel the complex molecular networks that regulate macrophage activation states.

Open in a separate window

Fig 2

The transcriptional response in BMDM is regulated by stimulation-specific trans-loci.

Expression quantitative trait loci (eQTL) in the RI mice were mapped in A) non-stimulated, B) IFNG+TNF-stimulated, C) Toxoplasma-infected, and D) CpG-Stimulated BMDM. Each dot represents a single eQTL (transcript). Significant eQTL located ≤ 10 Mb from the start of the physical location of the corresponding gene were designated as cis mapping (diagonal lines). All other eQTL were designated as trans-mapping (vertical lines). eQTL significance was calculated after 1000 permutations and reported at genome-wide thresholds corresponding to FDR ≤ 10%. Red spots identify genes mapping to a trans-eQTL hotpot (trans-band).

QTL hotspots are enriched for biological processes and putative regulators that modulate the expression of multiple genes

Variable gene expression can be due to sequence or structural variations close to (cis) or further removed from (trans) the gene itself, such as polymorphism in the promoter regions or at a distal transcription factor, respectively. Consequently, relative to the physical location of the corresponding gene, an eQTL can be categorized as either cis or trans. Thus, we designated eQTL that co-localized within a 10 Mb genomic window with the corresponding gene as cis and all other eQTL as trans [72]. Except for the Toxoplasma-infected BMDM (407 cis vs. 601 trans), most of the eQTL were located in cis in: resting (99 cis vs. 32 trans), IFNG+TNF-stimulated (194 cis vs. 173 trans), and CpG-stimulated (482 cis vs. 206 trans) BMDM (S1A–S1D Dataset).

Suppose a common locus was to modulate the expression of multiple genes in trans, then in linkage analysis, we should expect the eQTL for these genes to co-localize in the vicinity of the common locus forming a trans-eQTL hotspot (trans-band). Indeed, similar to previous studies [23, 24, 73], and indicative of a common variant regulating the expression of multiple genes in trans, we detected trans-eQTL hotspots, within a 10 Mb window, in all the samples (Fig 2A and 2B and S1A–S1D Dataset). Because it is possible for eQTL to co-localize by chance alone, we used Bonferroni-corrected p-values and Poisson distribution to compute the number of trans-eQTL that can co-localize in a 10 Mb genomic window by chance. Using these cutoffs, we identified 2, 3, 5, and 15 trans-eQTL hotspots in the resting, IFNG+TNF-stimulated, CpG-stimulated, and Toxoplasma-infected BMDM, respectively (S1A–S1D Dataset).

Previously, it was reported that eQTL that localize close to the physical location of the relevant gene (cis-eQTL) are a consequence of single nucleotide polymorphisms (SNPs) [72] or larger structural variants (SV), such as insertions and deletions. Thus, we investigated the nucleotide sequence in a 2000 bp window upstream and downstream of the transcription start site (TSS) of all significant eQTL. Consistent with these reports, we found that most cis-eQTL were associated with genes reported to have structural or sequence variations within 2000 bp upstream or downstream of their TSS [74, 75]. That is, in control BMDM we found 89 out of 99 (Hypergeometric test P≤2.3e-7); in the IFNG+TNF-stimulated BMDM we found 156 out of 196 (Hypergeometric test P≤3.3e-10); in the CpG-stimulated BMDM we found 422 out of 482 (Hypergeometric test P≤3.8e-23); and in the Toxoplasma-infected BMDM we found 369 out of 407 (Hypergeometric test P≤1.6e-49) cis-eQTL with polymorphisms within 2000bp upstream or downstream of TSS. These included genes with known immunological functions such as Gbp1, Gbp2, Irak4 and Srebf1.

As indicated above, trans-eQTL hotspots can be due to transcriptional regulation of several genes by a common genetic variant [76], such as a polymorphic or differentially expressed transcription factor, enhancer or repressor. Alternatively, a trans-band can be a result of a differentially expressed or polymorphic signaling protein, such as a cell surface receptor, that could lead to differential activation of transcription factors and the genes downstream of these transcription factors. Consequently, the eQTL localizing at the trans-eQTL hotspot are likely to be enriched for binding sites for a common transcription factor(s) physically located at the trans-eQTL-hotspot, biological function, or signaling pathway. Therefore, we used gene ontology (GO) [77] and rVISTA [78] to functionally characterize and search for transcription factor binding site enrichment in each of the eQTL in the different trans-hotspots (Table 3). We found that the three largest trans-bands in the Toxoplasma infected BMDM mapped to loci containing genes with either known or putative roles in macrophage inflammatory or metabolic processes, which are important host responses against intracellular pathogens. For instance the trans-band on chromosome 13 (117.1–127.8 Mb) overlapped docking protein 3 (Dok3), chemokine (C-X-C motif) ligand 14 (Cxcl14), and AU RNA binding protein/enoyl-coenzyme A hydratase (Auh), which are known to regulate various aspects of cellular inflammatory and metabolic processes [79–81]. Indeed, this trans-band was enriched for, among others, “natural killer cell mediated immunity”. Similarly, the trans-band in the IFNG+TNF-stimulated BMDM on chromosome 15 (89.6–96.6 Mb), which was enriched for “regulation of leukocyte mediated immunity”, overlapped the interleukin-1 receptor-associated kinase 4 (Irak4), known to regulate the immune response to a variety of infectious and inflammatory stimuli [18, 82, 83]. Furthermore, Irak4 contains genetic insertions and deletions in AJ relative to the reference B6 mouse strain [84]. Thus, the trans-bands are functionally enriched in biologically relevant processes and can potentially reveal novel regulators and insight in the complex gene interactions that modulate macrophage response to exogenous stimuli.

To identify putative regulators for the trans-bands, which can be variable transcription factors or signal transducers, we identified genes that were expressed (average FPKM ≥5), were physically located within 10 Mb on either side of the trans-eQTL hotspot, and were differentially expressed or had non-synonymous (NS) polymorphisms in AJ versus B6 mice. As an example, the chromosome 15 trans-band (between 79.7–99.7 Mb) contained 81 expressed coding and non-coding genes, 5 of which contained non-synonymous SNPs and 8 exhibited differential expression that mapped in cis. Of these genes, we considered Plxnb2, Irak4, and Apobec3 to be good candidates since they are known to be involved in immune response pathways [82, 85, 86], similar to the functional enrichment observed for the chromosome 15 trans-band. Additionally, these genes are polymorphic in AJ compared to B6. Due to insertions and deletions [87] and its function as an immune signaling adaptor, we considered Irak4 to be a strong candidate regulator for the chromosome 15 trans-band. Hence, we used shRNA to knockdown Irak4 in the IFNG+TNF-stimulated macrophages. Enrichment analysis on the perturbed genes following Irak4 knockdown revealed an overrepresentation (p = 0.005) of several chromosome 15 trans-eQTL (Table 4). shRNA-knockdown of the other putative candidate genes identified in this study did not perturbed most of the chromosome 15 trans-band eQTL (S1E Dataset), indicative of a specific effect of Irak4 knockdown on this trans-band.

Primary bone marrow derived macrophages and parasites

Bone marrow-derived macrophages (BMDM) were obtained from 6–8 weeks old AJ, C57BL/6J and 26 female AXB/BXA recombinant inbred mice (Jackson Laboratories). Marrow cells were obtained from each mouse by flushing the femur and tibia with cold phosphate buffered saline (PBS; GIBCO-Invitrogen). The cells were then centrifuged at 500 x g for 5 minutes at 4°C and re-suspended in 4 ml red cell lysis buffer (Sigma) and incubated on ice for 5 minutes. Next, the cells were passed through a 70 μm cell strainer (BD Biosciences) and centrifuged at 500 x g for 5 minutes at 4 °C. The cells from each mouse were subsequently grown on four 10 cm non-tissue culture petri dishes (Corning) in Dulbecco's modified Eagle's medium (DMEM; GIBCO-Invitrogen) supplemented with 10% heat-inactivated fetal bovine serum (FBS; HyClone), 2 mM L-glutamine, 1 mM sodium pyruvate, 1X MEM nonessential amino acids, and 50 μg/ml each of penicillin and streptomycin. To differentiate the cells into macrophages, the DMEM was conditioned with 20% L929 cell supernatant (containing GM-CSF (40 ng/ml), hereafter 20% L929). After incubating the cells at 37°C and 5% CO₂ for 3 days, the non-adherent cells were pipetted into 50 ml tubes and centrifuged at 500 x g for 5 minutes at 4°C and seeded in new 10 cm petri dishes. Simultaneously, the old 10 cm petri dishes were topped up with fresh media supplemented with 20% L929 i.e. in the end, for each mouse there were 8 petri dishes. After further incubation at 37°C and 5% CO₂ for 4 days, the BMDM were harvested and stored in liquid nitrogen (5 million cells/ aliquot). The BMDM yield for each mouse strain ranged between 60–150 million cells, with the B6 mouse consistently producing more cells. This protocol has previously been shown to yield pure (>99%) macrophages [12]. A Pru (type II) Toxoplasma gondii strain engineered to express firefly luciferase and GFP (Pru ΔHXGPRT A7) [109], maintained in the laboratory by serial passage on Human Foreskin Fibroblasts (HFF), was used for all infections.

Immortalized macrophages

To immortalize macrophages, we used J2 recombinant retrovirus [110] produced from ψCREJ2 cells (a generous gift from John MacMicking, Yale University School of Medicine). The J2-expressing cells were grown to confluency in DMEM medium supplemented with 10% FBS (D10). The medium containing retroviral particles was collected and passed through 0.45 μM low protein-binding filters (Millipore). In parallel, 5 x 10⁶ primary BMDM from AJ and C57BL/6J, obtained as described above, were thawed and grown for 2 days in DMEM medium supplemented with 10% FBS and 20% L929. After 2 days, the medium was replaced with the filtered J2-retrovirus-containing medium supplemented with 50% L929. After 24 hrs, the media was replaced with fresh D10 medium supplemented with 25% L929. Media was subsequently changed after every 24hrs with concomitant reduction in L929 until 10% L929 concentration was reached. Immortalized cells were harvested and stored in liquid nitrogen until use.

In vitro measurements

Unless otherwise stated, a total of 10⁴ /well immortalized BMDM (iBMDM) or 10⁵ /well primary BMDM were used in all the in vitro assays. Unless stated otherwise, before stimulation or infection, the primary BMDM or iBMDM were seeded overnight in D10 supplemented with 20% L929. Cellular nitric oxide levels were measured using the Griess reagent procedure on supernatant from non-stimulated or stimulated cells. Arginase activity was measured by quantifying the amount of urea as previously described [111]. Briefly, L-arginine was added to cell lysates and incubated at 37°C. After 1 hour, 175μl of an acid mixture containing sulfuric acid/phosphoric acid/water (H₂SO₄/H₃PO₄/H₂O) in a 1:3:7 ratio, was added to each well to stop the enzymatic reaction. Urea was quantified calorimetrically at 540 nm after adding 1.25 μl of 1-phenyl-1,2-propanedione-2-oxime (ISPF) and heating at 95°C for 30–60 minutes. This procedure abrogates interference from other metabolites generated, such as L-citrulline [111]. IL10, and IL12 were measured on the relevant cell supernatants using ELISA kits as previously described [12]. For parasite growth assay, cells were either left unstimulated (control) or stimulated with recombinant mouse IFNG (100 ng/ml, Peprotech) and TNF (100 ng/ml, AbD serotec) for ~18 hr. The supernatant was removed for nitric oxide assay and replaced with D10 containing Toxoplasma at an MOI ~1. The parasites were allowed to infect and replicate for 24 hrs before luciferase activity was measured using a luciferase assay kit (Promega) according to the manufacturer recommendations.

RNA sequencing

Primary BMDM were plated (3 x10⁶) overnight before stimulation or infection. For the stimulated samples, IFNG (100 ng/ml) and TNF (100 ng/ml) were added to each well for 18 hrs, while for the infected samples, a type II strain of Toxoplasma (Pru) was added to the confluent BMDM at an MOI of 1.3 for 8 hrs. Total RNA (Qiagen RNeasy Plus kit) was then isolated from the non-stimulated and non-infected cells (controls plated overnight), stimulated, and infected cells and the integrity, size, and concentration of RNA checked (Agilent 2100 Bioanalyser). The mRNA was then purified by polyA-tail enrichment (Dynabeads mRNA Purification Kit; Invitrogen), fragmented into 200–400 base-pairs, and reverse transcribed into cDNA before Illumina sequencing adapters were added to each end. Libraries were barcoded, multiplexed into 4 samples per sequencing lane in the Illumina HiSeq 2000, and sequenced from both ends resulting in 40 bp reads after discarding the barcodes. Our preliminary RNA-seq experiments with infected BMDM have shown that with 4 samples per lane, we still obtain enough read coverage for reliable gene expression analysis.

Gene expression measurements

Reads were initially mapped to the mouse genome (mm9) and the Toxoplasma (ME49 v8.2) genome using Bowtie (2.0.2) [112] and Tophat (v2.0.4) [66]. We then estimated gene expression levels in cufflinks (v2.0.0) [113] using the Illumina iGenomes refseq genome annotation (NCBI build 37.2) with the multi-read, compatible-hits corrections and upper quantile normalization options enabled. Because the reference genome to which we mapped the RNA-seq reads is based on the C57BL/6J genomic sequence, and due to the known polymorphisms between the AJ and C57BL/6J, we suspected that biases introduced at the read mapping stage might affect our expression results. To mitigate this potential bias towards the reference allele, we created a copy of the mouse genome in which all the known single nucleotide polymorphisms (SNPs) between AJ and B6, as annotated by the Wellcome Trust Sanger Institute sequencing (ftp://ftp-mouse.sanger.ac.uk/current_snps/), were converted to a third (neutral) nucleotide that is different from both the reference and AJ allele [67]. However, this did not substantially change the average proportion of uniquely mapped reads or expression profiles of individual genes in all the samples. In the end we used the mapping data generated from the synthetic genome to quantify gene expression levels.

Quantitative Trait Locus (QTL) mapping

To map QTL, we used 934 AXB/BXA genetic informative markers obtained from http://www.genenetwork.org. For all the in vitro measurements and gene expression linkage analysis, a genome-wide scan was performed using R/qtl [57]. Significance of QTL logarithm-of-odds (LOD) scores was assessed using 1000 permutations of the phenotype data [114] and the corresponding p-values reported. For the cellular phenotypes, QTL significance was reported at a genome-wide threshold corresponding to p < 0.05. However, for eQTL mapping, we further corrected for multiple testing on the multiple transcripts by using the p-values to estimate false discovery rate (FDR) in the qvalue package [115] and reported significant eQTL at FDR ≤ 10%. To identify cis- and trans- eQTL, we computed the distance from the position of the eQTL and the start of the physical location of the corresponding gene and designated any eQTL located <10 Mb from the corresponding gene as cis, otherwise trans eQTL. The procedure used to determine trans-bands has previously been described [25].

Validation of candidate regulators

We used shRNA to probe for functional or regulatory significance of some of the candidates identified in our analysis. To do this, we used C57BL/6J immortalized bone marrow-derived macrophages, described above. One day after plating, we added shRNA constructs containing a puromycin resistance marker (RNAi Platform, Broad Institute) in the presence of 8 μg/ml polybrene, followed by centrifugation at 800 x g for 2 hrs at 37°C. At the end of the spinfection, the cells were incubated for an additional 24 hrs at 37°C in 5% CO₂. The cells were then grown in fresh cell culture medium for an additional 24 hrs before adding 4 μg/ml puromycin. Transcript knockdown was measured by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) using the KAPA SYBR FAST Universal 2X qPCR Master Mix (KAPA Biosystems) on a LightCycler qPCR instrument (Roche). Fold knockdown was measured using the 2 delta-delta method [116] relative to LacZ-shRNA transduced cells. The puromycin-selected cells were either left stimulated or stimulated with IFNG+TNF, as described above, and the cell supernatant collected for nitric oxide assay.

The authors wish to thank the MIT BioMicroCenter for Illumina library preparation and sequencing. We also thank Daniel Gold and Markus Meissner for discussions and comments. We would also like to thank John MacMicking and Bae-Hoon Kim for providing the ψCREJ2 cells and the protocol for immortalizing the macrophages.