Genes and pathways of regulatory T cells regulated by adenosine A2A receptor: A bioinformatics study

Abstract

Objective:

We used a preliminary bioinformatics method to analyze the effects of the adenosine A2A receptor on regulatory T cells.

Methods:

GSE34006 data were retrieved from the GEO database, and the differential expression genes (DEGs) were analyzed. DAVID Bioinformatics Resources and Cytoscape software were used to analyze the GO and KEGG. The STRING database was used to construct the PPI networks. Finally, two kinds of algorithms were used to identify the hub genes of the PPI networks. Targetscan and AnimalTFDB were used to predict the regulatory miRNA and transcription factors (TFs) of the hub genes.

Results:

Two hundred and forty-nine DEGs involved in Treg cells regulated by the adenosine A2A receptor, with 163 genes upregulated and 86 genes downregulated. These DEGs were enriched into 75 GO terms and 5 pathways. Fos, Stat5a, Cebpb and Il10 were the hub genes and were upregulated. 13, 2 miRNAs, 1 and 3 miRNAs were predicted to regulate Fos, Stat5a, Cebpb and Il10, respectively. 1 TF was predicted to regulate Fos, Stat5a and Cebpb.

Conclusion:

This study presents a platform to discover the potential mechanism for adenosine A2A receptor on regulatory T cells. However, the functional studies are further needed to validate the hub genes and pathways.

KEYWORDS:

• Adenosine A2A receptor
• regulatory T cells
• gene microarray
• hub genes
• miRNA
• transcription factors

Introduction

Immune system function networks are regulated by interactions between the variable parts of lymphocytes and their secreted molecules (Andersen 2017; Zhang et al. 2017; Sharabi et al. 2018). The immunosuppressants or immunoregulators, known as suppressor T cells or regulatory T (Treg) cells, inhibit antigen–antibody reactions (Andersen 2017; Sharabi et al. 2018; Zheng and Wang 2020). The word ‘regulatory’ is used to describe the important Treg cell functions of maintaining immunological self-tolerance and immune homeostasis (Andersen 2017; Li et al. 2018).

However, regulatory functions are not only in the immune system; the functions of Treg cells in nonimmunological processes have also been found in recent years in bone repair (Lei et al. 2015; Zhang et al. 2017; Proto et al. 2018; Zheng and Wang 2020). Treg cells accelerate bone healing not just by reducing inflammatory responses in direct or indirect ways (Lei et al. 2015; Proto et al. 2018; Sharabi et al. 2018; Schlundt et al. 2019), but also by participating in cellular crosstalk and recruitment actions within mesenchymal stem cells (MSCs), osteoblasts, osteoclasts, and vascular endothelial cells (Sun et al. 2017; Zhang et al. 2017; Li et al. 2018; Schlundt et al. 2019; Zheng and Wang 2020).

The adenosine A2A receptor, as a G-protein coupled receptor, can form self-association homo-oligomers, like A2A–A2A, or heterodimers, like the dopamine D2 receptor, dopamine D3 receptor, cannabinoid CB1, glutamate mGlu5 receptor, and so on (Guerrero 2018; Jamwal et al. 2019). Thus, the functions of the adenosine A2A receptor include self-functions and polymer-induced functions, such as inflammatory regulation, oxidative stress buffering, ischemia reperfusion protection, and tissue healing (Li et al. 2012; Leone et al. 2015; Mediero et al. 2015; Klein and Bopp 2016; Maj et al. 2017; Guerrero 2018; Lopez et al. 2019; Zheng and Wang 2020).

ATL-1222 or CGS-21680 are substituted n-ethylcarbazole (NECA) derivatives used as classical A2A receptor agonists in many studies (Li et al. 2012; Bao et al. 2016b; Mediero et al. 2018; Zheng and Wang 2020; Zhou et al. 2020). We previously studied the effects of the adenosine A2A receptor agonists on Treg cells in rat tibial fracture models and found that adenosine A2A receptor agonists increased Treg cells in fractured rat blood. The phenomenon was found to be related to accelerated bone repair, but the exact mechanisms were unclear (Zheng and Wang 2020).

Treg cells activated by A2A receptor agonists also play important roles in many diseases (Lei et al. 2015; Mehta et al. 2018; Proto et al. 2018; Sharabi et al. 2018; Fischer et al. 2019). In a study using Treg cells isolated from septic mice, A2A receptor agonist-upregulated CD39 and CD73 expression was associated with the protection of the Treg cells (Bao et al. 2016b). In allergic asthma or chronic obstructive pulmonary disease (COPD), A2A receptor agonists could increase the ratio of Treg cells in blood and lung tissue to alleviate the illness (Guerrero 2018; Sharabi et al. 2018; Wang et al. 2018). When A2A receptor knockout occurs in Treg cells, Treg cells lose the protective effects in kidney ischemia-reperfusion injury (IRI) (Kinsey et al. 2012). Thus, the A2A receptor agonists or the A2A receptor are vital for Treg cells.

The Gene Expression Omnibus (GEO) database is a public gene expression database from the National Center of Biotechnology Information (NCBI) (Zhang et al. 2019). Bioinformatics is the discipline involving the use of computers for the collection, processing, storage, dissemination, analysis, and interpretation of biological information. We used databases that have been widely used to analyze the function and the regulation of genes, such as GO, KEGG and STRING. However, bioinformatics studies can provide only a preliminary analysis of each mechanism, which then needs to be verified experimentally.

Some studies utilized the computational pipelines to identify the hub genes from RNA-Seq or Microarray data. Fu et al. (2020) used the RNA-seq raw count data from TCGA database to evaluate the immune related genes and pathways for predicting the prognosis of lung squamous cell carcinoma. Udhaya Kumar et al. (2020a) gathered the data from the obesity and co-morbid diseases database to explore the essential signaling cascades and gene networks in diabesity. Wan et al. (2020) analyzed the data from TCGA database to investigate the different expressed lncRNAs, miRNAs and mRNAs in prostate cancer. Udhaya Kumar et al. (2020b) analyzed the GSE13985 dataset from GEO database to detect the differentially expressed genes and molecular pathways in familial hypercholesterolemia involved in atherosclerosis. They used bioinformatic strategies in their current research to help the experimental biologists to further carry forward the outcomes to treat patients.

In this study, after analyzing the genes and pathways of Treg cells regulated by the adenosine A2A receptor as an immunomodulator, we proposed the hypothesis that adenosine A2A receptors regulate multiple Treg cell signal pathways and genes that play an anti-inflammatory function. We retrieved the gene microarray data of Treg cells regulated by adenosine A2A receptor from the GEO database and then used the bioinformatics method to screen the DEGs among wild-type Treg cells stimulated by selective adenosine A2A receptor agonists, Treg cells with adenosine A2A receptor knockout, and wild-type Treg cells. We analyzed the regulation effects of the adenosine A2A receptor on Treg cells and identified key genes and pathways. This research contributes to understanding the mechanisms of adenosine A2A receptor agonists in Treg cells, but the results must be validated by analytical approaches and larger groups before a valid conclusion could be taken.

Materials and methods

Microarray data

GEO DataSets Database series GSE34006 was retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34006) (Barrett et al. 2012). The annotation platform was the GPL1261 Affymetrix Mouse Genome 430 2.0 Array [Mouse430_2]. The data experiment type was expression profiling by array and mainly focused on the function of the adenosine A2A receptor expression in Treg cells. The data experiment design was as follows: CD4⁺CD25⁺ cells were isolated from the spleens and lymph nodes of adenosine A2A receptor^+/+FOXP3^GFP and adenosine A2A receptor^–/–FOXP3^GFP mice using magnetic bead–based kits. Next, FOXP3^GFP Treg cells were sorted by an iCyte Reflection Cell Sorter. Adenosine A2A receptor^+/+ Treg cells were incubated for 18 h with 10 nM ATL-1222 (agonist group) or vehicle (DMSO)(control group). Adenosine A2A receptor^–/– Treg cells were incubated for 18 h with vehicle (DMSO) (knockout group).

Data processing

The microarray data were analyzed by R. A grey-scale map of the chip data, a normalized unscaled standard errors (NUSE) box graph, an RNA degradation map and an MA plot were generated to evaluate the quality of the data by affyPLM, RColorBrewer, affy, and oligo packages of the R software. Second, the data were preprocessed by the RMA method, and the GPL1261 platform was used to convert the probe IDs into the appropriate gene symbols. The missing values were supplemented by the k-Nearest Neighbor (KNN). Last, the DEGs among the groups were analyzed by limma packages of R. The criterions for DEGs were an adjusted P-value of less than 0.05 and a | log₂fold-change(FC) | greater or equal to 1. The code was deposited in the Mendeley data repository (https://doi.org/10.17632/t8b467z34v.1).

GO analysis and KEGG pathway analysis

We used two methods to make GO annotations and KEGG pathways for DEGs: DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/, Version 6.8), an online analysis database, and the BiNGO apps of Cytoscape 3.8.0 software (https://cytoscape.org/) (Huang et al. 2009; Huang da et al. 2009). P-values of less than .05 and gene counts greater than 2 were the criteria for both.

DAVID bioinformatics steps: click ‘start analysis’, paste the list of DEGs, and select the official gene symbol. Then, choose GO and KEGG to analyze. BiNGO apps bioinformatics steps: Open the app in Cytoscape 3.8.0 software, paste the list of DEGs and select Mus musculus. Then, choose GO to analyze.

PPI network construction

The STRING version 11.0 database (https://string-db.org/) was used to construct the DEGs’ PPI networks and to create networks for visualization (Szklarczyk et al. 2017). The confidence score for selection was greater than or equal to 0.4. We hid the disconnected nodes in the network. The hub genes in the PPI networks were determined according to the degree and the betweenness, as calculated by the cytoHubba 0.1 and MCODE of Cytoscape 3.8.0 software (Shannon et al. 2003).

miRNA-hub genes and TFs-hub gene regulatory network construction

The Targetscan online software (http://www.targetscan.org/) was used to predict the regulatory miRNA of hub genes (Garcia et al. 2011). Mouse species was selected and the hub genes was imported, then we got the regulatory miRNA. AnimalTFDB (v3.0) (http://bioinfo.life.hust.edu.cn/AnimalTFDB/) was used to predict the transcription factors (TFs) of hub genes (Hu et al. 2019). The protein sequences of hub genes were acquired in NCBI databases as FASTA format. The protein sequences were imported in AnimalTFDB (v3.0) and we got the TFs.

Results

Differentially expressed genes

In our study, GSE34006 was retrieved and analyzed by R. With the criteria for the DEGs being an adjusted P-value less than .05 and a | log₂fold-change(FC) | greater than or equal to 1, in all, 249 DEGs were selected among groups. In these DEGs, 163 genes were upregulated and 86 genes were downregulated (Figure 1). The top ten upregulated genes were: Tnfrsf9, Gls2, Ramp3, Odc1, Ctla2a, Ctla2b, Rgs2, LOC73899, Pdcd1 and Penk; The top ten downregulated genes were: 4930422I22Rik, Ighv14-2, Ppp1r3fos, Hrh3, 1700031F05Rik, 4930556N13Rik, Plekho2, Wfdc2, Slc9b2 and Srrm3os (Table 1).

Figure 1. Heat map of DEGs (the top 60 genes). Agonist: Adenosine A2A receptor^+/+ Treg cells incubated 18 h with 10 nM ATL-1222; knockout: Adenosine A2A receptor^–/– Treg cells were incubated 18 h with DMSO; wt: Adenosine A2A receptor^+/+ Treg cells were incubated 18 h with DMSO. Colors from white to blue mean increasing expression of DEGs among groups. DEGs, differentially expressed genes.

Figure 2. PPI network of DEGs. Hub genes are labeled by triangles. Red indicates upregulated genes, and green indicates downregulated genes. DEGs, differentially expressed genes.

Table 1 Ten most upregulated and downregulated DEGs.

Gene	logFC	adjusted P-value
Tnfrsf9	3.28	7.72E–04
Gls2	3.18	1.05E–03
Ramp3	2.99	1.05E–03
Odc1	2.70	1.38E–03
Ctla2a	2.63	8.82E–03
Ctla2b	2.61	1.49E–02
Rgs2	2.46	3.53E–04
LOC73899	2.42	3.53E–03
Pdcd1	2.39	1.05E–03
4930422I22Rik	−2.53	1.05E–03
Ighv14-2	−1.68	1.36E–02
Ppp1r3fos	−1.67	2.21E–03
Hrh3	−1.59	2.04E–03
1700031F05Rik	−1.42	1.99E–03
4930556N13Rik	−1.42	3.53E–03
Plekho2	−1.42	1.85E–03
Wfdc2	−1.38	3.34E–03
Slc9b2	−1.36	6.15E–03

GO analysis and KEGG pathways

DAVID Bioinformatics Resources 6.8 and BiNGO apps of Cytoscape 3.8.0 software were used to analyze GO annotation and the enrichment of KEGG pathways for the 249 DEGs. At last, the 249 DEGs were enriched to 75 Go terms (Biological Process (BP) 55; Cellular Component (CC) 7 and Molecular Function (MF), 13) with a P-value less than .05 and gene counts greater than 2 as the criteria. The top 10, or all significant terms with the lowest P-values of BP, CC, and MF are shown in Table 2. Five pathways were enriched for 249 DEG genes, with a P-value less than 0.05 and a gene count greater than 2 as criteria. All DEG gene pathways are shown in Table 3, including six genes involved in TNF signaling pathway (p = 5.70E–03), five genes involved in the TGF-beta signaling pathway (p = 1.20E–02), six genes involved in transcriptional misregulation in cancer (p = 0.03), seven genes involved in cytokine–cytokine receptor interaction (p = 4.20E–02) and four genes involved in pertussis (p = 4.30E–02).

Table 2 GO terms for DEGs (top 10 or all).

Biological process	Gene symbol	Count	P-value	Cellular component	Gene symbol	count	P-value	Molecular function	Gene symbol	count	P-value
response to lipopolysaccharide	Cebpb, Fos, Ccr7, Cxcl10, Dusp10, Junb, Ly96, Penk, Tnfrsf18, Tnfrsf9.	10	1.00E–04	external side of plasma membrane	Cd69, Cd86, Antxr2, Ccr7, Cxcl10, Ly75, Pdcd1, Tnfrsf18, Tnfrsf9.	9	8.10E–03	sequence-specific DNA binding	Bcl6, Cebpb, Dmrtb1, Fos, Jdp2, Pou2f3, Ar, Creb3l2, Crem, Epas1, Foxj1, Fosl2, Junb, Nr4a2, Rel, Stat5a, Sp5, Zfp449.	18	1.50E–04
response to cAMP	Fos, Areg, Crem, Dusp1, Junb, Per1.	6	1.10E–04	Lysosome	Abcb9, 3110043O21Rik, Got1, P2ry2, Ramp3, Slc15a3, Tmbim1, Tmem175, Tnfaip3.	9	1.00E–02	MAP kinase tyrosine/serine/threonine phosphatase activity	Dusp1, Dusp10, Dusp4, Dusp6.	4	2.40E–04
inactivation of MAPK activity	Dusp1, Dusp10, Dusp4, Dusp6.	4	3.40E–04	extrinsic component of cytoplasmic side of plasma membrane	Errfi1, Rgs1, Snx18, Snx9.	4	2.30E–02	RNA polymerase II regulatory region sequence-specific DNA binding	Bcl6, Cebpb, E430018J23Rik, Ascl3, Fosl2, Irf8, Junb, Nr4a2, Rel, Sp5.	10	3.20E–04
inflammatory response	Bcl6, Cd180, Ccr7, Ccrl2, Cxcl10, Il10, Ly96, Rel, Tnfrsf18, Tnfrsf9, Tnfaip3.	11	1.50E–03	endosome membrane	Rab15, Snx18, Tmbim1, Tmem175, Vps37b.	5	2.50E–02	RNA polymerase II core promoter proximal region sequence-specific DNA binding	Cebpb, Fos, H2afz, Jdp2, Pou2f3, Ar, Creb3l2, Foxj1, Fosl2, Junb, Nr4a2, Per1.	12	8.10E–04
response to light stimulus	Fos, Bhlhe40, Dusp1, Per1.	4	1.60E–03	membrane	Abcb1a, Abcb9, Bambi, Cd180, Cd209a, Cd248, Cd69, Cd86, Dnajb14, Errfi1, Fos, Gpr171, Gprc5a, Gem, Mrgprg, Gnptg, Nat8, Nlrx1, Rab15, Ap3b1, Areg, Antxr2, Creb3l2, Casq1, Crbn, Ccr7, Ccrl2, Dgat1, Eea1, Elovl3, Emb, Ern2, Efnb3, Fam187b, Gxylt2, Gypc, Gna15, Hrh3, Has1, Katnb1, Kifc5b, Ly75, Mreg, Mtnr1a, Pth1r, Pdcd1lg2, Pdcd1, Pcsk4, P2ry10, P2ry2, Ramp3, Rgs1, Rgs2, Rnf125, Sctr, Smim23, Smim3, Slc15a3, Slc22a15, Slc25a33, Slc5a3, Slc7a3, Slc9b2, Sstr2, Snx18, Snx9, Syne3, Sytl2, Tmbim1, Tmem151a, Tmem175, Tnfsf8, Tnfrsf18, Tnfrsf9, Upk1a, Vps37b, Vopp1.	77	2.50E–02	transcription factor activity, sequence-specific DNA binding	Bcl6, Cebpb, Dmrtb1, Fos, Jdp2, Pou2f3, E430018J23Rik, Ar, Creb3l2, Crem, Epas1, Foxj1, Fosl2, Id3, Irf8, Junb, Nr4a2, Rel, Stat5a, Zfp449.	20	1.00E–03
negative regulation of B cell activation	Samsn1, Foxj1, Tnfaip3.	3	1.80E–03	transcription factor complex	Fos, Pou2f3, Ascl3, Crem, Epas1, Junb, Rel.	7	3.40E–02	transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	Cebpb, Fos, Pou2f3, Ar, Creb3l2, Epas1, Foxj1, Junb, Nr4a2, Stat5a.	10	1.30E–03
positive regulation of transcription from RNA polymerase II promoter	Cebpb, Fos, H2afz, Pou2f3, Ar, Bmp4, Creb3l2, Crem, Epas1, Foxj1, Fosl2, Il10, Junb, Nr4a2, Per1, Rgcc, Rel, Serpinf2, Stat5a, Tnfsf8.	20	2.40E–03	extracellular exosome	Abcb1a, Cd248, Cd86, Cobll1, Eps8l1, Gprc5a, Gipc2, H2afz, Gnptg, Nat8, Rab15, Wfdc2, Wdr60, Akr1b8, F12, Eea1, Fth1, Got1, Klk1b27, Ly75, Myo1d, Oaf, Pon3, Pth1r, Pdcd1lg2, Serpinf2, Scpep1, Snx18, Snx9, Txndc8, Tmbim1, Tnfaip3, Upk1a, Vps37b.	34	3.60E–02	signal transducer activity	Gpr171, Gprc5a, Mrgprg, Ccr7, Ccrl2, Gna15, Hrh3, Irs2, Mtnr1a, Pth1r, P2ry10, P2ry2, Sctr, Stat5a, Sstr2, Vopp1.	16	1.70E–03
circadian rhythm	Bhlhe40, Crem, Id2, Id3, Mtnr1a, Per1.	6	3.50E–03					RNA polymerase II core promoter sequence-specific DNA binding	Cebpb, Fos, H2afz, Nr4a2, Stat5a.	5	4.10E–03
cellular response to lipopolysaccharide	Cebpb, Cd180, Cd86, Cxcl10, Irf8, Il10, Ly96, Tnfaip3.	8	3.60E–03					G-protein coupled purinergic nucleotide receptor activity	Gpr171, P2ry10, P2ry2.	3	6.90E–03
cellular response to estradiol stimulus	H2afz, Il10, Ramp3, Sstr2.	4	3.60E–03					protein tyrosine/serine/threonine phosphatase activity	Dusp1, Dusp10, Dusp4, Dusp6.	4	6.90E–03

Note: Count, the number of DEGs involved in GO terms.

Table 3 KEGG pathway analysis for DEGs (all).

KEGG Term	Gene symbol	Count	P-value
TNF signaling pathway	Cebpb, Fos, Creb3l2, Cxcl10, Junb, Tnfaip3.	6	5.70E–03
TGF-beta signaling pathway	Bambi, Bmp4, Cdkn2b, Id2, Id3.	5	1.20E–02
Transcriptional misregulation in cancer	Bcl6, Cebpb, Cd86, Id2, Rel, Supt3.	6	3.00E–02
Cytokine-cytokine receptor interaction	Ccr7, Cxcl10, Il10, Gm13308, Tnfsf8, Tnfrsf18, Tnfrsf9.	7	4.20E–02
Pertussis	Fos, Irf8, Il10, Ly96.	4	4.30E–02

Note: Count, the number of DEGs involved in pathway.

PPI network construction

Two hundred and forty-nine DEGs were input into the STRING database, and the disconnected nodes in the network were hidden. Then, a PPI network was constructed, which included 206 nodes and 219 edges (Figure 2). The average node degree was 2.13, and the average local clustering coefficient was 0.323. The P-value of PPI enrichment was less than 1.0E–16. Then, a.tsv file from the STRING database was imported into Cytoscape 3.8.0 software to find the hub genes. We used two kinds of computing methods: degree and betweenness. Then the top ten genes were sorted out by each computing method (Table 4). Finally, we used the intersection of the top ten genes as the hub genes. The four hub genes, Fos, Stat5a, Cebpb and Il10, were upregulated (Table 5).

Table 4 The top 10 genes in PPI network by two kinds computing methods.

Gene	Betweenness	Gene	Degree
Fos	2809.94513	Fos	21
Bmp4	1332.34286	Il10	21
Stat5a	1213.41455	Ccr7	19
Akr1b8	1008	Cd86	17
Il10	993.73108	Cxcl10	16
Gls2	876	Irf8	13
Penk	695.59576	Stat5a	13
Aldh4a1	600	Cebpb	12
Wisp1	600	Cd69	12
Cebpb	555.15629	Bcl6	11

Table 5 The expression levels of these hub genes from the studied dataset.

Hub genes	Agonist 1	Agonist 2	Knockout 1	Knockout 2	wt1	wt2	q value
Fos	6.181985	5.885081	4.340302	4.169498	5.353466	5.145954	0.001845
Stat5a	7.174003	6.679347	5.546582	5.587561	5.983301	5.773283	0.005603
Il10	4.770814	4.229043	2.752944	2.958154	3.008476	2.748858	0.003899
Cebpb	8.015737	7.93418	5.911602	5.774646	6.897759	6.329701	0.002011

miRNA-hub genes and TFs-hub gene regulatory network construction

Thirteen miRNAs were predicted to regulate Fos. 2 miRNAs were regulated Stat5a, 1 miRNA was regulated Cebpb, and 3 miRNAs were regulated Il10. 1 TF was predicted to regulate Fos, 1 TF was predicted to regulate Stat5a, 1 TF was predicted to regulate Cebpb. however, 0 TF was predicted to regulate Il10. (Table 6)

Table 6 miRNA-hub genes and TFs-hub gene regulatory network.

hub gene	miRNAs	TFs
Fos	miR-221-3p, miR-222-3p, miR-29-3p, miR-101-3p.1, miR-142-5p, miR-802-5p, miR-101a-3p.2, miR-101b-3p.1, miR-101b-3p.2, miR-139-5p, miR-181-5p, miR-144-3p, miR-383-5p.1	AAH29814.1
Stat5a	miR-141-3p, miR-200a-3p	CAA88419.1
Cebpb	miR-155-5p	AAA37374.1
Il10	let-7-5p, miR-98-5p, miR-27-3p	no

Discussion

Treg cells are a special type of T cell that inhibits antigen–antibody reactions acting as immunosuppressants or immunoregulators (Andersen 2017; Zhang et al. 2017; Mehta et al. 2018; Sharabi et al. 2018). These cells are vital in regulating inflammatory, immune functions, tissue healing, angiopoiesis and tumor metastasis (Zhang et al. 2017; Lee 2018; Sharabi et al. 2018; Wang et al. 2019; Zheng and Wang 2020). In addition, they play an indispensable role in autoimmune disease, cardiovascular diseases, tumors, COPD, asthma, bone healing, and so on (Lee 2018; Sharabi et al. 2018; Wang et al. 2018; Zheng and Wang 2020). The exploration of the regulation of Treg cells is crucial to advance the treatment of many diseases.

The adenosine A2A receptor is a regulator in various cells, such as endothelial cells, macrophages, osteoclasts, osteoblasts, cancer cells, and Treg cells (Kinsey et al. 2012; Li et al. 2012; Leone et al. 2015; Maj et al. 2017; Guerrero 2018; Zheng and Wang 2020). However, there are few studies that discussed the involved genes and pathways of adenosine A2A receptor on regulatory T cells through the adenosine A2A receptors of Treg cells are essential for cell function (Sitkovsky et al. 2008).

The biological effects of adenosine A2A receptors on regulatory T cells have been studied in many diseases. Kinsey et al. (2012) found that adenosine A2A receptor-deficient Treg cells inhibited the kidney IRI protection of Tregs, while receptor activation by ATL-1222 could enhance Treg cells to suppress kidney IRI. Klein and Bopp (2016) reported that adenosine could strongly improve the suppressive functions of Treg in the course of inflammation by the adenosine A2A receptor. However, the specific genes and pathways of Treg cells regulated by the adenosine A2A receptor have not been reported.

In this integrated study, we used bioinformatic analysis to explore the potential genes and pathways participated the biological effects of adenosine A2A receptor on regulatory T cells as the previous studies (Udhaya Kumar et al. 2019; Yan et al. 2019; Udhaya Kumar et al. 2020c; Mishra et al. 2021; Udhaya Kumar et al. 2021). We analyzed the gene expression changes in Treg cells that are influenced by the adenosine A2A receptor. To accomplish this, GSE34006 was used in the GEO DataSets database. Using the criteria for DEGs (P-value <.05 and | log2fold-change(FC) | ≥1), 249 DEGs were selected. In these DEGs, 163 genes were upregulated, and 86 genes were downregulated.

The regulated biological processes of adenosine A2A receptor on regulatory T cells were performed by GO enrichment analysis. We found that response to lipopolysaccharide, response to cAMP, inactivation of MAPK activity, inflammatory response, response to light stimulus, negative regulation of B cell activation, positive regulation of transcription from RNA polymerase II promoter, circadian rhythm, cellular response to lipopolysaccharide and cellular response to estradiol stimulus were the top 10 most significant biological processes with the lowest P-values. The results were shown that immune and inflammatory regulation were the main biological processes of adenosine A2A receptor on regulatory T cells. However, only five pathways were enriched for the 249 DEGs genes: the TNF signaling pathway, the TGF-beta signaling pathway, transcriptional misregulation in cancer, cytokine–cytokine receptor interaction, and pertussis.

TNF signaling pathway and the TGF-beta signaling pathway were the two pathways that we thought were the main pathways regulated by adenosine A2A receptor on regulatory T cells and these were needed to be verified by ELISA test and knockout cell model.

Tumor necrosis factor (TNF) is an important cytokine in inflammatory reactions and immune system response (Sun et al. 2017; Mehta et al. 2018). It produces pro-inflammatory and anti-inflammatory effects in autoimmune diseases, cancers, tissue healing, and infectious diseases (Mehta et al. 2018). Mehta et al. (2018) summarized that TNF could exert strong protective or pathogenic activity on Treg cells. With co-stimulatory molecules, such as OX40 L and GITRL, TNF can promote the suppressive capacities of Treg cells; furthermore, by itself, TNF can downregulate can downregulate the immunosuppressive functions of Treg cells (Mehta et al. 2018). The adenosine A2A receptor enhances Treg cells’ protection functions, such as anti-inflammation or immunosuppression, possibly by activating TNF and co-stimulatory molecule interactions, such as cytokine–cytokine receptor interaction, in GO terms. Transforming growth factor-beta (TGF-beta) is essential for Treg cells (Bao et al. 2016a; Wang et al. 2019). The application of TGF-beta-blocking antibodies or chemicals leads to the impairment of the Treg cells’ suppressive functions (Toomer and Malek 2018). Treg cells secrete TGF-beta to accelerate tissue healing and alleviate inflammation damage (Maj et al. 2017; Salamon et al. 2017; Lee 2018b; Toomer and Malek 2018). In TGF-beta deficient mice, the expression of forkhead/winged-helix transcription factor 3 (FOXP3) within the Treg cells was low, leading to Treg cell dysfunction (Toomer and Malek 2018). Our findings suggest that the adenosine A2A receptor plays an important regulatory role by affecting the expression TGF-beta, TNF and co-stimulatory molecules.

To assess the degree of interaction between genes, we used the STRING database to build a PPI network. We then used cytoHubba 0.1 and MCODE of the software Cytoscape 3.8.0 to calculate the hub genes using degree and betweenness computations. Finally, we extracted the intersection of the top 10 genes as the hub genes and identified four hub genes, Fos, Stat5 a, Cebpb, and Il10, which were upregulated. The Fox gene regulates Treg cell functions. Bao et al. found that the knockdown of c-Fos (production of the Fox gene) reduced FOXP3 expressions in Treg cells and that induction of c-Fos would elevate the transcription of the FOXP3 gene (2016a, 2016b). Lee (2018a) found that the IRF8 gene controlled the type I immune responses and that IRF8-deficient Treg cells expressed aberrant cytokines or caused CXC chemokine receptor 3 (CXCR 3) deficiency; in addition, they found that IRF8-deficient Treg cells could not migrate to inflammatory sites. Interleukin-10 (IL-10) is an important anti-inflammatory factor in tumor immunity and inflammatory immunity. Treg cells secrete IL-10 to play an immunosuppressive role (Maj et al. 2017; Proto et al. 2018; Wang et al. 2019). The Stat5a gene is critical for peripheral Foxp3 expression and Treg cell maintenance; Stat5a binds to the murine Foxp3 gene, implying a direct role for Stat5 in regulating Foxp3 transcription (Yao et al. 2007).

These results confirm our hypothesis that adenosine A2A receptors regulate multiple signal pathways and genes in Treg cells to perform anti-inflammatory functions. Overall, the results must be validated by analytical approaches and larger groups before any valid conclusion could be taken.

Conclusions

In summary, we used bioinformatics to analyze the DEGs of Treg cells regulated by the adenosine A2A receptor. We obtained 249 DEGs. We then performed GO enrichment, KEGG pathway analysis, and a PPI network of these DEGs. Four hub genes, Fos, Stat5a, Cebpb and Il10, were found to be upregulated. These hub genes may play important roles in Treg cells regulated by the adenosine A2A receptor. This study is based on gene transcription data and must be confirmed by rigorous, molecular biology experimental research.

Abbreviations

Treg cells	=	Regulatory T cells
DEGs	=	differential expression genes
FC	=	fold-change
DAVID	=	database for annotation, visualization and integrated discovery
KEGG	=	Kyoto encyclopedia of genes and genomes
PPI	=	protein–protein interaction
GO	=	gene ontology
TNF	=	tumor necrosis factor
TGF-beta	=	transforming growth factor-beta
MSCs	=	mesenchyma stem cells
NECA	=	n-ethylcarbazole
COPD	=	chronic obstructive pulmonary disease
IRI	=	ischemia-reperfusion injury
GEO	=	gene expression omnibus
NCBI	=	National Center of Biotechnology Information
DMSO	=	dimethyl sulfoxide
NUSE	=	normalized unscaled standard errors
KNN	=	k-nearest neighbor
BP	=	biological process
CC	=	cellular component
MF	=	molecular function
CXCR 3	=	C-X-C chemokine receptor 3
IL-10	=	interleukin-10
TFs	=	transcription factors

Consent for publication

All authors agreed to the publication of this research.

Availability of data and materials

The raw data for this study was obtained from the GEO database [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34006]. The supplemental data used to support the findings of this study have been deposited in the Mendeley Data repository [https://doi.org/10.17632/n4k6wywftg.1].

Code availability

R language: https://www.r-project.org/. (The code was deposited in the Mendeley data repository (https://doi.org/10.17632/t8b467z34v.1)). DAVID Bioinformatics Resources 6.8: https://david.ncifcrf.gov/ Cytoscape 3.8.0 software: https://cytoscape.org/ STRING version 11.0 database: https://string-db.org/ Targetscan software: http://www.targetscan.org/ AnimalTFDB (v3.0): http://bioinfo.life.hust.edu.cn/AnimalTFDB/

Acknowledgements

The authors thank Kinsey GR, okusa MD for uploading database series GSE34006 to(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34006) for public use. The authors are also grateful for the continuous updating of R language and the development of related software used in this research. Jingyi Wang and Xi Zheng designed this study; Dong Wang analyzed GSE34006 data from the GEO database; Dong Wang wrote the article, and Xi Zheng revised the article. All the authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).