Introduction

The inflammatory component in neurodegenerative diseases

Age is the main risk factor in neurodegenerative disorders, which are the cause of approximately 6.8 million deaths each year (https://www.who.int/news/item/27-02-2007-neurological-disorders-affect-millions-globally-who-report; accessed on September 20, 2022). Alzheimer’s disease (AD) and Parkinson’s disease (PD) are the two most prevalent. Common features of neurodegenerative diseases include deposition of protein aggregates, progressive neuronal death (cortical and hippocampal neurons in AD, and nigrostriatal dopaminergic neurons in PD), and expression of activation markers in microglia surrounding neurodegenerative structures. The mechanism of neuronal death is still under intense investigation; on the one hand, factors such as oxidative stress, mitochondrial dysfunction, and glutamate excitotoxicity are attributed to it. Neuroprotective approaches generally target receptors on neurons. Several preclinical studies have suggested, over the years, that neuroprotection can be afforded via modulation of the functionality of adenosine receptors expressed in neurons [1,2,3]. Four are the receptors for the autocoid adenosine, A1, A2A, A2B, and A3; they all belong to the G-protein-coupled receptor (GPCR) superfamily and are widely distributed in the mammalian body. The A2A and the A2B are coupled to Gs, and their activation leads to increases in intracellular cAMP, whereas the A1 and the A3 are coupled to Gi and their activation leads to decreases in intracellular cAMP. In the last two decades, the adenosinergic system, and particularly the A2AR, has gained relevance mainly due to two factors. A first in class A2AR antagonist, istradefylline, has been approved in Japan and the USA for the adjuvant therapy of patients with PD. This development has raised hopes in Academia and in the Pharmaceutical Industry that there will be more drug approvals targeting adenosine receptors for a wide range of diseases, from obesity to cancer [4,5,6,7,8,9,10].

As suffering neurons may not be able to expand their life on their own, there is the alternative of targeting receptors on glial cells. In diseases that occur with neuroinflammation, microglia acquire interest. In fact, microglial activation, previously considered harmful, is now becoming relevant due to the possibility of directing these cells towards a neuroprotective phenotype [11,12,13,14,15]. The nomenclature: M0 for resting, M1 for pro-inflammatory, and M2 for neuroprotective are used as a first approximation to understand the phenotypic status of a microglial cell in a given environment [16]. Based on the expression of cannabinoid receptors in activated microglia and in primary cultures of a transgenic model of Alzheimer’s disease (APPSw,Ind), it was hypothesized that microglia help to delay cognitive decline by several months. In fact, these animals, which overexpress the Swedish and Indiana mutated form of the human amyloid precursor protein, do not present gross alterations until later in life and among the few differential trends with control animals is a “pre-activated” phenotype in the microglia isolated from 3-day-old pups [17].

A3R agonists have shown promising anti-inflammatory effects in preclinical assays. They have been considered for the therapy of serious diseases coursing with marked immune system dysfunction, such as rheumatoid arthritis and psoriasis [18], among others. Also, A3R agonists have been suggested as having potential to combat retinal neurodegeneration [19], post-ischemic brain damage [20], or neuroinflammation after subarachnoid hemorrhage [21]. Even though there are several pharmacological and functional studies related to the adenosinergic system, the effect of extracellular adenosine on gene expression is poorly understood. Here, we aimed to determine by a transcriptomic approach whether (i) activation of the A3R in activated microglia modulates gene expression, with special attention to genes encoding factors involved in immune responses, and (ii) microglia are biased towards the neuroprotective phenotype.

Methodology

Reagents

Trizol (Ambion Life Technologies, 15696026), chloroform (Sigma_Aldrich, C2432-500ML), isopropanol (PanReac AppliChem, 131090.1211), lipopolysaccharide (Sigma_Aldrich, L4391-1MG), interferon-γ human (Sigma Aldrich; I3265-1MG), 2-chloro-N6-(3-iodobenzyl)-adenosine-5'-N-methyluronamide (2-Cl-IB-MECA; ToCris Bioscience; 1104).

Isolation and activation of microglia

To prepare primary microglia, the brain of C57BL/6J wild-type mice was removed at postnatal days 2 to 4. Microglial cells were isolated and plated at a confluence of 40,000 cells/0.32 cm2. Briefly, the tissue was dissected, and after carefully removing the meninges, the brain samples were digested with 0.25% trypsin for 20 min at 37 °C. Trypsinization was stopped by adding an equal volume of culture medium (Dulbecco’s modified Eagle medium-F-12 nutrient mixture, Invitrogen). Cells were brought to a single cell suspension by repeated pipetting followed by passage through a 100-μm pore mesh. Pelleted (7 min, 200 g) cells were resuspended in supplemented DMEM and seeded at a density of 3.5 × 105 cells/mL in 6-well plates. Twenty four hours later, the medium was replaced by Dulbecco’s modified Eagle medium supplemented with 2 mM L-glutamine, 100 U/mL penicillin/streptomycin, MEM Non-Essential Amino Acid Solution (1/100) and 10% (v/v) heat inactivated Fetal Bovine Serum (FBS) (all supplements were from Invitrogen, Paisley, Scotland, United Kingdom). Cells were maintained in a humid atmosphere of 5% CO2 at 37 °C and were used for RNA extraction after 15 days of culture. Cells were activated for 48 h with 0.01% (v/v) LPS and 0.002% (v/v) IFN-γ in Dulbecco’s modified Eagle medium and, at 24, 32, and 40 h, vehicle or 200 nM 2-Cl-IB-MECA was added.

After this, we proceeded to extraction of total RNA with the Trizol-based method using chloroform and isopropanol to purify the RNA. Quality control was assessed by measuring the 280/260 absorption ratio and the RIN in the services of the University of Barcelona (CCiTUB: Scientific and Technological Centers of the Universitat de Barcelona). Sequencing was done in the German facilities of Novogene using the following protocol: Messenger RNA was purified from total RNA using poly-T oligo-attached magnetic beads. After fragmentation, the first-strand cDNA was synthesized using random hexamer primers, followed by the second-strand cDNA synthesis using either dUTP for directional library or dTTP for non-directional library. For the non-directional library, it was ready after end repair, A-tailing, adapter ligation, size selection, amplification, and purification. For the directional library, it was ready after end repair, A-tailing, adapter ligation, size selection, USER enzyme digestion, amplification, and purification.

The library was checked with Qubit and real-time PCR for quantification and bioanalyzer for size distribution detection. Quantified libraries are pooled and sequenced on Illumina platforms, according to effective library concentration and data amount. Novogene NovaSeq 6000 was the platform, pair-end 150 (PE150) the sequencing strategy, and 151+8+8+151 the sequencing cycles.

The clustering of the index-coded samples was performed according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina platform, and paired-end reads were generated.

For quality control, raw data (raw reads) in “fastq” format were firstly processed through Perl scripts. In this step, clean data (clean reads) were obtained by removing, from raw data, the reads containing adapter, reads containing poly-N, and low-quality reads. At the same time, Quality scores, Q20 and Q30, and GC content were also considered. All the downstream analyses were done using the clean data.

To map the reads, reference genome and gene annotation files were downloaded from the Ensembl genome browser. Hisat2 v2.0.5 was used for both indexing and alignment of clean paired-end reads to the reference genome. We selected Hisat2 as the mapping tool because it can generate a database of splice junctions based on the gene model annotation file. Therefore it provides better mapping results than other non-splice mapping tools.

RNAseq data processing

Feature counts v1.5.0-p3 software was used to count the reads numbers mapped to each gene. And then fragments per kilobase of each gene was calculated based on the length of the gene and reads count mapped to this gene. The expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) considers the effect of sequencing depth and gene length for the reads count at the same time and is currently the most commonly used method for estimating gene expression levels.

Differential expression analysis

Differential expression analysis of two conditions/groups (three biological replicates per condition) was performed using the DESeq2R package (1.20.0). DESeq2 with biological replicates provides statistical means for determining differential expression in gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the “edgeR R” package (version 3.22.5).

Only those genes whose expression comparison had a false discovery rate (FDR) < 0.05, and a fold change (FC) > |1.5| were selected and classified into those whose expression was upregulated or downregulated in each specific treatment.

Gene set enrichment analyses

The STRING online tool (defined as a “database of known and predicted protein-protein interactions” https://string-db.org/ accessed on April 15, 2022) was used for obtaining the GOs of the DEGs whose expression decreased (downregulated) or increased (upregulated) with the following settings: no additional shells (only provided proteins were considered), medium confidence (0.4), and full STRING network, which includes direct (physical) and indirect (functional) interactions. These analyses provide a list of overrepresented GOs in each group of DEGs (those corresponding to the genes with increased expression, and those corresponding to the genes with decreased expression). Then, these GOs were clustered using the REVIGO online tool (http://revigo.irb.hr/), which reduces the number of variables and makes possible to group the GOs by similarity in a 2D plot. Finally, for improving the informative capacity of the resulting graphs, the Cytoscape software (v.3.9.1) (https://cytoscape.org/) was used.

Also, the Enrichr online tool “interactive and collaborative HTML5 gene list enrichment analysis tool” (https://maayanlab.cloud/Enrichr/ accessed on Abril 30, 2022) was used, to assess which transcription factors (TFs) were overrepresented; of the different options offered within Enrichr, we selected TRRUST (Transcriptional Regulatory Relationships Unraveled by Sentence-based Text-mining v.2), which is a manually curated database of transcriptional regulatory pathways and TFs. Afterwards, the STRING online tool was used again for assessing for different interactions between the overrepresented TFs. The STRING settings were the same described in the previous paragraph.

Results

RNAseq in activated microglia in the presence of an agonist of the A3R

The aim of the article was not to assess how adenosine receptors might affect microglial activation but rather how an A3R agonist might affect the phenotype of activated microglia. Accordingly, we designed an experiment in which microglial cells were activated using a standard protocol, i.e., using LPS and IFN-γ and leaving the cells for 48 h before obtaining RNA. Comparison of differentially expressed genes was performed on total RNA isolated from activated primary microglia treated or not with 2-Cl-IB-MECA. To be sure that the A3R agonist would be able to regulate gene transcription, we thought that (i) at least 8 h after addition of the agonist was required for detecting gene activation/repression and that (ii) various doses, spaced over time, were required to mimic the continuous presence of adenosine in the extracellular milieu when microglia are activated in a physiological setting.

RNA isolation was carried out 48 h after microglial activation and 2-Cl-IB-MECA was added to the culture at 24, 32, and 40 h. mRNA sequencing and counts were performed as indicated in methods. Data analysis showed a large number of differentially expressed genes when activated microglia were treated with 2-Cl-IB-MECA. Comparing the two conditions (plus/minus A3R agonist) and with these two criteria: false discovery rate (FDR) < 0.05 and fold change (FC) > |1.5|, the number of genes whose expression was downregulated by agonist treatment (1414) was greater than the number of genes whose expression was upregulated (502). The heat maps are shown in Fig. 1 (see complementary information in supplementary Table S1).

Fig. 1
figure 1

Heatmap of the differentially expressed genes in activated microglia. Columns indicate the sample (quadruplicates for each condition: control in the 4 columns on the right and treated with 2-Cl-IB-MECA in the 4 columns on the left). Only the name of genes related to inflammation and pathways related to cytokines are shown. The heat map shows only those genes that, when comparing control and agonist treatment, met the criteria: FDR < 0.05 and FC > |1.5|. Darker colors indicate downregulated gene expression in microglia treated with 2-Cl-IB-MECA, while lighter colors indicate upregulated gene expression

Full size image

The gene set enrichment analysis (GSEA) was subsequently performed taking into account the genes whose expression was downregulated genes upon A3R agonist treatment. When only downregulated gene expression for transcription factors were considered 3, clusters appeared (Fig. 2A). In fact, using the genes whose expression decreases, several (145) gene ontologies (GOs) were found (Fig. 2). They were mainly clustered within developmental-related and immune-related processes. The rest corresponded to GOs from cell cycle control to regulation of cell adhesion and cell-to-cell communication regulation. Figure 2B shows the detail of the several GOs altered upon A3R agonist activation within immunological-related events. Table 1 shows GOs when genes with increased expression were considered and only those that were significant (FDR < 0.05). The most significant were those concerning regionalization and pattern specification.

Fig. 2
figure 2

Microglial events that are regulated by transcription factors whose gene expression is downregulated by treatment with 2-Cl-IB-MECA, the receptor agonist. A 3 clusters reflecting processes in which genes downregulated upon 2-Cl-IB-MECA treatment are involved. B Immune-related microglial events that are regulated by the genes are downregulated by 2-Cl-IB-MECA treatment

Full size image
Table 1 Gene ontologies overrepresented when considering genes whose expression was increased upon 2-Cl-IB-MECA treatment
Full size table

Analysis of known and predicted interactions, taking into account the protein products of the downregulated genes, reveals cAMP-responsive element-binding protein 1 (Creb1) as a node that connects two networks. One of the networks contains the regulatory factor X and associated proteins, and in the other network, the proto-oncogene ETS 1 (Ets1) and the member of the SMAD family 3 (Smad3) stand out (Fig. 3).

Fig. 3
figure 3

STRING analysis of interactions between the product of genes whose expression is decreased upon A3R agonist treatment. The colors of the connection lines represent the different types of making associations, either from known or predicted interactions: from text mining (lime), experimentally determined (magenta), gene co-occurrence (blue), and co-expression (black). Creb1, cAMP-responsive element-binding protein 1; Ets1, ETS proto-oncogene 1; Smad3, SMAD family member 3; Ep300, E1A-binding protein p300; Srf, serum response factor; Rfx5, regulatory factor X5; Rfxank, regulatory factor X associated ankyrin-containing protein; Rfxap, regulatory factor X associated protein; Stat3, signal transducer and activator of transcription 3; Hmga2, high-mobility group AT-hook 2; Snai1, snail family transcriptional repressor 1; Myb, myeloblastosis family of transcription factors

Full size image

We finally addressed the expression of genes for biomarkers of M1 and M2 microglia. Treatment of activated microglia with 2-Cl-IB-MECA led to both upregulation and downregulation of the expression of pro-inflammatory genes, with the chemokine ligand CXCL11 being the product of the gene whose expression was most increased and the chemokine ligand CCL19 being the product of the gene whose expression decreased in greater magnitude. Following treatment with the A3R agonist, the expression of several genes encoding M2 biomarkers was decreased with the exception of the peroxisome proliferator-activated receptor γ (PPAR-γ) gene, whose expression was increased. Especially relevant was the decrease in the expression of the gene that encodes the chemokine ligand CXCL13. In general, there is no particular trend with respect to microglial polarization, although under the conditions of the A3R activation assay here performed, increases in the expression of genes known to be involved in M2 polarization are not obtained (with the exception of PPAR-γ).

Discussion

The transcriptomics data here presented confirms both A3R expression in activated microglia and a relevant role of the receptor in mediating the regulation of microglial activation by adenosine. One of the relevant findings is the large number of genes whose expression is altered by the activation of A3R. The number of genes in the mouse genome with protein sequence data is reportedly 25,059 [22]; therefore, the number of genes whose expression is significantly regulated by 2-Cl-IB-MECA was about 8% of the total. A limitation of the study is the assay conditions, namely, robust activation using both LPF and IFN- γ, and 3 doses of the agonist at 24, 32, and 40 h after initiation of activation. Our results are in contrast to the reported beneficial effects of 2-Cl-IB-MECA after subarachnoid hemorrhage in aged rats. The anti-inflammatory effects in this study are seemingly mediated by the P38/STAT6 pathway; however, details on the nature of microglial cells are lacking [21]. Differential activation and/or different agonist treatment protocols could lead to quantitative and/or qualitative differences in results. An example is the reduced expression of inducible nitric oxide synthase (iNOS) by 2-Cl-IB-MECA in retinal microglia under elevated hydrostatic pressure mimicking elevated intraocular pressure, that is, glaucomatous conditions [23, 24]. The BV2 microglial cell line treated with LPS expresses the A3R whose activation suppresses TNF-α production by a mechanism that depends on Akt and PI 3-kinase [25]. In a report by Koscsó et al. [26], activation of the A2B receptor in microglia leads to increase in IL10 production, whereas activation of the A3R leads to downregulating the expression of the gene coding for IL10 (Fig. 4). Interestingly, the mechanisms mediated by A2B receptor activation depends on Creb [26], a gene whose expression was found to be regulated upon A3R activation (Fig. 3).

Fig. 4
figure 4

Histogram showing the relative expression of genes for microglial phenotype-related biomarkers. A Degree of variation (log2FoldChange-FC-) of genes for M1 biomarkers. B Degree of variation (log2FoldChange) of genes for M2 biomarkers. Increases in expression due to agonist treatment are in green, and decreases are in red. For all these data, the FDR was < 0.05; that is, genes whose expression was not significantly altered upon A3R treatment are not shown. Only genes relevant to M1/M2 polarization are shown

Full size image

A second relevant finding was that increases in gene expression occurred much less than decreases in gene expression. It is reasonable to hypothesize that adenosine, via A3R, is moderating microglial activation, thus reducing the burden associated to keeping a high degree of gene expression. Transcription factors whose gene expression was regulated by 2-Cl-IB-MECA in activated microglia (Fig. 2, Table 1) are involved in various cellular functions, from the most general to the most specifically related to regulation of immunological responses. Gene ontologies detected after data analysis showed significant clustering with developmental processes, and this fits with the finding that, in microglia, A3R agonists promote ADP-induced process extension and migration [27].

The M0 for resting, M1 for pro-inflammatory, and M2 for homeostasis/restoration nomenclature was first coined for macrophages (see [28]), cells that in the periphery perform a function equivalent to that of microglia in the CNS. It is considered that adenosine receptors in macrophages have relevant regulatory functions, thus being potential targets in the therapy of inflammatory diseases (see [29, 30] for review). As in microglia, A3R activation in macrophages inhibits TNF-α production [31,32,33] and increases the LPS-stimulated IL-10 production [33]. Agonists of the receptor exert anti-inflammatory effects by reducing the production of macrophage inflammatory protein-1α (MIP-1α) in immunostimulated RAW macrophages [34]. In summary, in both microglia and macrophages, the activation of A3R presumably leads to anti-inflammatory actions.

The transcriptomics study presented here does not provide a basis for a M2 polarization unless there are M2 biomarkers that are still unknown. In regard to M1 polarization, the results are conflicting as there are some genes that are upregulated and others that are downregulated. This fact could indicate that there are intermediate phenotypes [35, 36]. The M0/M1/M2 nomenclature has been instrumental in pointing out that microglial activation can progress to neuroprotective phenotypes [16], but the actual phenotypes may not be clearly established based on M1/M2 markers. Transcriptomics studies may help to elucidate differential activation modes, for instance, comparing results coming from cells activated in a more physiological way. To assess the potential of A3R in microglia as a target for neuroprotection, it is also important to play with the doses and find the right time to act, neither too early nor too late. Also relevant is to consider pathways that have not been previously considered. Our results suggest that for better assessing both microglial activation and the potential of A3R agonists in neuroprotection, Creb1, and Smad3 and the RFX transcription factor family must be considered.