Identification of key genes and pathways contributing to artery tertiary lymphoid organ development in advanced mouse atherosclerosis.

Atherosclerosis is a leading cause of mortality worldwide. Artery tertiary lymphoid organ (ATLO) neogenesis is affected by abdominal aorta atherosclerosis, which may lead to an immune response. The present study obtained microarray data to investigate the gene expression differences underlying the potential pathogenesis of atherosclerosis and to elucidate the mechanisms underlying ATLO development. Microarray studies of the aorta, plaques, adventitia, blood, spleen, renal lymph nodes and ATLO were downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) were identified in aorta clusters and ATLO clusters. Kyoto Encyclopedia of Genes and Genomes enrichment and Gene Ontology (GO) analyses were conducted to predict the biological functions of DEGs. The results demonstrated that interleukin 7 receptor (Il7r), C‑X‑C motif chemokine ligand (Cxcl)16, Cxcl13, Cxcl12, C‑C motif chemokine receptor 2, C‑C motif chemokine ligand (Ccl)8, Ccl5 and Ccl12 may function through pathways associated with 'cytokine‑cytokine receptor interaction' and 'chemokine signaling pathway' in ATLO. Gene expression alterations were validated by reverse transcription‑quantitative polymerase chain reaction. Il7r appeared to be the central gene involved in these events, and chemokines and/or chemokine receptors were visualized by GO enrichment. A protein‑protein interaction network was constructed, which suggested that Il7r had a core function in all clusters. Taken together, the results indicated that Il7r upregulation may serve an important role in ATLO development via 'cytokine‑cytokine receptor interaction' and 'chemokine signaling pathway'. This may provide novel perspectives for understanding ATLO development and the regulation of the immune response in atherosclerosis.

Introduction

Atherosclerosis is a chronic inflammatory disease with high morbidity and mortality worldwide (1,2). Rupture-prone vulnerable plaques are a key characteristic of atherosclerosis, caused by disordered lipid metabolism and the inflammatory response (3,4). However, at present, no effective therapeutic strategies are available to treat atherosclerosis. Increasing evidence suggests that aortic adventitia serve a critical role in the immune response in atherosclerosis (5–9). Artery tertiary lymphoid organs (ATLOs), which were initially identified in aged apolipoprotein E (ApoE)−/− mouse adventitia, may be involved in the regulation of the immune response to atherosclerosis (6,7,10–13). Advanced ATLOs consist of distinct immune cell compartments organized into T-cell areas, B-cell follicles and plasma cell niches (7,10). It has been identified that ATLOs recruit naive T cells and generate cluster of differentiation (CD)4+ and CD8+ T cells, consequently inducing Treg and memory B cell production. Similar to lymph nodes, lymph vessels and high endothelial venules in ATLOs specifically facilitate lymphocyte recruitment and migration. This evidence suggests that ATLOs may mediate the immune response to atherosclerosis (11–13). Therefore, it is necessary to focus on the molecular mechanisms of ATLO development for further understanding of ATLOs and their potential impact on atherosclerosis immunity.

Recent advances in microarray methods and data mining tools offer the possibility to simultaneously analyze numerous genes associated with the complex mechanisms of ATLOs. The microarray dataset GSE40156 was used to compare mRNA expression in the aorta, plaque, adventitia and ATLO, in addition to the blood and secondary lymphoid organs (sLOs; spleen and renal lymph nodes). Microarray-based transcriptional profiling was used to evaluate ATLOs and the associated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The identified genes provided insights into the biological processes underlying ATLOs. However, reliably deploying differentially expressed genes (DEGs) into functionally relevant mechanisms is challenging. In the present study, the array data of GSE40156 were downloaded and DEGs associated with KEGG enrichment were verified by bioinformatics methods. Furthermore, a protein-protein interaction (PPI) network was constructed and dissected. Currently, only a limited number of studies concerning ATLOs have been published, to the best of our knowledge. The present study highlighted the DEGs and KEGG pathways involved in the progression of ATLOs, enhancing the understanding of immune responses in atherosclerosis and identifying the candidate genes that may be used for the development of novel diagnostic and therapeutic strategies.

Materials and methods

Affymetrix microarray data

The GSE40156 dataset (12) was downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) to determine and investigate the candidate DEGs that potentially contribute to the progression of ATLOs. Samples of the ApoE−/− and wild-type (WT) aorta, ATLO, plaque, adventitia, blood, spleen and renal lymph node were included in the present study. Additionally, the raw data and annotation files were collected with reference to the GPL1261 (Affymetrix Mouse Genome 430 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and GPL8321 platforms for downstream analysis (Affymetrix Mouse Genome 430A 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc.).

Data processing

The original data were pre-processed with background correction, normalization and expression calculation using the ‘affy’ package (14) in R (version 3.4.2, http://www.R-project.org/). During probe selection for further microarray analysis, non-gene-matched probes were removed, and the probe with the lowest adjusted P-value was selected where one gene matched multiple probes.

Identification of DEGs

The ‘limma’ package (15) in R was used to identify DEGs between the different groups. Following probe annotation, gene expression profiling data were extracted and log2 transformed. The ‘eBayes function’ based on classical Bayesian algorithm (15) was applied to determine DEGs, with P<0.05 and |log2 fold change (FC)|≥1.5 as cut-off criteria.

Functional and pathway enrichment analysis

Gene Ontology (GO) and KEGG pathway analyses were conducted for the DEGs using DAVID (https://david.ncifcrf.gov/). GO classifications consist of molecular function (MF), biological process (BP) and cellular component (CC) terms (16,17). KEGG is a database used to allocate gene sets to their relevant pathways (18). DEGs were further classified by GO term and KEGG pathway enrichment, and visualized using the ‘GOplot’ package in R (16). Adjusted P<0.05 was set as the cut-off criterion.

PPI network construction and module selection

Information on predicted or experimental protein interaction was acquired by means of PPI network construction and analysis, using the STRING database (19,20). This database integrates experiments, databases and text mining data for prediction. The value of interactions between protein pairs was given as a combined score. DEGs were mapped in the database with a cut-off value of a combined score of >0.9. Cytoscape software version 3.6.0 was used to construct the PPI networks (21). A gene network module of the interleukin 7 receptor (Il7r) was constructed to determine its exact biological importance.

Animals

Male ApoE−/− mice (n=10; 3 6 weeks, 21.4±0.5 g; 3 32 weeks, 30.2±0.2 g; and 4 56–60 weeks, 32.6±0.7 g) were obtained from the laboratory animal center of Huazhong University of Science and Technology (Hubei, China) and housed in a controlled environment (20±2°C; 12 h light/dark cycle; 50±10% humidity) with free access to water and food. The animal protocol was reviewed and approved by The Institutional Animal Research Committee of Tongji Medical College (Wuhan, China). All the mice were sacrificed, and the aorta, renal lymph node (RLN), spleen and blood were collected and fresh frozen in liquid nitrogen until required.

Laser capture microdissection (LCM) procedure

Following the harvesting of the tissues, the aorta was perfused in situ using EDTA/PBS for 10 min to remove blood and preserve the adventitia. A dissection microscope was used to remove adipose tissue and para-aortic ganglia but not lymph nodes. Later, the aorta was embedded in Tissue-Tek and stored at −80°C before being cut into sections (10 µm) onto cooled membrane slides. The slides were dried 10 min at 45°C and stored in a vacuum desiccator to avoid moisture uptake. LCM were performed with the Palm MicroBeam laser system (P.A.L.M. Microlaser Technologies AG; Carl Zeiss AG, Oberkochen, Germany) to microdissect the adventitia, media and plaque according to their position of appearance under the LCM microscopy (10). Samples obtained from LCM were later used for reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

RT-qPCR validation

Total RNA was isolated by using RNAiso Plus (Takara Bio, Inc., Otsu, Japan) according to the manufacturer's instructions. Reverse transcription of RNA was performed using PrimeScript™ RT reagent Kit (Takara) in 37°C for 15 min, 85°C for 5 sec and kept in 4°C. Quantitative real-time PCR was performed using PrimeScript™ RT reagent kit (Takara Bio, Inc.) on a StepOnePlus™ Real-time PCR System (Applied Biosystems, Foster City, CA, USA). Thermal cycling conditions were: 30 sec predenaturation at 95°C, followed by 40 cycles of 5 sec denaturation at 95°C and 30 sec annealing at 60°C. Primers were selected using PrimerBank (https://pga.mgh.harvard.edu/primerbank/) and are listed in Table I. All PCR was performed in triplicate, and mRNA fold changes were calculated by the 2− ΔΔCq method (22) using GAPDH as internal reference.

Table I.

Reverse transcription-quantitative polymerase chain reaction primer information.

Gene	Forward	Reverse
Il7r	5′-GCGGACGATCACTCCTTCTG-3′	5′-AGCCCCACATATTTGAAATTCCA-3′
Cxcl12	5′-TGCATCAGTGACGGTAAACCA-3′	5′-CACAGTTTGGAGTGTTGAGGAT-3′
Cxcl13	5′-GGCCACGGTATTCTGGAAGC-3′	5′-ACCGACAACAGTTGAAATCACTC-3′
Cxcl16	5′-ACCCTTGTCTCTTGCGTTCTT-3′	5′-CAAAGTACCCTGCGGTATCTG-3′
Ccr2	5′-ATCCACGGCATACTATCAACATC-3′	5′-TCGTAGTCATACGGTGTGGTG-3′
Ccl5	5′-TTTGCCTACCTCTCCCTCG-3′	5′-CGACTGCAAGATTGGAGCACT-3′
Ccl8	5′-TCTACGCAGTGCTTCTTTGCC-3′	5′-AAGGGGGATCTTCAGCTTTAGTA-3′
Ccl12	5′-ATTTCCACACTTCTATGCCTCCT-3′	5′-ATCCAGTATGGTCCTGAAGATCA-3′
GAPDH	5′-AGGTCGGTGTGAACGGATTTG-3′	5′-TGTAGACCATGTAGTTGAGGTCA-3′

Statistical analysis

Statistical calculations were performed using GraphPad Prism 7 (GraphPad Software, Inc., La Jolla, CA, USA). Experiments were performed in triplicate and data analyzed using an unpaired Student's t-test or one-way analysis of variance with Sidak's multiple comparisons test, and are expressed as the mean ± standard error of the mean. DEGs analyses were performed using one-way ANOVA with Benjamini-Hochberg correction. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of KEGG pathways and differentially expressed genes in 78-week-old mouse aortas

The GSE40156 dataset was acquired for downstream analysis. Hierarchical clustering analysis was used to assess the overall quality of the group data (data not shown). Furthermore, the DEGs between hyperlipidemic mouse groups were verified (78-week ApoE−/− aorta samples vs. 78 week normal WT aorta samples). Following preprocessing, log2-transformed gene expression data from all samples were extracted for analysis. The DEGs between the 78-week ApoE−/− aorta and 78-week wild-type aorta were selected via the classical Bayesian algorithm with |log2FC|≥1.5 and P<0.05. Gene-unmatched probes were removed, and if one gene was matched by multiple probes, only the probe with the lowest adjusted P-value was included. In total, 417 probes (395 upregulated and 22 downregulated, excluding ApoE) were identified (data not shown). As expected, the 78 week ApoE−/− and wild-type aorta was able to be distinguished by hierarchical gene clustering. Subsequently, these 417 probe-matched DEGs were subjected to KEGG enrichment analysis to determine clusters of co-expressed genes sharing the same pathway. The results suggested that the top 10 KEGG pathways significantly enriched by DEGs were ‘chemokine signaling pathway’, ‘phagosome’, ‘staphylococcus aureus infection’, ‘osteoclast differentiation’, ‘tuberculosis’, ‘rheumatoid arthritis’, ‘cytokine-cytokine receptor interaction’, ‘leishmaniasis’, ‘hematopoietic cell lineage’ and ‘B cell receptor signaling pathway’. However, ‘phagosome’, ‘staphylococcus aureus infection’, ‘osteoclast differentiation’, ‘tuberculosis’, ‘rheumatoid arthritis’, ‘leishmaniasis’ and ‘hematopoietic cell lineage’ had no direct association with immunity in atherosclerosis. In addition, ‘B cell receptor signaling pathway’ has previously been investigated (11). Therefore, the present study focused on ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ for further analysis. In total, 41 DEGs were included in these two KEGG pathways (Table II).

Table II.

‘Cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ DEGs.

Group	DEGs
ApoE−/− vs. WT (78-week aorta)	Vav1, Rac2, Ptk2b, Prkcb, Prex1, Pik3r5, Pik3cg, Ncf1, Hck, Gngt2, Dock2, Arrb2, Tnfrsf1b, Tnfrsf13b, Tnfrsf11b, Ltb, Lepr, Il7r, Il6, Il10ra, Cxcr4, Cxcl9, Cxcl16, Cxcl14, Cxcl13, Cxcl12, Cxcl10, Cxcl1, Csf2rb, Csf2ra, Ccr5, Ccr2, Ccr1, Ccl9, Ccl8, Ccl7, Ccl6, Ccl5, Ccl3, Ccl2, Ccl12
ATLO vs. plaque	Itk, Ccr9, Adcy3, Tnfsf13b, Tnfrsf17, Tnfrsf13c, Tnfrsf13b, Tnfrsf12a, Tnfrsf11b, Tgfb2, Pdgfc, Pdgfa, Ngfr, Ltb, Lifr, Lep, Inhbb, Il7r, Il18r1, Hgf, Ghr, Egfr, Cxcr6, Cxcr5, Cxcl9, Cxcl16, Cxcl13, Cxcl12, Csf2ra, Ccr6, Ccr2, Ccl8, Ccl5, Ccl4, Ccl3, Ccl2, Ccl12, Ccl11
ATLO vs. ApoE−/− adventitia	Vav1, Rac2, Prkcb, Pik3cg, Itk, Hck, Dock2, Arrb2, Tnfsf13b, Tnfrsf18, Tnfrsf13c, Tnfrsf13b, Ppbp, Ltb, Il7r, Il6, Il21r, Il1r1, Il1b, Il18r1, Il17ra, Il13ra1, Il10ra, Cxcr6, Cxcr5, Cxcr3, Cxcl16, Cxcl14, Cxcl13, Cxcl12, Cxcl10, Cxcl1, Csf2rb, Clcf1, Ccr6, Ccr5, Ccr2, Ccr1, Ccl8, Adcy7, Ccl5, Ccl12
Common genes	Tnfrsf13b, Ltb, Il7r, Cxcl16, Cxcl13, Cxcl12, Ccr2, Ccl8, Ccl5, Ccl12

DEGs, differentially expressed genes; ATLO, artery tertiary lymphoid organ; WT, wild-type; ApoE, apolipoprotein E.

Common KEGG pathways and DEGs in aorta and ATLO clusters

ATLOs formed in the plaque-adjacent aorta adventitia in aged ApoE−/− mice with atherosclerosis. During atherosclerosis progression and ATLO formation, the adventitia is reorganized, recruiting lymphocytes and other immune cells, including B cells, to generate a dense immune cell aggregate (10). The effects of ATLOs in atherosclerosis progression have previously been demonstrated in hyperlipidemic mice (11–13). Furthermore, to assess the associated KEGG pathways in ATLOs, microarray data of ATLO, plaque and adventitia, separated from the aorta via laser capture microdissection technique (10), were pooled to analyze transcript atlases. The same DEG and KEGG analysis methods were utilized in ATLO plaque and adventitia clusters to determine which genes and pathways were co-expressed (data not shown). In the final analysis, the top 15 KEGG pathways of aorta clusters (78-week ApoE−/− aorta and WT aorta clusters) and ATLO clusters (ATLO plaque and adventitia clusters) were selected, and it was verified that they shared four common KEGG pathways (Table III). Among them, ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ were the top two pathways. Therefore, it was hypothesized that they were involved in the core mechanism of ATLO development. Sets of candidate DEGs were obtained from these two pathways (Table II). Heatmaps of these DEGs in aorta clusters, ATLO plaque clusters and ATLO adventitia clusters are presented in Fig. 1. In total, 10 DEGs were shared between these groups: TNF receptor superfamily member 13B, lymphotoxin-β, Il7r, C-X-C motif chemokine ligand (Cxcl)16, Cxcl13, Cxcl12, C-C motif chemokine receptor 2 (Ccr2), C-C motif chemokine ligand (Ccl)8, Ccl5 and Ccl12. Lötzer et al (23) previously reported that tumor necrosis factor receptor superfamily member 1A/lymphotoxin β-receptor cross-talk is involved in lymphorganogenic chemokine protein secretion and tertiary lymphoid organogenesis. Il7r, Cxcl16, Cxcl13, Cxcl12, Ccr2, Ccl8, Ccl5 and Ccl12 were identified as candidate genes and may help to elucidate the mechanism of ATLO formation.

Table III.

Aorta cluster and ATLO cluster KEGG enrichment.

Group	KEGG pathways
ApoE−/− vs. WT (78-week aorta)	‘Chemokine signaling pathway’; ‘Phagosome’; ‘Staphylococcus aureus infection’; ‘Osteoclast differentiation’; ‘Tuberculosis’; ‘Rheumatoid arthritis’; ‘Cytokine-cytokine receptor interaction’; ‘Leishmaniasis’; ‘Hematopoietic cell lineage’; ‘B cell receptor signaling pathway’; ‘Natural killer cell mediated cytotoxicity’; ‘Leukocyte transendothelial migration’; ‘Lysosome’; ‘Antigen processing and presentation’; ‘TNF signaling pathway’.
ATLO vs. plaque	‘Cytokine-cytokine receptor interaction’; ‘Hematopoietic cell lineage’; ‘Rheumatoid arthritis’; ‘Complement and coagulation cascades’; ‘Intestinal immune network for IgA production’; ‘Phagosome’; ‘PI3K-Akt signaling pathway’; ‘Metabolic pathways’; ‘Lysosome’; ‘ECM-receptor interaction’; ‘Primary immunodeficiency’; ‘Focal adhesion’; ‘Staphylococcus aureus infection’; ‘Malaria’.
ATLO vs. ApoE−/− adventitia	‘Metabolic pathways’; ‘Microbial metabolism in diverse environments’; ‘Carbon metabolism’; ‘Citrate cycle (TCA cycle)’; ‘Propanoate metabolism’; ‘Cytokine-cytokine receptor interaction’; ‘Valine, leucine and isoleucine degradation’; ‘Chemokine signaling pathway’; ‘Hematopoietic cell lineage’; ‘Non-alcoholic fatty liver disease (NAFLD)’; ‘Staphylococcus aureus infection’; ‘Fatty acid metabolism’; ‘B cell receptor signaling pathway’; ‘Osteoclast differentiation’; ‘Leishmaniasis’.
Common pathways	‘Cytokine-cytokine receptor interaction’; ‘Chemokine signaling pathway’; ‘Hematopoietic cell lineage’; ‘Staphylococcus aureus infection’.

KEGG, Kyoto Encyclopedia of Genes and Genomes; ATLO, artery tertiary lymphoid organ; WT, wild-type; ApoE, apolipoprotein E.

Figure 1.

Transcript maps reveal the DEGs in the Kyoto Encyclopedia of Genes and Genomes enrichment terms ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ in aorta and ATLO clusters. (A) Heatmaps of DEGs in the aorta clusters. (B) ATLO-plaque clusters. (C) ATLO-adventitia clusters. The color scale of the raw Z-scores represents high (red) and low expression (green). DEGs, differentially expressed genes; ATLO, artery tertiary lymphoid organ; ApoE, apolipoprotein E.

Furthermore, Il7r, Cxcl16, Cxcl13, Cxcl12, Ccr2, Ccl8, Ccl5 and Ccl12 expression in the aorta, ATLO, plaque and adventitia samples was determined by RT-qPCR (Fig. 2A and B). All of these eight DEGs displayed stable, gradual upregulation with age in ApoE−/− mouse aortas. Expression was significantly different between the 78-week ApoE−/− and 78-week WT aortas, and between ATLO and plaque or between ATLO and ApoE−/− mouse adventitia. To the best of our knowledge, notable atherosclerosis formation occurs at 32 weeks (7). Meanwhile, ATLO neogenesis is affected by atherosclerosis (10). Therefore, Il7r, Cxcl16, Cxcl13, Cxcl12, Ccr2, Ccl8, Ccl5 and Ccl12 may be key genes involved in atherosclerosis-associated ATLO development.

Figure 2.

Reverse transcription-quantitative polymerase chain reaction analysis of eight identified DEGs in the aorta and ATLO clusters. (A) DEGs in the aorta from WT and ApoE−/− mice at 6, 32 and 78 weeks. *P<0.05, **P<0.01 vs. WT adult; #P<0.05, ##P<0.01 vs. WT aged. NS, not significant. (B) DEGs in the ATLO with plaque and adventitia. *P<0.05, **P<0.01 vs. ATLO. (C) Il7r in ATLO with RLN, spleen and blood. NS, not significant. DEGs analyses were performed using one-way ANOVA with Benjamini-Hochberg correction. DEGs, differentially expressed genes; ATLO, artery tertiary lymphoid organ; Il7r, interleukin 7 receptor; WT, wild-type; RLN, renal lymph node; Cxcl, C-X-C motif chemokine ligand; Ccl12, C-C motif chemokine ligand 12; Ccr, C-C motif chemokine receptor; ApoE, apolipoprotein E.

GO enrichment analyses and PPI network

To explore the biological functions involved in ATLOs, GO analysis was performed (Fig. 3). GO terms with an adjusted P<0.05 were considered significantly enriched. DEGs from aorta clusters, ATLO plaque clusters and ATLO adventitia clusters were significantly enriched in the GO terms ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’. The top 15 biological process (BP), eight molecular function (MF) and 12 cellular component (CC) GO terms were visualized in Fig. 3A, C and E.

Figure 3.

GO term analysis illustrating DEGs in the aorta and ATLO clusters. (A) The bubble plot displays GO research based on DEGs of ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ between 78 week ApoE−/− and wild-type aorta. The y-axis indicates the significance of the term (-log10 adjusted P-value) and the x-axis indicates the Z-score. Bubbles indicate the GO terms, with green indicating BP terms, red indicating cellular component terms and blue indicating molecular function terms. The bubble size indicates the gene numbers in the GO terms. (B) Circular plots illustrating interleukin 7 receptor-associated BP GO terms in aorta clusters. The genes are bridged by ribbons to their assigned BP terms. The dark-to-light blue color indicates the logFC. (C) The bubble plot displays GO research based on DEGs of ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ between ATLO and plaque. (D) Circular plots illustrating interleukin 7 receptor-associated BP GO terms in ATLO-plaque clusters. (E) The bubble plot displays GO research based on DEGs of ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ between ATLO and adventitia. (F) Circular plots illustrating interleukin 7 receptor-associated BP GO terms in ATLO-adventitia clusters. DEGs, differentially expressed genes; ATLO, artery tertiary lymphoid organ; GO, Gene Ontology; BP, biological process; FC, fold change.

To the best of our knowledge, the association between Il7r and tertiary lymphoid organogenesis has not been previously described. The Il7r gene encodes the interleukin 7 receptor, which is a protein found on the surface of various cells and has been suggested to serve a critical role in development of Peyer's patches and certain peripheral lymph nodes (24,25). Thus, the present study focused on Il7r-associated GO terms. The GO categories of BP associated with Il7r in aorta clusters, ATLO-plaque clusters and ATLO-adventitia clusters are presented in Fig. 3B, D, F and Table IV. The common BP GO terms were: ‘cytokine-mediated signaling pathway (GO:0019221)’ and ‘cellular response to chemical stimulus (GO:0070887)’; thus, it was hypothesized that these were the principal Il7r-associated biological functions in ATLO development.

Table IV.

Interleukin 7 receptor-associated BP GO terms.

Group	GO term
Aorta clusters	‘Cytokine-mediated signaling pathway’; ‘Immune system process’; ‘Cellular response to chemical stimulus’.
ATLO-plaque clusters	‘Cytokine-mediated signaling pathway’; ‘Cell surface receptor signaling pathway’; ‘Response to cytokine’; ‘Cellular response to chemical stimulus’; ‘Cellular response to cytokine stimulus’; ‘Cellular response to organic substance’; ‘Signal transduction’.
ATLO-adventitia clusters	‘Cytokine-mediated signaling pathway’; ‘Cellular response to cytokine stimulus’; ‘Response to cytokine’; ‘Cellular response to chemical stimulus’; ‘Cell surface receptor signaling pathway’; ‘Positive regulation of immune system process’; ‘Cellular response to organic substance’; ‘Immune system process’.
Common BP GO terms	‘Cytokine-mediated signaling pathway’; ‘Cellular response to chemical stimulus’.

BP, biological process; GO, Gene Ontology; ATLO, artery tertiary lymphoid organ.

PPI network analysis (Fig. 4) was also used for the purpose of characterizing core processes at the protein level. As illustrated in Fig. 4A, C and E, the majority of the nodes were included in ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ in three clusters. For further analysis, Il7r modules were selected from the PPI network. The results demonstrated that the protein-coding genes in the network interacted with Il7r (Fig. 4B, D and F). In the aorta clusters, Il7r interacted with interleukin (Il)6, interleukin 10 receptor subunit-α (Il10ra), protein tyrosine kinase 2β, Ccr5, colony stimulating factor 2 receptor α-subunit (Csf2ra), Cxcr4, Ccl5 and Cxcl13; in ATLO plaque clusters, Il7r interacted with leptin, Cxcr5, IL2 inducible T cell kinase (Itk), Csf2ra, Ccr9, Ccl5, Cxcl13 and Cxcr6; and in ATLO-adventitia clusters, Il7r interacted with Il6, Tnfrsf18, Cxcr5, Il10ra, Cxcr3, Ccr5, Itk, Ccl5, Cxcl13 and Cxcr6. Ccl5 and Cxcl13 were demonstrated to be involved in all of these three clusters with Il7r. Notably, it was also observed that Ccl5 and Cxcl13 were associated with the GO terms ‘cytokine-mediated signaling pathway’ and ‘cellular response to chemical stimulus’, important Il7r-associated BP GO terms. According to previous studies, Ccl5 and Cxcl13 are lymphorganogenic chemokines that trigger the development of elaborate bona fide ATLOs (10,23). It was subsequently deduced that Il7r may be involved, with Ccl5 and Cxcl13, in cytokine-mediated signaling pathways and the cellular response to chemical stimuli in the progression of ATLO development.

Figure 4.

PPI network of DEGs in the aorta and ATLO clusters. (A) PPI network showing differentially expressed genes in aorta clusters. The size of a node is proportional to the number of connections. (B) PPI network of Il7r in aorta clusters. (C) PPI network of DEGs in ATLO plaque clusters. (D) PPI network of Il7r in ATLO plaque clusters. (E) PPI network of DEGs in ATLO adventitia clusters. (F) PPI network of Il7r in ATLO adventitia clusters. The nodes in the network represent proteins (indicated by gene names). Lines demonstrate the associations between proteins. PPI, protein-protein interaction; ATLO, artery tertiary lymphoid organ; Il7r, interleukin 7 receptor.

Il7r is not involved in ATLO spleen/blood/RLN clusters

In contrast, spleen-, blood- and RLN-transcript atlases were noticeably larger, with a more pronounced level of DEGs in ATLOs (data not shown). Analysis of spleen-, blood- and RLN-regulated genes revealed that 599 common DEGs remained (Fig. 5A). The enrichment KEGG pathways of the 599 DEGs were ‘hematopoietic cell lineage’, ‘rheumatoid arthritis’, ‘complement and ‘coagulation cascades’, ‘intestinal immune network for IgA production’, ‘osteoclast differentiation’, ‘leishmaniasis’, ‘proteoglycans in cancer’ and ‘viral myocarditis’. ‘Cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ were not identified. The transcript atlas demonstrated no significant alterations in Il7r expression in ATLOs vs. spleen, blood, or RLN. To confirm these results, GO terms and KEGG enrichment were also separately applied for ATLO-RLN clusters (Fig. 5B-F). It is noteworthy that the ATLO-RLN clusters still performed a comparable enrichment in ‘cytokine-cytokine receptor interaction’ (adjusted P=6.34×10−10) and ‘chemokine signaling pathway; (adjusted P=4.54×10−6). However, they were less important in the ATLO-RLN clusters and Il7r was not identified in the KEGG pathways, indicating that the expression level of Il17r in ATLOs was no different compared with that in RLN, which was also validated by the RT-qPCR (Fig. 2C). As Il7r serves a vital role in lymph node development, it was hypothesized that Il7r may be involved in lymphoid organ neogenesis in ATLOs and RLN, meaning that there were no identified differences between ATLOs and RLN. Similar results were found in ATLO-spleen/blood clusters (data not shown).

Figure 5.

Functional analysis showing DEGs in ATLO-spleen/blood/RLN clusters. (A) Venn diagram representation of the comparative analysis of DEGs in ATLO spleen/blood/RLN clusters. (B) Heatmaps of DEGs of ‘cytokine-cytokine receptor interaction’ in ATLO RLN clusters. (C) Heatmaps of DEGs of ‘chemokine signaling pathway’ in ATLO RLN clusters. The color scale of the raw Z-scores represents high (red) and low (green) expression. (D) The bubble plot displays GO enrichment based on DEGs of ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ between ATLO and RLN. The y-axis indicates the significance of the term (-log10 adjusted P-value) and the x-axis indicates the Z-score. Bubbles indicate the GO terms, with green indicating the BP terms, red the cellular component terms and blue the molecular function terms. Bubble size indicates the gene numbers in the GO terms. (E) Circular plot showing Il7r-associated BP GO terms in ATLO RLN clusters. The genes are bridged by ribbons to their assigned BP terms. The dark-to-light blue color indicates the logFC. (F) PPI network based on these DEGs between ATLO and RLN. The size of a node is proportional to the connection degree. The nodes in network represent proteins (indicated by gene names). Lines demonstrate the associations between proteins. PPI, protein-protein interaction; DEGs, differentially expressed genes; ATLO, artery tertiary lymphoid organ; Il7r, interleukin 7 receptor; BP, biological process; FC, fold change; RLN, renal lymph node; spl, spleen.

Discussion

It is well known that the immune system is crucially involved in atherogenesis. However, there remain questions as to how, when and where cell subset interactions are organized in disease development. The outermost connective tissue that covers arteries is highly active and complex. It is termed adventitia and functions in artery maintenance and homeostasis (5,7). During aging, the quantity of inflammatory leukocytes, including monocytes, macrophages and T-cells, are progressively increased in the adventitia (26), which drives the development of ATLOs from the appearance of atherosclerotic plaques (10). ATLO neogenesis occurs from 32 weeks and is noted in 75% of 78-week-old mice with atherosclerosis (10). They are complex plaque-associated lymphoid organs associated with the immune response, and are crucial in atherosclerosis.

In the present study, ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ were identified from aorta and ATLO clusters. These two pathways are of major importance in guiding immune cells to and within lymphoid and non-lymphoid tissues (27–29). While the development and architectural organization of ATLOs requires homeostatic chemokines or lymphoid cytokines, Il7r appears to be a key gene linking these events with chemokines or chemokine receptors. Previous studies have demonstrated that Il7r is involved in lymph node development and inflammatory disease (24,30). Consistent with this, in the present study, Il7r exhibited significant differential expression in the aged ApoE−/− aorta and ATLOs, compared with WT mice. A limited number of studies have demonstrated the direct association between Il7r and atherosclerosis. However, it has been reported that Il7r is upregulated in atherosclerotic tissues, and Il7r may lead to atherosclerosis through the inflammatory pathway (10,31), and is also associated with a number of other inflammatory diseases such as ulcerative colitis (30), rheumatoid arthritis (32) and multiple sclerosis (33). This evidence suggests that Il7r may also be required for the progression of atherosclerosis. It is known that the formation of ATLOs is influenced by plaques and that their sizes increase with plaque size (10). Thus, it was hypothesized that Il7r may be involved in atherosclerosis-associated ATLO development.

In addition to the aforementioned Il7r gene, the chemokines Cxcl16, Cxcl13, Cxcl12, Ccr2, Ccl8, Ccl5 and Ccl12 were observed to be candidate genes in ATLO formation, and a number of these have already been identified in previous studies (10,23,34). Furthermore, it was identified that they were enriched with Il7r in ‘cytokine-cytokine receptor interaction’ and ‘chemokine signaling pathway’ in the present study. In the PPI network analysis, it was demonstrated that Il7r was also co-expressed with Ccl5 and Cxcl13 at the core of cellular processing at the protein level. Ccl5 is chemotactic for T cells, eosinophils and basophils, triggering the recruitment of leukocytes into inflammatory sites (35). Chemokine Cxcl13 is essential for lymph node initiation (36). In the GO enrichment analysis, Ccl5 and Cxcl13 were also involved in the two major Il7r-associated BP GO terms, namely ‘cytokine-mediated signaling pathway’ and ‘cellular response to chemical stimulus’. These two GO terms may be associated with ATLO development. Combining the results from our analysis, it was concluded that Il7r may contribute to chemokine recruitment in atherosclerosis-associated ATLO development.

The present research was a re-analysis based on available microarray data, rather than the reproduction of conclusions presented by the original publications. The original study also used immunostaining and paid more attention to gene and biological process-associated T/B cell immune processes (10). In the present study, the DEGs and KEGG pathways in the gene expression profile throughout the development of ATLOs were focused upon. Il7r acted as a core component in ATLO lymphogenesis in a hyperlipidemic state. This may contribute to the understanding of atherosclerotic immunity in the adventitia.

Nevertheless, certain limitations, such as a lack of verification at the experimental level, exist in the present study. Complementary experiments in ATLOs will be performed in the future for further confirmation.