RETRACTED ARTICLE: Selected TLR7/8 agonist and type I interferon (IFN-α) cooperatively redefine the microglia transcriptome.

Introduction

Microglial cells, the resident immune sentinels of the central nervous system (CNS), are thought to participate in the pathogenesis of many neurological disorders (Crotti and Ransohoff 2016). These cells express multiple Toll-like receptors (TLRs), which are a phylogenetically conserved diverse family of sensors that drive innate immune responses following interactions with pathogen-associated molecular pattern (PAMPs), including TLR3 (viral double-stranded RNA (dsRNA)), TLR4 (lipopolysaccharide (LPS)), TLR7/8 (single-stranded RNA (ssRNA)) and TLR9 (DNA) (Glass et al. 2010; Holtman et al. 2015; Crotti and Ransohoff 2016). The interaction of TLRs with PAMPs, such as TLR7/8 and ssRNA, triggers several signaling cascades that lead to the induction of numerous target genes involved in inflammation (O’Neill et al. 2013). Among the TLRs, TLR7/8 is one of the main receptors expressed by glial cells in response to ssRNA (Butchi et al. 2008; Butchi et al. 2011). Thus, several TLR7/8 agonists including R848 (Resiquimod), an imidazoquinoline, have been discovered, and have been used to study not only immune cell (e.g., macrophage) activation but also molecular signatures associated with TLR7/TLR8-dependent signaling pathway (Butchi et al. 2008; Butchi et al. 2011; O’Neill et al. 2013).

It has been evident as after ligand recognition by TLR7/TLR8, it activates the intrinsic signaling pathways and induces type I interferon (e.g., IFN-α) to mediate innate immune responses (Uematsu and Akira 2007). IFN-α is also a pleiotropic cytokine that can either augment or suppress the expression of genes related to chronic infections and multiple sclerosis (Ivashkiv and Donlin 2014), suggesting both stimuli are able to activate different pathways and that these pathways are not redundant. Therefore, microglial cells exposed to either R848, IFN-α alone or combination can be an efficient system to investigate neuroinflammatory conditions. The interaction of IFN-α and R848 triggers synergistic IL-6 production in mouse monocyte-derived dendritic cell (DC), as measured by ELISA (Kreutz et al. 2015). However, the expression levels of a large group of immune genes using genome-wide approaches (e.g. RNA sequencing) and the molecular mechanisms underlying this crosstalk between IFN-α and TLR7/8 responses in murine PM are largely unaddressed.

Since microglial activation in vitro constitutes a valuable tool to investigate their role in normal and pathological conditions and to model neuroinflammation, here, we aim to elucidate the effects of R848, IFN-α and R848 with IFN-α on PM. To analyze the expression pattern of genes affected by these stimuli, we performed RNA sequencing (RNA-seq) analysis in primary microglial cells stimulated with R848, IFN-α or R848 with IFN-α. RNA-seq is increasingly used to determine gene expression, as it provides an unbiased digital readout and improved detection at the extremes of the transcriptome of any mammalian cell and is extremely accurate compared to microarrays (Ozsolak and Milos 2011). Here, we used RNA-seq to detect the genes that are significantly up- or downregulated according to pairwise comparisons, termed differentially expressed genes (DEGs), in murine primary microglia treated with R848, IFN-α, or R848 with IFN-α. The outcomes of these studies allowed us to identify microglial transcriptional signatures for R848, IFN-α, or R848 with IFN-α. We also identified trans-regulatory elements (e.g., altered transcription factor expression, activation or motif specificity) that may drive distinct gene expression in R848-, IFN-α- or R848 with IFN-α-treated PM. The outcome of these studies confirmed that R848, IFN-α and R848 with IFN-α stimuli generates different gene expression patterns and can constitute useful tools to study neuroinflammation.

Materials and methods

Isolation, culture and stimulation of primary microglial cells

Primary microglial (PM) cells were isolated from 3-day-old ICR mice as previously described (Witting and Moller 2011) with minor modifications. All experiments were performed in accordance with Institutional Animal Care and Use Committee (IACUC) guidelines and approved by the IACUC committee of Chung-Ang University (IACUC Number: 2016-00009). Briefly, whole brains of neonatal mice were dissected out of the skull, and blood vessels and meninges were carefully removed. Then the tissues from whole brains of 12 mice were pooled together, finely minced with sterile surgical blade, and digested using a Neural Tissue Dissociation Kit-Postnatal Neurons (Miltenyi Biotec, Germany, 130-094-802). Next, the digested cells were passed through a 70-μm nylon cell strainer (BD Bioscience, Franklin Lakes, NJ) and seeded in l-lysine-coated T-75 flasks in DMEM/nutrient 122 mixture F-12 (DMEM/F12, 1:1) containing 20% FBS (catalog # 26140; Gibco, Waltham, MA), 100 IU/ml penicillin and 10 μg/ml streptomycin (catalog # 15140) obtained from Invitrogen (Waltham, MA). The cells were maintained in a humidified incubator with a 95% air and 5% CO2 atmosphere at 37 °C. The medium was changed every 2–3 days. After 2 weeks of culture, the mixed glial cell cultures were shaken on incubating shaker at 110 rpm at 37 °C for 55 min, wherein microglial cells detached from flasks. Then the detached microglial cell suspension was collected and seeded on poly-l-lysine-coated cell culture plates. Next, microglial cells were sub-plated and used for further experiments. Microglial cells attach to the culture plate bottom much more efficiently than oligodendrocytes. At 2 h after seeding before addition of fresh media, aspiration of the old media removes unattached contaminating oligodendrocytes. More than 92% of cells obtained were PM as quantified by CD11b (rat monoclonal immunoglobulin G2b (IgG2b), clone: M1/70.15.11.5, Miltenyi Biotec Germany) FACS analysis (Suppl. Fig. 1a). We also used an anti-rat secondary antibody to avoid auto-immunofluorescent labeling (also known as background staining) in FACS analysis. Primary microglial cells were plated in a humidified incubator with 95% air and a 5% CO2 atmosphere at 37 °C for 24 h before stimulation. The cells were treated with 1 μM R848 (Resiquimod) (Hemmi et al. 2002), 100 U/ml IFN-α (Zimmermann et al. 2016) and their combination for the specified times under normal culture conditions. The cells were stimulated after diluting the stock solution to the mentioned concentration in DMEM for experiments. The morphology of PM at 4 h with and without (control) treatment with R848, IFN-α and R848 with IFN-α was observed under microscope for each independent experiment. There were no noticeable changes in morphology for these stimuli, and a representative image has been shown in Supplementary Fig. 1b. R848 were purchased from Invivogen Inc., San Diego, CA.; and pure (> 95%), and recombinant mouse IFN-α (Cat#12100-1) with activity > 1X105 U/ml were purchased from PBL Assay Science, 131 Ethel Road West, Suite 6 Piscataway, NJ 08854, USA.

Total RNA isolation and cDNA library preparation for transcriptome sequencing (RNA-seq)

Total RNA was extracted using RNAiso Plus (Takara Bio Inc., Shiga, Japan) and a QIAGEN RNeasy® Mini Kit (QIAGEN, Hilden, Germany). PM cells were completely lysed using RNAiso Plus and then 200 μl of chloroform was added. The tubes were then inverted for 5 min. The mixture was centrifuged at 12,000 × g for 15 min at 4 °C, and the upper phase was collected and transferred to a new tube. Same volume of isopropanol alcohol was added into it and was inverted 5–6 times and was kept on ice fully emerged for 10 min. Then the mixture was passed through an RNeasy mini column. The column was washed with wash buffer. To elute the RNA, RNase-free water (30 μl) was added directly onto the RNase mini column, which was then centrifuged at 12,000 × g for 3 min at 4 °C. To deplete ribosomal RNA (rRNA) from the total RNA preparations, a RiboMinus Eukaryote kit (Life Technologies, Carlsbad, CA) was used according to the manufacturer’s instructions. RNA libraries were prepared using a NEBNext® Ultra™ directional RNA library preparation kit for Illumina® (New England Biolabs, Ipswich, MA). The obtained rRNA-depleted total RNA was fragmented into small pieces using divalent cations at elevated temperatures. First-strand cDNA was synthesized using reverse transcriptase and random primers, and second-strand cDNA synthesis was then performed using DNA polymerase I and RNase H. The cDNA fragments were processed using an end-repair reaction after the addition of a single ‘A’ base, followed by adapter ligation. These products were purified and amplified using PCR to generate the final cDNA library. The cDNA fragments were sequenced using an Illumina HiSeq 2000. Biological triplicate RNA sequencing for each condition was performed on independent RNA samples from either R848, IFN-α or combination stimulated PM: control 4 h (3 samples); R848 4 h (3 samples), IFN-α 4 h (3 samples), and R848 with IFN-α 4 h (3 samples).

Differentially expressed gene analysis using RNA-seq data

FASTQ files from RNA-seq experiments were clipped, trimmed of adapters, and the low-quality reads were removed by the trimming algorithm “Trimmomatic” (Bolger et al. 2014). Quality-controlled FASTQ files were aligned to Mus musculus UCSC mm10 reference genome sequence using the STAR (version 2.5.1) aligner software (Dobin et al. 2013). To measure differential gene expression, DESeq 2 (Love et al. 2014) was used. A subset of condition-specific expression was defined as showing a log2 fold change ≥ 2 and P ≤ 0.01 in expression between controls, R848-, IFN-α-, and R848 with IFN-α-treated samples. The RNA-seq experiments were visualized using HOMER (Heinz et al. 2010) after custom tracks were prepared for the UCSC Genome Browser (http://genome.ucsc.edu/). The RNA-seq data sets were deposited in the NCBI Gene Expression Omnibus database under dataset accession numbers GSE79898 and GSE104056.

Quantitative real-time PCR (qRT-PCR)

The reverse transcription of the RNA samples was performed as per distributor protocol using 2 μg of total RNA, 1 μl of oligo(dT)-primer (per reaction) and a Prime Script 1st strand cDNA Synthesis Kit (Takara Bio Inc., Shiga, Japan). The oligo(dT) primer and RNA templates were mixed and denatured at 65 °C for 5 min and then cooled for 2 min on ice. Prime Script buffer (5×), RTase and RNAse inhibitor were added to the cooled template mixture and incubated for 1 h at 50 °C before an enzyme inactivation step was performed at 70 °C for 15 min. qRT-PCR was performed using SYBR Green PCR Master Mix (Takara Bio Inc., Shiga, Japan) and a 7500 Real-Time PCR System (Applied Biosystems, Waltham, MA). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control. Complementary DNA samples were diluted 1.5-fold, and qRT-PCR was performed using an ABI-7500 Real-Time PCR System (Applied Biosystems, Waltham, MA) with SYBR Premix Ex-Taq II (Takara Bio Inc., Shiga, Japan) according to the manufacturer’s instructions. The reactions were performed in a total volume of 20 μl that contained 0.4 mM of each primer (Table 1). Each PCR series included a no-template negative control that contained water instead of cDNA and reverse transcriptase-negative control for each gene. Triplicate measurements were performed for all reactions. Different samples were evaluated using 96-well plates in the gene expression experiments, and all samples were analyzed on a single plate for the endogenous control experiments. The results were analyzed using the critical threshold (ΔCT) methods in the ABI-7500 software program with the Norm finder and geNorm-plus algorithms. The primers were designed using Primer Express software (Applied Biosystems, Waltham, MA).

Table 1

List of primers used in qRT-PCR studies

Gene	Forward (5′→3′)	Reverse (5′→3′)
GAPDH	AAGGTCGGAGTCAACGGATT	CTCCTGGAAGATGGTGATGG
TNFA	CAGGCGGTGCCTATGTCTC	CGATCACCCCGAAGTTCAGTAG
TNFSF10	CCAACCACCAGGCTACAGG	GCGTCACACTCAAGCTCTG
IL1A	ACTGCACCCAAACCGAAGTC	TGGGGACACCTTTTAGCATCTT
IL1B	GAAATGCCACCTTTTGACAGTG	CTGGATGCTCTCATCAGGACA
CCL2	TGTACCATGACACTCTGCAAC	CAACGATGAATTGGCGTGGAA
CCL3	TTCCTGCTGTTTCTCTTACACCT	CTGTCTGCCTCTTTTGGTCAG
CCL4	CTGGGCCAGATAAGGCTCC	CATGGGGCACTGGATATTGTT
CXCL11	ATTTCCACACTTCTATGCCTCCT	ATCCAGTATGGTCCTGAAGATCA
CXCL10	TGCTGGGTCTGAGTGGGACT	CCCTATGGCCCTCATTCTCAC
TLR7	TTGCATCTGGCGTCTACACT	GGTTTAGGAGGGCAAGGGTG
IFNA	GGCACAGAAGTGTTCCATAAAGT	GAGGCAGGGCTTCCGATAG
IFNB	AGCTCCAAGAAAGGACGAACA	GCCCTGTAGGTGAGGTTGAT
IFNAR1	ACATCGACCCGTCCACAGTAT	CAGAGGGGTAGGCTTGTCTC
IFNAR2	TGGGTCTGCCACAAATGGAG	TCCAGTGTTTGCGTGTTACTC
ISG20	GAGGGCTGTTGGTTCTTGACT	CCTCGGGTCGGATGTACTTG
IFIT1	CCCTGACGACGTGGACTATG	GCCGACAGAGTGATCTTGGT
IFIT2	CAGATGGTCAATTGGTGCCA	TGCAAGAACCCCTGGATCTC
IFIT3	AGCAGGAGGTCTCTGACAATG	GGCTTCCTCTAAACTGTTGAGC
IL6	TGCCCGAACAAGGCTCTTC	CAGCCAGTTGATGCTCTGC

Functional annotation

To functionally annotate the most significant genes, gene ontology (GO) analysis was performed by DAVID (Database for Annotation, Visualization and Integrated Discovery), version 6.8 (da Huang et al. 2009). GO was analyzed using a modified Fisher’s exact P value in the DAVID program. P values less than 0.001 were considered greatly enriched in the annotation category.

Canonical pathway and upstream regulator analysis of datasets

An ingenuity pathway analysis (IPA) (Ingenuity Systems, http://www.ingenuity.com, CA) was performed to analyze the most significant canonical pathways and upstream regulator analysis in the datasets as previously described (Kramer et al. 2014). The genes from datasets associated with canonical pathways in the Ingenuity Pathways Knowledge Base (IPAKB) were considered for literary analysis. The significance of the associations between datasets and canonical pathways was measured in the following manner: (1) the ratio of the number of genes from the dataset that mapped to a canonical pathway was divided by the total number of genes that mapped to the same canonical pathway and (2) Fisher’s exact test for a P value indicating the probability that the association could be explained by chance. After uploading the datasets, upstream regulator analysis was used to predict the upstream transcriptional regulators on the literature and compiled in the IPAKB. Gene networks were algorithmically generated based on connectivity. The analysis examines how many known targets of the upstream regulators are present in the dataset and also the direction of change. The graphical representation of molecular relationships between upstream regulator and gene products was based on the biological relationship between two nodes was represented as an edge (line). All edges were supported by at least one reference from the literature, text book or canonical information in IPAKB. The intensity of node color represented the degree of up-regulation (red). The nodes were displayed using shapes to represent functional classes of gene products.

Transcription factor-binding motif enrichment analysis

NCBI reference sequence mRNA accession numbers were subjected to transcription factor-binding motif analysis using the web-based software Pscan (Zambelli et al. 2009). The JASPAR (Portales-Casamar et al. 2010) database of transcription factor-binding sequences was analyzed using enriched groups of − 950 base pair (bp) sequences to + 50 bp of the 5′ upstream promoters. The range − 950 to + 50 was selected from the range options in Pscan to obtain the best coverage for a − 1000 to + 50 bp range.

Statistical analysis

The data were analyzed using Origin Pro 8 software (Origin Lab Corporation, Northampton, MA). Each value is expressed as the mean ± standard error of the mean (SEM). All qRT-PCR data were analyzed with SPSS 17.0 software (SPSS Inc., Chicago, IL). The data were tested using one-way ANOVA followed by Tukey’s HSD post hoc test. *P < 0.01 and **P < 0.001 were considered significant.

Results

Time point determination and concentration optimization of R848 and IFN-α required for the production of selective inflammatory cytokines and ISGs

First, to determine the time point for optimum acute stimulation, we incubated highly purified PM with ultrapure R848 (1 µM), IFN-α (100 U/ml) and a combination of the two for up to 24 h. Then we measured the mRNA expression of a selective inflammatory cytokine TNF-α and an ISG IFIT2 by qRT-PCR. These two genes often show upregulation in inflammatory conditions associated with TLRs and type I interferon receptor activation. Our qRT-PCR results showed an undetectable expression of these two genes in resting cells, while they were expressed optimally and synergistically at 4 h and declined onward (Fig. 1a). We again attempted to optimize the concentration by stimulating cells with varying concentrations of R848 and IFN-α alone, as shown in Fig. 1b. As with increasing duration of exposure, we found increased expression of TNF-α with increasing concentration of R848 and increased expression of IFIT2 with increasing concentration of IFN-α (Fig. 1b). Notably, both TNF-α and IFIT2 showed high response to R848 and IFN-α, respectively, with time and concentration (Fig. 1a, b). Based on the abovementioned results, it seems likely that in comparison to R848 and IFN-α alone, R848 with IFN-α is more efficient to activate genes comprising cytokines and ISGs in PM in a time- and concentration-dependent manner.

Fig. 1

Effect of R848 and IFN-α on stimulation of selective inflammatory cytokines and ISGs in microglial cells. a Primary microglial (PM) cells were stimulated with R848 (1 μM), IFN-α (100 U) and R848 with IFN-α for up to 24 h. b PM cells were stimulated with R848 (0.1–5 μM) and IFN-α (10–500 U) for 4 h. Quantitative mRNA expression of TNF-α and IFIT2 was measured by qRT-PCR. qRT-PCR data were pooled from three independent experiments, each in triplicate. Data presented as mean ± SEM; *P ≤ 0.01 and **P ≤ 0.001, and determined using one way ANOVA followed by Tukey’s multiple comparison test

Global gene transcription dynamics involved in R848-, IFN-α- or co-treated cells

To gain more insight into the genome-wide kinetics of gene transcription, we again stimulated the PM either with R848 (1 µM) or IFN-α (100 U/ml) or with a combination of the two for 4 h in three biological replicates. We chose this 4 h time point based on a preliminary experiment to determine the optimum stimulation time for immune gene expression in PM. After quality check, isolated mRNA from induced PM were subjected to RNA-seq assay as per previously described protocol (Pulido-Salgado et al. 2018). Approximately 3 million 100-nucleotide reads obtained from each sample were used to map to the mouse genome and to quantitate mRNA expression. Next, principal component analysis (PCA) examined the congruency among biological replicates. The PCA showed a good separation and a high level of consistency between biological replicates of the same population in PM (Suppl. Fig. 2a).

The analysis of the RNA-seq data distinguished the number of genes that were differentially expressed: 393, 394 and 836 genes were upregulated, whereas 61, 19 and 57 genes were downregulated, respectively, in R848-, IFN-α- and R848 with IFN-α-treated PM (log2 fold change ≥ 2, and P < 0.001) (Fig. 2a). Among the upregulated genes, we identified 332 and 362 genes in the R848 and IFN-α samples, respectively, that showed overlap with the set of genes upregulated in the co-treated samples, while only 62 and 32 genes, respectively, were upregulated uniquely by R848 and IFN-α (Fig. 2b). Importantly, we found 295 R848 with IFN-α-upregulated genes that were induced synergistically (log2 fold change ≥ 2 over R848 or IFN-α alone induction) (Fig. 2c). Of note, 25 R848- upregulated genes were repressed by IFN-α, while no IFN-α-upregulated gene was repressed by R848 (log2 fold change ≥ 2) (Fig. 2c). A large number of cytokine, chemokine and interferon response genes important for innate responses were non-synergistically induced upon co-activation of PM with R848 and IFN-α. Nonetheless, synergistically induced genes that showed robust expression included a substantial number of cytokines (TNF-α, IL1A, IL1B, IL12B, IL6, etc.), chemokines (CCL2, CCL3, CCL4, CCL7, CCL9, CCL12, CXCL10, etc.), and interferon and ISGs (MX1, MX2, IFIT1, IFIT2, IFIT3, IFITM2, IFITM3, OASL1, OASL2, etc.) (Fig. 2d). These results suggest that there are gene expression signatures with a combination of cytokines, chemokines and ISGs that represent optimum PM activation.

Fig. 2

Effect of R848, IFN-α and R848 with IFN-α on global gene transcription in PM cells. a Pie charts displaying the number of genes altered by R848, IFN-α and R848 with IFN-α in PM (log2 fold change ≥ 2, P < 0.0001). b Venn diagram depicting the overlap of the number of upregulated genes of R848- and IFN-α-treated PM with R848- with IFN-α-treated PM. c Upper pie chart displaying the genes synergistically/non-synergistically upregulated in R848 with IFN-α-treated PM. d Heat map representation depicting cytokine, chemokine and interferon response genes induced in the global RNA-seq experiments (P ≤ 0.01, and log2 fold change ≥ 2 over R848 and IFN-α upregulation). e Heat map representation depicting R848-upregulated and IFN-α-suppressed genes identified in the global RNA-seq experiments (P ≤ 0.01, and log2 fold change ≥ 2). Heat maps were generated with the Multi Experiment Viewer (version 4.8) software. Data represent three biological replicates of single isolation

Among the R848-upregulated/IFN-α-suppressed genes, selective members of the CC and CXC chemokine families (CXCL1, CXCL2, CXCL3, CXCL5 and CCL22), interleukin-1 receptor kinase (IRAK3), and immunoglobulin superfamily (IGSF6) were noticeably repressed by IFN-α (Fig. 2e). Remarkably, from the overall expression data, we observed that R848 treatment induced significantly higher cytokine and chemokine gene expression levels, and IFN-α treatment induced the expression of interferons and ISGs, while R848 with IFN-α induced the expression of cytokines, chemokines and ISGs in activated PM.

To examine the downregulated genes in all three conditions, we found that some genes previously known to be involved with immune regulation or to be essential for PM homeostasis, such as chemokine receptor (CXCR4), colony-stimulating factor 1 receptor (CSF1R), G-protein coupled receptor (GPR157), prostaglandin synthase (PTGS1) and retinoid X receptor α (RXRA), were downregulated significantly in co-treated cells (Suppl. Fig. 2b, c). Among them, CSF1R, the cell surface receptor for CSF1, highly expresses during amyotrophic lateral sclerosis, and inhibition of CSF1R slows the progression of ALS in mouse models (Martinez-Muriana et al. 2016). RXRα, a key nuclear receptor, attenuates host antiviral responses by suppressing IFN genes (Ma et al. 2014). These data suggest that alterations in the expression levels of these transcripts during pathologic conditions not only reflect unique functional capabilities but also can be used as potential targets to identify these cells in distinct physiologic conditions.

Divergence of the expression of TFs and paradox of promoter conservation in R848, IFN-α or R848 with IFN-α-stimulated microglial cells

We next considered that this synergistic responsiveness of PM might be controlled by TFs, as they initiate and regulate the transcription of genes. To investigate this possibility, we found a group of approximately 40 TFs that showed differential expression (log2 FC ≥ 2) in co-treated PM. Among these TFs, we found that several members of the activator of transcription factor (ATF), nuclear factor kappa B (NF-κB), interferon regulatory factor (IRF), and signal transduction and transcription activation (STAT) families were significantly upregulated in activated PM (Fig. 3a). Noticeably, higher expression levels of NF-κB family members were observed in R848-treated and co-treated cells, while higher expression levels of IRF and STAT family members were observed in IFN-α-treated and co-treated cells, and a similar level of expression of ATF family members was observed in all three conditions (Fig. 3a). However, the expression levels of a large group of TFs, such as IRF3, Kruppel-like factors [KLFs (3, 7, 9, 10, 13 and 16)], ATFs (1, 2, 6 and 6B), CREB3, and STATs (5B and 6), were marginally or mostly unaffected in all three conditions (Fig. 3a). In particular, ATFs (3 and 4), AT-rich interactive domain-containing protein 5A (ARID5A), basic leucine zipper ATF-like TF (BATF), CCAAT/enhancer-binding protein (CEBPB), cAMP-responsive element-binding proteins (CREB5 and CREM), forkhead box protein C1 (FOXC1), KLF6, suppressor of cytokine signaling (SOCS7), nuclear TF X-box-binding-like 1 (NFXL1), JUNB, and selective members of the NF-κB (RELA), and IRF (1, 2, 7, 8 and 9) families showed robust synergistic expression in co-treated PM (Fig. 3b). However, we found that STAT6, KLF11 and three members of the IRFs (2, 4 and 6) were downregulated not only in R848-treated and INF-α-treated cells but also in co-treated cells (Fig. 3b).

Fig. 3

Effect of R848, IFN-α or R848 with IFN-α on key TFs in activated PM cells. Heat map representation depicting a the expression levels of the TFs that were dysregulated by R848, IFN-α and R848 with IFN-α (P ≤ 0.0001, and log2 fold change ≥ 2). Data represent three biological replicates of single isolation. b Genes commonly upregulated or suppressed by R848, IFN-α and R848 with IFN-α (P ≤ 0.0001, and log2 fold change ≥ 2). Heat maps were generated with the Multi Experiment Viewer (version 4.8) software. c Patterns of TFs motif enrichment within the promoters of the commonly induced genes (P ≤ 0.01, and log2 fold change ≥ 2). d Venn diagrams displaying the commonly induced genes regulated by IRF1, IRF3/7, NF-κB1 and STAT1 in R848 with IFN-α-treated cells. e The IRF1- and IRF7-connected immune response genes as defined by the IPA molecule activity predictor of commonly induced genes (P ≤ 0.01, and log2 fold change ≥ 2). The graphs represent the mean fold values of enrichment relative to IgG/control from three independent experiments, each in triplicate. Data are mean ± SEM *P ≤ 0.01 and **P ≤ 0.001 compared with control

To evaluate the functional role of TFs for the transcription of expressed genes during co-stimulation, we conducted motif analysis to identify TF-binding motifs within − 950 to +50 bp windows relative to genomic loci TSS using a computational approach with Pscan software (Zambelli et al. 2009, Portales-Casamar et al. 2010). Importantly, our analysis results showed that the most significant motif enrichment was in IRF1 (P value: 1.03E−51) and IRF7 (P value: 2.50E−63) rather than in NF-κB1 (P value: 2.87E−16) and STAT1 (P value: 1.97E−17) within the synergistic gene set (Fig. 3c). Next, we determined the number of gene promoters having IRF1- and IRF7-binding motifs, as such results predicted that 71.25% and 77.96% of co-stimulation-induced upregulated genes met the promoter occurrence distribution score (Fig. 3d; Table 2). In contrast, only 49.15% and 53.89% of co-stimulation-induced upregulated genes met the promoter occurrence distributions score for NF-κB1 and STAT1, respectively (Fig. 3d). These results were consistent with our transcriptomic data for IRF1, 7, NF-κB1 and STAT1 (Fig. 3e). In addition to DNA-binding factor analysis, we also applied IPA software (Kramer et al. 2014) to identify target genes that were directly or indirectly activated by the identified TFs. The assessment of upstream regulators by IPA (Kramer et al. 2014) similarly revealed that the expression levels of most synergistically upregulated genes were also directly regulated by the identified TFs, including IRF1 (activation z score: 5.066) and IRF7 (activation z score: 7.236) (Fig. 3e; Table 3).

Table 2

Top 50 IRF1- and IRF7-binding motif sequences of co-induced genes in PM cells (log2 FC ≥ 2, and P ≤ 0.001)

Gene symbol	Score	Position	Sequence	Strand
IRF1 occurrence position distribution (score ≥ 0.791)
CD274	0.964021	− 350	TTTCACTTTCACTTTTAGTTT	+
MX2	0.95806	− 944	ACTTAGTTTCACTTTCATTTC	−
PARP10	0.929686	− 4	AGTCAGTTTCACTTTTGTTTT	+
TRIM21	0.923749	8	TTTCACTTTCAGTTTCCTCTC	−
CDK6	0.90797	− 344	TTCTACTTTCAGTTTTTCTAC	−
IFIT1	0.906538	− 110	CTTCAGTTTCACTTTCCAGTC	+
IFIT1BL2	0.906538	− 56	CTTCAGTTTCACTTTCCAGTC	−
PARP14	0.904022	− 59	AACTTCTTTCGCTTTCATTTC	−
IFIT3	0.901236	− 162	GGTAAGTTTCACTTTCCTCTT	+
IFIT3B	0.901236	− 194	GGTAAGTTTCACTTTCCTCTT	+
TRIM26	0.900241	− 166	TTCCGATTTCACTTTCCTTTT	+
IL27	0.899216	− 84	GCCCAGTTTCACTTTCTGTCC	−
IL15	0.897907	− 316	GGGCTCTTTCTCTTTCACTTT	+
OAS3	0.895126	− 80	CTTCACTTTCGTTTTCTCCTC	−
MX1	0.882159	− 918	CTCTGGTTTCCGTTTCATTTC	−
CCRL2	0.870084	− 656	TTATAGTTACACTTTCCGTTT	+
CXCL11	0.856754	− 390	CTTTACTTTTTTTTTTCCTTC	−
PARP9	0.853857	− 252	GTTTGGTTTTGGTTTTGGTTT	−
CCL8	0.853573	− 85	TCTTGCTTTCATTTCCCATTA	+
OASL2	0.850481	− 250	GTTTGGTTTTGTTTTTGTTTT	−
CASP12	0.849561	− 798	TTTTATTTTTATTTTTTATTT	+
CCL5	0.847906	− 156	TTTCAGTTTTCTTTTCCATTT	+
IL6	0.844323	− 292	TGTGAATTTCAGTTTTCTTTC	+
OASL1	0.840788	− 655	TATGAGTTTCTCTTTTCCTCG	+
TRIM13	0.837021	− 374	GTTTTGTTTTGTTTTTTGTTT	−
CXCL10	0.836646	− 218	GGTAAGTTTCACTTTCCAAAG	−
CCL2	0.836297	− 117	TTCAACTTCCACTTTCCATCA	+
CASP1	0.832613	− 561	GTTGGCTTTTTTTTTTTTTTT	+
CD47	0.83234	− 575	TATCTGTTTTTCTTTCTTTGT	−
CD40	0.831735	− 505	ACCCAGTTTCTCTTTCTTGAG	−
IFIT2	0.831186	− 543	GTTGAGTCTCAATTTCAATTT	+
CD180	0.828598	− 504	GATTACTTTCTCTCTCACCCT	−
CCL12	0.82757	− 305	TGTGACTTCTAGTTTCCTTTC	+
MMP13	0.817673	− 158	AGATGCCTTCATTTTCCATTT	+
IL1A	0.81631	− 330	TCCTTGTTTGGCTTTCACTCT	+
PARP8	0.812007	− 55	ACGGAGTCTCACTTTCTCCCT	−
IFITM3	0.81144	14	CGGCAGTTTCGGTTTCTCAGA	−
IL11	0.806313	− 847	TATTCCTTTTTCTTTTGGTCC	+
CCL7	0.806135	− 241	TTTTTTTTTTTTTTTTTTTTT	−
OAS2	0.805292	− 611	GTGAGGTTTCTTTTTGTGTTT	+
MMP10	0.804619	− 518	ACCTACTCTCTGTTTCAGAAT	−
MMP8	0.804484	− 204	ACAGTCTTCCAGTTTCTGTCT	−
CASP4	0.802083	− 4	AGTAACTTTCATTTTACTCTG	+
MMP28	0.801113	− 210	GCTTGGTTCCAGTTTCCCAAA	−
IL12B	0.797428	− 72	TTCTACTTTGGGTTTCCATCA	+
TRIM34B	0.797188	− 334	ACCAAGTCTCACTTTTCGTCC	−
CCL4	0.796758	− 810	CCTTACTTTGAGTTTGACTGT	+
TNF	0.796345	− 155	CCTCTGTCTCGGTTTCTTCTC	−
CASP7	0.795638	− 492	TGTCTCGTTCTGTTTTTGTTT	−
CX3CL1	0.793004	− 464	TCTGGGTCTCAGTTTCCCCAC	+
IRF7 occurrence position distribution (score ≥ 0.799)
IL27	0.93915	− 81	CAGAAAGTGAAACT	+
IFIT3	0.930494	− 158	AGGAAAGTGAAACT	−
IFIT3B	0.930494	− 190	AGGAAAGTGAAACT	−
CCRL2	0.927443	− 771	CAGAAAATGAAACT	−
CXCL10	0.917424	− 215	TGGAAAGTGAAACT	+
IFIT1	0.917424	− 106	TGGAAAGTGAAACT	−
IFIT1BL2	0.917424	− 53	TGGAAAGTGAAACT	+
TRIM21	0.912534	11	AGGAAACTGAAAGT	+
CD274	0.903844	− 340	ACGAAACTAAAAGT	−
GBP7	0.903346	− 55	CTGAAACTGAAACT	−
ISG15	0.902317	− 68	CCGAAACAGAAAAT	+
DUSP28	0.896727	− 894	ATGAAAGTGAAACC	+
IFITM3	0.896219	17	GAGAAACCGAAACT	+
CD40	0.879632	− 502	AAGAAAGAGAAACT	+
CASP4	0.878122	0	AGTAAAATGAAAGT	−
TRIM26	0.876245	− 162	AGGAAAGTGAAATC	−
NOS2	0.875728	− 885	ATGAAAGTGAAATA	−
CASP12	0.870545	− 2	TCAAAACCGAAAGC	−
IRG1	0.870451	− 67	ACAAAAGTGAAAGG	+
CXCL11	0.862913	− 123	ACAAAAGAGAAACT	+
TRIM14	0.855474	− 4	CAGAAATCGAAACC	−
CCL5	0.853544	− 132	CATAAAATGAAAAC	−
TLR2	0.852267	− 773	GAGAAAGAGAAAAT	+
IFI44	0.84368	− 8	CGAAAACTGAAACT	−
IL15	0.842898	− 306	AGAAAAGTGAAAGA	−
IFIT2	0.83864	− 640	GGGAAAGTAAAAAT	−
CCL2	0.836435	− 113	TGGAAAGTGGAAGT	−
DUSP2	0.830229	− 430	TCGATAGCAAAAAT	−
CXCL16	0.828173	− 310	CCTAAAGTGAGATT	+
ISG20	0.827865	− 131	TCCAAAATGACAGT	−
MMP8	0.82728	− 780	ACGAAAACTAACAT	−
CD38	0.826755	− 25	AAGCAAGTGAAAAA	+
CXCL1	0.826514	− 164	CAAAAAGCAAAAAT	+
CCL4	0.825604	− 929	CAGAAACAGAAAAC	−
IL6	0.824672	− 288	AGAAAACTGAAATT	−
TLR1	0.824026	− 880	ATCAAAGTGAAATC	+
MMP3	0.8225	− 510	ACAAAAATAAAAGA	+
IL12B	0.82091	− 68	TGGAAACCCAAAGT	−
CD69	0.817776	− 637	AGGAAACAGAAAGC	−
TRIM25	0.816272	− 12	TCGAAACTGAACAG	−
CCL12	0.812674	− 187	TAGACAGCGAAACA	−
MMP10	0.812058	− 349	TGCAAAGTGAATGT	−
CXCL13	0.811598	− 937	TCCAAATCAAAAGT	+
IFI204	0.811567	− 163	GGGAAATTGAAAGC	+
TNF	0.808402	− 152	AAGAAACCGAGACA	+
IL10	0.807648	− 250	GCTAAAAAGAAAAA	+
TLR3	0.806525	− 536	ACAGAAGTGAAAGC	−
IL1R1	0.806397	− 711	AAAAAACCAAAAAT	+
IL1A	0.803145	− 870	TGGGAACTGAAACT	+
CASP1	0.798994	− 876	CATAAAATGACAGT	−

Table 3

Genes predicted to be regulated by IRF1 and IRF7 as identified by IPA analysis in R848- and IFN-α treated PM cells

R848-induced cells		IFN-α-induced cells
IRF1 predicted to be activated (35 genes) (P = 2.23E−26)	IRF7 predicted to be activated (56 genes) (P = 1.26E−31)	IRF1 predicted to be activated (39 genes) (P = 3.23E−29)	IRF7 predicted to be activated (42 genes) (P = 4.44E−41)
CXCL1, CCL3, CCL2, CXCL5, S100A8, C3, CXCL3, TLR1, CXCL2, TLR2, CCL5, CXCL11, CCL4, CCL7, CXCL10, CASP4, MYD88, CHIL1, IL1B, ZC3H12A, TNIP1, IL1A, GBP5, OLR1, LYN, C4B, IL27, HCK, NLRP3, IFI202B, CCL12, CXCL13, CLEC7A, TNFAIP3, CD14	C3, CASP1, CASP4, CASP8, CCL2, CCL5, CCRL2, CD274, CD40, CD86, CMPK2, CXCL10, CXCL11, CXCL9, EIF2AK2, FAM26F, FCGR1A, FGL2, GBP2, GBP3, GBP5, GBP6, HERC6, ICAM1, IFI16, IFI35, IFI47, IFIT1, IFIT1B, IFIT2, IFITM3, IGTP, IL15, IL15RA, IL1R1, IL6, ISG15, KLRK1, MX1/MX2, OAS1, OASL, PTGS2, RSAD2, SLFN1, SLFN13, SLFN2, SLFN5, SP110, TAP1, TGTP1, TGTP2, TLR9, TNF, TNFSF10, TRAFD1, USP18	CASP1, CCL5, CCL7, CD14, CD40, CEBPB, CEBPD, CLIC4, CXCL10, DTX2, FCGR2, GBP5, GCH1, GFAP, IFI47, IKB, IL12B, IL1B, IRAK2, LCN2, LILRB4, NLRP3, NOS2, MAPK, PLSCR1, PRKAA, PSME2, SERCA, SLC16A10, SLFN2, SPATA13, SUSD6, TANK, TNFA, TNFAIP2, TNFAIP3, TNFSF15, TNIP1, TRIP10	APOBEC3, LY86, TLR2, TLR3, TLR7, TLR9, ISG20, NLRC5, NOD2, TMEM173, CASP4, MYD88, NOD1, OASL2, OASL1, MX1, MX2, LYN, BST2, IFI202B, OAS1B, RIPK2, OAS1A, RNF135, EIF2AK2, IFIH1, C3, CSF1, IFITM3, PML, KLRK1, OAS3, RSAD2, SP110, OAS2, NAIP6, IRGM1, TRIM25, AIM2, DDX58, IFIT3, TRIM56, IFIT2, IFIT1

R848 with IFN-α combination enhances transcription of epigenetic modifiers, G protein-coupled receptors, nuclear receptors and matrix metalloproteinases

In addition to cytokines, chemokines, ISGs and TFs, our RNA-seq transcriptomes revealed synergistic expression of transcripts encoding protein for ligands, receptors, and enzymes, which strongly suggests new microglial functions in neuroinflammatory diseases. Microglia express genes of interest such as several epigenetic modifiers, including histone methyltransferases (SETDB2, MLLT3, MLLT6 and SUZ12), histone demethylases (KDM6B and KDM1A), histone deacetylases (HDAC3 and SIRT1) and histone bromodomains (BRD2 and BRD4). When evaluating the synergistic genes, we unexpectedly identified 35 previously undetected transcripts that encode protein for G protein-coupled receptors (GPR84, GPR18, GPR137, GPR180, GPR107 and GPR37L1), nuclear receptors (NR1D1, NR1D2, NR4A1, RARA and NR3C1), and matrix metalloproteinases (MMP10, MMP12, MMP13, MMP3, MMP8 and MMP9) (Fig. 4a–d). Surprisingly, several of these genes were frequently overexpressed not only in R848-treated but also in IFN-α-treated cells. Previous reports demonstrated that epigenetic modifiers, G protein-coupled receptors, nuclear receptors and matrix metalloproteinases exhibit significant links with transcriptional activation and that result in the synthesis and secretion of inflammatory factors and, in some cases, molecules that suppress immune responses. In addition, several of these modifiers and receptors have also been described as therapeutic targets to modify immune-related diseases (Parks et al. 2004; Huang and Glass 2010, Diehl et al. 2011; Insel et al. 2015; Raghuraman et al. 2016). For example, GPR84 is a pro-inflammatory receptor of microglial cells in a neuropathic pain mouse model (Bouchard et al. 2007), while MMPs (e.g., MMP12, MMP13, MMP3, MMP8) have been described as a therapeutic target for inflammatory and vascular diseases (Yong 2005; Hu et al. 2007).

Fig. 4

Effect of R848, IFN-α or R848 with IFN-α on epigenetic modifiers, G protein-coupled receptors, nuclear receptors and matrix metalloproteinase in PM cells. a Heat map representation depicting the expression of epigenetic modifiers, G protein-coupled receptors, nuclear receptors and matrix metalloproteinases selectively dysregulated by R848, IFN-α and R848 with IFN-α treated (P ≤ 0.01, and log2 fold change ≥ 2) in the global RNA-seq experiments. Heat maps were generated with the Multi Experiment Viewer (version 4.8) software. Data represent three biological replicates of single isolation. b Bar graph displaying commonly upregulated selective epigenetic regulators, G protein-coupled receptors and MMPs (log2 FC ≥ 2). c The transcript abundance of KDM6B, GPR84, MMP13 and NR1D1 genes in control, R848-, IFN-α- and R848 with IFN-α-treated cells. The read count was represented by measuring the average read obtained from triplicate RNA-seq experiments. d UCSC Genome Browser images representing the normalized RNA-seq read density of KDM6B, GPR84, MMP13 and NR1D1 genes in control, R848-, IFN-α- and R848 with IFN-α-treated cells

Functional and pathway analyses of R848 with IFN-α-induced common genes

To gain insight into the biological processes enriched in commonly induced IRF1- and IRF7-targeted genes, we used DAVID Bioinformatics Informatics Resources (da Huang et al. 2009) to classify the results into gene ontology (GO) categories. Based on the molecules present in the dataset, the DAVID GO analysis revealed that the functions most associated with commonly upregulated genes were related to the immune system process and response to stimulus (Fig. 5a). Interestingly, cell signaling and detoxification processes were also significantly enriched in synergistic gene sets (Fig. 5a). To determine the potential canonical pathways of these induced genes, we utilized the IPA (Kramer et al. 2014), a powerful analysis tool that represents the relevant molecular functions based on functional knowledge inputs. The major categories of the canonical pathways were the communication between the innate and adaptive immune cells, the role of pattern recognition receptors in the recognition of bacteria and viruses, activation of IRF by cytosolic pattern recognition receptors, and interferon signaling (Fig. 5b). To corroborate these functional findings, we analyzed the influence of co-action on molecular signaling networks in co-induced PM cells using knowledge-based IPA (Kramer et al. 2014). Commonly induced gene sets revealed signaling networks related to innate and inflammatory responses and to immunological diseases. In particular, in the co-treated group, IRF1 and IRF7 were identified as the central modulator hubs (Fig. 5c). Together, these data imply that the IRF1 and IRF7 pathways may induce the gene expression related to the inflammatory response in PM during R848 with IFN-α exposure.

Fig. 5

Functional annotation and canonical pathways associated with co-induced genes. a GO term enrichment analysis for the “biological process” category of the IRF1- and IRF7-targeted commonly induced genes in the PM cells. The top GO terms are ranked by the gene ontology enrichment. b The most highly represented canonical pathways of the IRF1- and IRF7-targeted commonly induced genes in the PM cells. Pathways ranked by Bonferroni–Hochberg-corrected − log(PB–H) calculated by Fisher’s exact test with the threshold set to 0.05. The line graph shows the ratio of commonly induced genes enriched in each canonical pathway relative to the deposited GO terms in IPA. c Gene networks top 1 displaying interactions among the commonly induced genes at different cellular levels as determined by IPA gene network analysis. The activity of genes highly connected to this network, namely IRF1 and IRF7 in top-1 network function hubs, as assessed using the IPA molecule activity predictor

Identification of unique transcripts of R848 or IFN-α-inducible genes

Both TLR7/8 and IFNAR1 are expressed simultaneously in numerous infections, but individual activation also manifests distinct biological responses in immune cells (Kreutz et al. 2015). To identify TLR7 or IFNAR1 ligation-specific global gene expression in PM, we further examined our RNA-seq data. We found two subsets of 215 and 240 genes specifically expressed (log2 fold change ≥ 2) in cells treated with R848 and IFN-α, respectively. Of note, these genes were unexpressed by cells treated with the alternate agent. Genes that had their expression significantly enhanced by R848 treatment include cytokines (TNF-ΑIP2, TNFRSF1B, TNFSF15, TNFSF9, IL12B, IL1A, IL1B, IRAK1BP1, IRAK3, etc.), chemokines (CX3CL1, CXCL1, CXCL2, CXCL3, CXCL5 and CXCR3), and some TFs (SOCS3, FOXP4, NF-κB2, BATF, CEBPD, CEBPB, REL, RELB, SPI1, etc.) (Fig. 6a). Similarly, genes that had a log2 fold change ≥ 2 in expression with IFN-α treatment include interferons and ISGs (IFI203, IFI27L2A, IFI44, IFI44L, IFIH1, OAS1A, OAS1B, OAS1C, OAS1G, OAS2, OAS3, etc.) and TFs (ARID5A, IRF2, IRF5, NFXL1, POU3F1, STAT1, STAT2, ZNFX1, etc.) (Fig. 6b). We again performed GO analysis of R848- and IFN-α-specific gene sets using DAVID (da Huang et al. 2009). Interestingly, our results suggest that both gene sets were significantly enriched in immune responses and immune system processes (Fig. 6c, d). These results suggest that the functions of either R848- or IFN-α-induced genes were related to the inflammatory response in PM.

Fig. 6

Effects of R848 and IFN-α alone on the induction of genes in PM cells. Heat map representation depicting a R848 upregulated genes unexpressed in IFN-α-treated PM and b IFN-α upregulated genes unexpressed in R848-treated PM as identified in the RNA-seq experiments (P ≤ 0.01 and log2 fold change ≥ 2). Heat maps were generated with the Multi Experiment Viewer (version 4.8) software. c, d GO term enrichment analysis for the “biological process” category of the R848- and IFN-α-induced unique genes. The top GO terms are ranked by the gene ontology enrichment. Data represent three biological replicates of single isolation

Confirmation of R848- and IFN-α-dysregulated genes as identified by RNA-seq

RNA-seq analysis provides highly accurate expression results; however, to gain additional ratification, we attempted to validate a selective set of differentially expressed genes. Most were selected for validation according to their distinct induction upon R848, IFN-α and R848 with IFN-α treatment of PM. Here, we again incubated highly purified PM with these stimuli for 4 h and measured the mRNA abundance by generating single-stranded template cDNA from the mRNA. The cDNA template was then amplified in the quantitative step, during which the fluorescence emitted by labeled hybridization probes or intercalating dyes changed as the DNA amplification process progressed. With a carefully constructed standard curve, qPCR produced an absolute measurement of the number of copies of original mRNA. There is very good agreement between the RNA-seq and qRT-PCR results in terms of the direction of change as well as its magnitude. Among the 13 genes selected for verification, 11 genes (IL1A, IL6, TNF-α, CCL2, CCL3, CCL4, CXCL11, IFIT1, IFIT2, IFIT3 and ISG20) were induced synergistically in co-treated PM (Fig. 7). However, one was insignificant (TNFSF10) in the qRT-PCR analysis compared with the RNA-seq experiments, and another was undetected. Together, these results support the robustness of our RNA-seq data.

Fig. 7

Confirmation of co-induced selected genes. The bar graphs are representing the TNF-α, IL1A, IL6, TNFSF10, IFIT1, IFIT2, IFIT3, ISG20, CCL2, CCL3, CCL4 and CXCL11 gene expression in R848-, IFN-α-, R848 with IFN-α treated cells over the control. The gene expression levels were normalized to the GAPDH transcript levels and compared with the control. qRT-PCR data are pooled from three independent experiments, each in triplicate. Data are mean ± SEM; *P ≤ 0.01 and **P ≤ 0.001 compared to control

Discussion

In the present study, using high-resolution transcriptome analysis, we established the transcriptional profile mostly involved in PM activation in response to R848, IFN-α and R848 with IFN-α stimuli. A few studies have attempted to use qRT-PCR to determine gene expression changes associated with TLR7/8 agonists and IFN-α treatment in macrophages, B cells, and brains (Siren et al. 2005; Pirhonen et al. 2007; Severa et al. 2007; Butchi et al. 2008; Butchi et al. 2011; Poovassery and Bishop 2012; Kreutz et al. 2015; Zimmermann et al. 2016). However, thus far, a genome-wide search for the similarities and differences between the effects of these two treatments on microglial gene expression has not been conducted. Given the ever-increasing importance of microglia to the field of neuroinflammation research, the ability to isolate high yield of primary microglia is the preliminary to be measured to investigate the role of microglial modulation of inflammation. Here, we used 92% pure microglia obtained from 3-day-old ICR mice brain. In support of our study, 90–95% pure were also used to perform highly accurate experiments to interrogate microglial functions in vitro, including cellular phenotyping, transcriptome analysis, cytokines/chemokines release and neuroinflammatory disease modeling. The strength of our analysis, which aimed to provide comprehensive and comparative transcriptional profiles of responses to inflammatory stimulation, was enhanced by the use of RNA-seq to analyze IFN-α-mediated TLR7/8 cross-regulation in murine microglia. Importantly, we evaluated DNA-binding factors that may drive distinct gene expression in R848-, IFN-α- or R848 with IFN-α-primed primary microglia. Our data support the contention that signaling crosstalk occurs between R848 and IFN-α to cross-regulate transcriptional responses that are critical components of the innate immune system and may lead to neuroinflammatory processes.

IFNs could have either immunostimulatory or immunosuppressive functions in inflammation/antiviral responses (Ivashkiv and Donlin 2014). In our study, we found that a set of cytokines/chemokines, antiviral genes, and IRGs associated with inflammation (Holtman et al. 2015, Srinivasan et al. 2017) was synergistically upregulated in response to R848 with IFN-α compared to their upregulation by single treatment with R848 or IFN-α in PM (Fig. 2d). Both the number of cytokines/chemokines, antiviral genes and IRGs and the extent of the fold changes in the synergistically altered genes were significantly higher in R848 with IFN-α compared to R848- or IFN-α-treated PM. Synergistic inflammatory gene expression was also observed upon stimulation with inflammatory stimuli such as TLR agonist combination alone or with IFNs or combination of other cytokines in dendritic cells, macrophages and brain cells (Napolitani et al. 2005; Qiao et al. 2013; Suet et al. 2013; Kreutz et al. 2015; Goldstein et al. 2017). In contrast, IFN-α-suppressed R848 induced significantly higher level of several specific gene families involved in immune responses, including CXCL1, CXCL2, CXCL3, CXCL5, and CXCL9, among others, compared with only IFN-α- or R848-treated PM (Fig. 2e). These data suggest that alterations in the expression levels of these proinflammatory transcripts during pathologic conditions not only reflect unique functional capabilities but also can be used as potential marks to identify these cells in distinct physiologic conditions.

Our RNA-seq data showed that transcripts of selective inflammatory cytokine/chemokine genes were expressed synergistically in PM treated with R848 and IFN-α (Fig. 2d). Inflammatory cytokines/chemokines from microglia/macrophages mediate defense of the host from various pathogens such as viruses (Murray and Wynn 2011; Klein and Hunter 2017). Although the initial immune response to pathogens is achieved by only a limited number of inflammatory cytokines/chemokines, the anti-pathogen effector programs triggered by cytokine/chemokine system are based on the concerted action of hundreds of ISGs (Schoggins et al. 2011). There is a report which demonstrated that ISG IFIT2 deficiency results in uncontrolled neurotropic coronavirus replication and enhanced encephalitis via impaired alpha/beta interferon induction in mouse brain (Butchi et al. 2014). IFIT2 deficiency also does not regulate mRNA expression of inflammatory cytokines and chemokines such as IL1B, IL6, TNF, CCL2, CCL5, CXCL9, CXCL10, and IFN-γ, or even many ISGs both in microglia and macrophages (Butchi et al. 2014), suggesting the importance of IFIT2 in limiting virus infection in the CNS. Only TLR7/8 ligation in microglia and subsequent ISG transcription were modest and context dependent (Pirhonen et al. 2007; Butchi et al. 2008; Butchi et al. 2011; Schoggins et al. 2011; Butchi et al. 2014). The result of concurrent signaling of TLR7/8 and IFNs induced ISGs of IFITs, ISGs, IFITMs is likely to be filler to molecular signatures of efficient inflammatory/defense responses (Fig. 2d).

Gene transcriptional events during cellular activation are largely controlled by designated TFs. Using this array, we also found that a set of TFs was largely affected either synergistically or co-repressed in PM, suggesting that this set of TFs might include important regulators of IFN-α-associated immunostimulatory and immunosuppressive effects (Fig. 3). More importantly, we identified several TFs, including IRFs (1, 7, 8 and 9) and RELA, that were synergistically upregulated in PM. Our results are consistent with those described in a previously published report showing that IRF1 and IRF8 are critical for microglia activation (Masuda et al. 2015). In our RNA-seq analysis, we also identified several other TFs (ARID5A, ATF3, ATF4, CEBPB, CREB3, CREB5, CREM, FOXC1, JUNB, KLF6, SOCS7, SOX9, NFXL1, and STAT4, among others) in synergistically induced PM. Each of these TFs (IRF1, IRF7 IRF8, and RELA, among others) is predicted to be central to some aspect of the synergistic responses and may represent candidates for experimental validation using knockout or overexpression models. We next found that the promoters of synergistically expressed genes were enriched for IRF1 and IRF7 but not for NF-κB1 and STAT1, as shown in Fig. 3.

Another interesting finding is that our RNA-seq analysis identified several important epigenetic regulators that were synergistically induced by IFN-α and R848 in PM. Recently, we found that the histone demethylases KDM1A and KDM4A, the histone methyltransferases NSD3 and SETDB2 and the DNA methyltransferase DNMT3L were strikingly differentially expressed in LPS-induced PM (data not shown). Importantly, our RNA-seq data revealed not only that those epigenetic regulators were strikingly synergistically expressed in IFN-α- and R848-induced PM but also that bromodomain and extra-terminal motif (BET) proteins BRD2 and BRD4, histone demethylase KDM6B, and the histone deacetylase SIRT1 were strikingly synergistically expressed in IFN-α- and R848-induced PM (Fig. 4). Previous studies demonstrated that SIRT1 and SETDB2 can potentially regulate proinflammatory gene expression in macrophages (Chen et al. 2015; Schliehe et al. 2015). However, the mechanism by which those important epigenetic regulators become synergistically activated remains unknown. Determining how these epigenetic regulators, in combination with modified TFs, can regulate inflammatory genes synergistically in microglial cells would be intriguing.

Our genome-wide analysis employing the major experimental uses of microglia, along with the integration of multiple gene sets and bioinformatics analysis, provides the most robust and comprehensive assessment to date of IFN-α-mediated TLR7/8 cross-regulation in murine microglia at the level of the microglial transcriptome. However, changes observed in these studies may reflect not only the gene expression profiles of microglia, but also those of other CNS cells possibly astrocyte or oligodendrocytes.

Conclusions

Conclusively, using global profiling with a bioinformatics approach, we have described herein the molecular signatures induced by R848, IFN-α alone or co-treatment in microglial cells. Our genome-wide, integrative analysis has revealed the integration of signaling crosstalk between TLR7/8 and IFN-α at the level of the transcriptome in association with changes in related TFs. These data may break new ground in the study of the role of microglia in neurological disorders. Our findings provide a better understanding of the complex activation of the IFN-α-induced TLR7/8 cross-regulation occurring in microglia, and this knowledge could be utilized in elucidating novel targets to modulate microglia activation by neuroinflammatory disorders.

Below is the link to the electronic supplementary material.

Supplementary material 1 (TIFF 6863 kb)

Supplementary material 2 (TIFF 3427 kb)