Comparative Systems Analyses Reveal Molecular Signatures of Clinically tested Vaccine Adjuvants

A better understanding of the mechanisms of action of human adjuvants could inform a rational development of next generation vaccines for human use. Here, we exploited a genome wide transcriptomics analysis combined with a systems biology approach to determine the molecular signatures induced by four clinically tested vaccine adjuvants, namely CAF01, IC31, GLA-SE and Alum in mice. We report signature molecules, pathways, gene modules and networks, which are shared by or otherwise exclusive to these clinical-grade adjuvants in whole blood and draining lymph nodes of mice. Intriguingly, co-expression analysis revealed blood gene modules highly enriched for molecules with documented roles in T follicular helper (TFH) and germinal center (GC) responses. We could show that all adjuvants enhanced, although with different magnitude and kinetics, TFH and GC B cell responses in draining lymph nodes. These results represent, to our knowledge, the first comparative systems analysis of clinically tested vaccine adjuvants that may provide new insights into the mechanisms of action of human adjuvants.

Scientific RepoRts | 6:39097 | DOI: 10.1038/srep39097 to mycobacterial cord factor (trehalose dimycolate), and acting through mincle receptor activation 8,9 . IC31 is another two-component adjuvant consisting of the cationic peptide KLK combined with the Toll-like receptor 9 (TLR9)-stimulatory oligodeoxynucleotide ODN1a 10 . Glucopyranosyl lipid adjuvant-stable emulsion (GLA-SE) is a synthetic TLR4 agonist in a squalene-based oil-in-water emulsion 11 . Aluminum-based mineral salts have been considered the adjuvant of choice for vaccination against infections that can be prevented by an antibody (Ab) response and have been successfully used in several licensed vaccines 12,13 .
As a model antigen for the present study we chose H56, a clinical-grade multistage tuberculosis (TB) vaccine candidate under clinical evaluation (NCT01865487, NCT02375698, NCT02375698, NCT02503839, NCT01967134, NCT02378207). H56 is a fusion protein consisting of two well-characterized early secreted antigens, antigen 85B (Ag85B) and early secreted antigen 6 (ESAT-6) and one latency-associated antigen Rv2660c from Mycobacterium tuberculosis 14 . This vaccine candidate has been shown to provide protective immunity in mice, in combination with the novel adjuvants tested in the present study (CAF01, IC31 and GLA-SE) 7 . It has furthermore provided protective immunity against TB in cynomolgus macaques in combination with IC31 15 , and CAF01 16 . H56 is in clinical development in combination with IC31 17 .
Using a genome-wide transcriptomics analysis of whole blood (WB) and draining lymph nodes (dLNs) at early time points after immunization combined with a systems biology approach, we identified both common and unique transcriptional profiles of the four adjuvants tested. Transcriptional changes associated with the induction of adaptive immune responses were commonly detected in WB for all adjuvants, whereas co-expression analysis identified distinct gene modules containing adjuvant-unique molecules that have previously been shown to be involved in germinal center (GC) reactions and/or the skewing of T helper cell responses. Guided by the co-expression analysis data, we therefore studied the ability of the four adjuvants to induce GC formation in the draining LNs following immunization in mice. These data showed that all of these adjuvants, when given in combination with H56 antigen, enhanced TFH and GC B cell responses compared with those of the H56 antigen alone, although kinetics and cell frequency differed among the different adjuvanted groups.
The results reported herein are, to our knowledge, the first comparative systems biology analysis of clinically tested vaccine adjuvants and thus could inform the rational development of vaccine adjuvants for human use.

Transcriptomic profiles in WB and dLNs of mice induced by different clinically tested adjuvants.
High-quality RNA samples extracted from WB and the dLNs of mice at four different time points (6 h, 24 h, 72 h and 168 h) after a subcutaneous (s.c.) injection of clinical-grade H56 with or without each clinical-grade adjuvant were subjected to genome-wide microarray analysis.
Mice receiving H56 alone were used as a control group for the groups receiving H56 in combination with any of the four adjuvants. Differentially expressed genes (DEGs) were thus defined as those with a significant (adjusted P value < 0.05 FDR) change in expression in an adjuvant group compared with H56 alone. Alum did not induce any appreciable changes in the transcripts from WB, except for the 72 h time point, where 89 DEGs were detected, most of which were shared by other adjuvants (Fig. 1a,b). Although the kinetics and magnitude of the transcriptional changes varied widely among the groups, the other three adjuvants (CAF01, IC31 and GLA-SE) induced more dramatic changes in the WB transcriptional profile than Alum. GLA-SE was the only adjuvant that had already induced significant transcriptional changes (5,858 DEGs) at 6 h, with the highest number of DEGs (11,460 DEGs) at 72 h (Fig. 1a). CAF01 induced a later onset (24 h) of gene expression in WB (1,152 DEGs), and the highest number of significantly changed transcripts was observed at 72 h (4,490 DEGs) (Fig. 1a). The changes in the WB transcripts induced by IC31 were first observed and peaked at 24 h (3,191 DEGs) and decreased at 72 h (1,031 DEGs) (Fig. 1a). A higher proportion of downregulated transcripts was observed at the two early time points (6 h and/or 24 h), and greater numbers of upregulated transcripts were detected at 72 h for all adjuvants (Fig. 1c).
Next, we studied global changes in transcript profiles of dLNs in response to the different adjuvants. Overall, the transcriptional changes detected in dLNs were modest compared to those in WB. Alum induced limited changes in dLNs at 6 h. This number gradually increased at 24 and 72 h before peaking at 168 h (171 DEGs) (Fig. 1d). GLA-SE induced early changes in gene expression in dLNs at 6 h (1,787 DEGs), which increased at 24 h (9,440 DEGs) and 72 h (10,415 DEGs), followed by a decline at 168 h (5,420 DEGs), although the number of DEGs remained relatively high (Fig. 1d). Similar to the transcriptomic profile of WB, both CAF01 and IC31 induced much smaller numbers of DEGs in dLNs, with a later onset compared with GLA-SE. CAF01 caused changes in gene expression mainly at 24 h (460 DEGs) and only limited changes at 6 h, 72 h and 168 h (Fig. 1d). Similarly, the IC31-induced gene expression in dLNs peaked at 24 h (624 DEGs), decreased at 72 h, and then increased again at 168 h (Fig. 1d). Few DEGs were observed in dLNs of all adjuvant-treated groups, although at 24 h, GLA-SE shared 288 DEGs and 315 DEGs with IC31 and CAF01, respectively (Fig. 1e). The ratio between the up-and downregulated DEGs in dLNs showed similar trends for the different adjuvants across the different time points (Fig. 1f).
Overall, the adjuvants studied herein differed in terms of the magnitude and kinetics of the transcriptional changes induced in both WB and dLNs. In dLNs, GLA-SE induced the highest levels and the earliest onset of transcriptional changes, whereas both CAF01 and IC31 induced later and less detectable changes in the transcriptional profile compared with GLA-SE. Alum induced the most subtle transcriptional changes in both WB and dLNs.
Gene functional patterns of WB revealed common and unique features of clinically tested adjuvants. To identify patterns of biological processes that were affected by the four adjuvants in WB, gene ontology (GO) enrichment mapping was performed for the DEGs of all adjuvants at the different time points. Three GO term groups were selected based on the clustering dendrograms, as explained in Materials and Methods. Distinct areas of enriched GO terms with upregulated DEGs were observed for all four adjuvants Scientific RepoRts | 6:39097 | DOI: 10.1038/srep39097 at 72 h ( Supplementary Fig. 1). Grouping of these GO terms revealed that most belonged to cellular, metabolic and immune system processes, as well as to responses to stimuli (Fig. 2a). Another distinct group of GO terms belonging to biological regulation and metabolic and cellular processes was significantly enriched in the majority of the upregulated DEGs at 72 h in the CAF01, IC31 and GLA-SE groups, but to a lesser extent in the Alum group (Fig. 2b). Many of the GO terms that were commonly changed in WB (mainly at 72 h) by the three novel adjuvants were related to the initiation of adaptive immune responses (e.g., positive regulation of adaptive immune responses, positive regulation of B cell proliferation, and antigen processing and presentation).
Large numbers of GO terms were significantly enriched at early time points in the CAF01 (24 h and decreasing by 72 h) and GLA-SE groups (6 h and gradually decreasing by 168 h), with the majority of the DEGs being upregulated. These GO terms were not significantly enriched in the IC31 group at 6 h or 24 h, and they mainly contained downregulated DEGs at 72 h. The majority of these GO terms belonged to biological regulation and response to stimuli associated with innate and inflammatory functions, such as positive regulation of monocyte chemotactic protein-1 production (MCP-1) and positive regulation of interleukin (IL)-6 and IL-1β production (Fig. 2c).
Induction of systemic responses of cytokines and chemokines by clinically tested adjuvants. Next, we compared the patterns of cytokines and chemokines in sera of mice after administration of the adjuvants using a 23-plex luminex assay, including many known pro-and anti-inflammatory mediators. In line with the whole blood transcriptomics data, we observed the most subtle and greatest responses of cytokines and chemokines by Alum and GLA-SE, respectively (   A directional gene set enrichment analysis of the gene ontology (GO) terms in the DEGs was performed using the PIANO package. At 72 h, a group of GO terms was identified as commonly upregulated by all four adjuvants (a), and another group of GO terms was commonly upregulated by CAF01, IC31 and GLA-SE (b). A third group of GO terms was commonly upregulated by CAF01 and GLA-SE at early time points, whereas it was heavily downregulated by IC31 at 168 h (c). The heatmap shows the value of the enrichment score (− log10(enrichment P value)). The GO terms are color-coded to illustrate the direction of the gene expression changes in the majority of genes included in the respective GO term, i.e., red for upregulation and blue for downregulation.
Scientific RepoRts | 6:39097 | DOI: 10.1038/srep39097 (FC 5.24) and was the only adjuvant inducing high levels of IL-9 (FC 13.57). Taken together, these results, in line with the transcriptomics data, show that GLA-SE induces early and strong pro-inflammatory responses, and that CAF01 and IC31 mount delayed pro-inflammatory cytokine responses of smaller magnitude compared to those of GLA-SE.

Co-expression analysis in WB revealed shared and unique features of clinically tested adjuvants.
Next, a weighted gene co-expression network analysis of DEGs across all time points for each adjuvant was conducted to identify the most highly affected gene modules (groups of highly co-expressed genes). Co-expression analyses utilize interaction patterns among significantly changed genes to reduce high-dimensional microarray datasets in an unbiased way. The co-expression analysis of DEGs in WB identified six different co-expressed modules for CAF01, IC31 or GLA-SE; each module was assigned a color code (Fig. 4a). The modules contained various numbers of genes, ranging from 77 genes (GLA-SE red module) to 566 genes (GLA-SE turquoise module). No co-expression module was identified for Alum due to its limited number of DEGs. Next, gene lists were generated from each module, and Ingenuity Pathway Analysis (IPA) was employed to identify enrichments in known biofunctions. The significantly changed biofunctions for each module and the associated DEGs are listed in Supplementary Table 1 Whereas biofunctions associated with innate immune responses were scattered across different modules of the IC31 group (Supplementary Table 1), both GLA-SE and CAF01 had distinct modules containing molecules involved in innate inflammatory and phagocytosis responses (brown for CAF01 and red for GLA-SE, Figs 4b and f). For both adjuvants, the majority of the biofunctions were associated with cellular movement, such as leukocyte migration, chemotaxis of myeloid cells, chemotaxis of phagocytes and cell movement of neutrophils.
Macrophage-inducible C-type lectin (Mincle) transcripts (Clec4e) were noted in many of the most highly significant innate biofunctions for both CAF01 (brown module) and GLA-SE (red module) (Supplementary Table 1). Mincle was first described as being induced upon lipopolysaccharide stimulation 18 and has recently been reported as a key receptor for the CAF01 19 .
The brown module with the most highly significant biofunctions in the IC31 group included categories of biofunctions that are directly linked to adaptive immune responses, such as humoral immune response and cell-mediated immune response. Furthermore, many of the other categories, such as cellular development and hematological system development, contained biofunctions associated with the initiation of adaptive immune responses, including proliferation of lymphocytes and quantity of lymphocytes, respectively. Pou class 2-associating factor 1 (Pou2af1), a molecule shown to play an important role in multiple stages of B cell development and GC formation in mice upon immunization with T cell-dependent antigen 20,21 , was one of the transcripts involved in most biofunctions of adaptive immunity (Fig. 5b). Other biofunctions associated with the activation of adaptive immune responses and lymphocyte differentiation were detected in several other modules in the IC31 group, particularly the turquoise module ( Fig. 4e and   For CAF01, biofunctions associated with lymphocyte activation were most abundant in the turquoise module, with the majority related to T cells and antigen-presenting cells (APCs) across different categories, such as differentiation of myeloid dendritic cells (DCs), differentiation of CD4 + T lymphocytes and class switching of B lymphocytes (Fig. 4c). Transcripts for the cytokine Tgfb and the transcription factor Batf were among the co-expressed molecules that were most commonly associated with the differentiation of APCs and lymphocytes in the CAF01 group (Fig. 5a). Batf is a transcription factor with important roles in the development of Th17, Th2 and TFH cells, as well as in the differentiation of Ab-producing B cells 22 .
For GLA-SE, most of the highly significant biofunctions associated with lymphocyte activation were detected in the yellow module. These biofunctions included differentiation of T lymphocytes, function of B lymphocytes and quantity of TFH cells ( Fig. 5c and Supplementary Table 1). The biofunctions associated with lymphocyte activation and differentiation in the yellow module of GLA-SE in WB included transcripts of nuclear factor of activated T cells 1 and 2 (Nfatc1 and Nfatc2) and IL-21 receptor (Il21r) (Fig. 5c), which have known functions in T cell activation, TFH activation and GC reactions.
Overall, co-expression analysis of DEGs in WB revealed that CAF01, IC31 and GLA-SE had distinct gene modules associated with biofunctions of adaptive immune responses, and that these gene modules contain molecules with documented roles in mounting characteristic T cell responses induced by each adjuvant.
Co-expression and biofunction analysis of dLNs. The co-expression analysis of DEGs in dLNs identified distinct modules for IC31 and GLA-SE. However, due to the limited numbers of significant DEGs in dLNs in the Alum and CAF01 groups, the co-expression analysis for these adjuvants revealed no distinct modules. Two modules were defined for IC31 and 11 modules for GLA-SE, with variable numbers of transcripts ranging from 70 DEGs in the blue module for IC31 to 774 DEGs in the turquoise module for GLA-SE (Supplementary Table 2). The small blue module of IC31 contained molecules involved in immune cell trafficking, such as the accumulation of granulocytes and recruitment of leukocytes, as well as the initiation of adaptive (especially cell-mediated) immune responses, including the density and recruitment of T lymphocytes. The turquoise module of GLA-SE contained the most highly significant biofunctions, many of which were associated with infectious diseases, inflammatory responses and hematological system development and function. Other significantly changed biofunctions for GLA-SE associated with the initiation of immune responses following vaccination, included the cellular immune response, quantity of APCs, and quantity of IgG. A few other immune-related biofunctions were spread across the different modules of the GLA-SE group, including the quantity of IgG2a and development of  Adjuvant-mediated induction of TFH cell population in dLNs. The unsupervised co-expression analysis described above identified adjuvant-specific gene modules that included molecules with known functions in TFH cells and/or the formation of GCs. This finding prompted us to study the impact of these human adjuvants on GC reactions and TFH cell responses in dLNs following immunization with H56 and the adjuvants. Mice were s.c.-immunized once with H56 alone or in combination with any of the four adjuvants. During the first 10 days a 10-fold expansion of the total dLN cell numbers was observed in the GLA-SE adjuvant-treated mice compared to the H56 group alone, whereas a more modest expansion ranging from 2-to 5-fold was induced by the other three adjuvants over the same time period (Fig. 6a). This suggests a very rapid and broad immune induction in the local dLN by GLA-SE, which was also indicated by the massive number of DEGs in both dLNs and blood within the first 72 h after immunization (Fig. 1). The TFH cells in dLNs on days 3, 5, 7, 10, 12 and 14 after immunization (Fig. 6b) were studied by evaluating the frequencies of CXCR5 hi , PD-1 hi CD4 + T cells (Fig. 6c). The frequencies and numbers of TFH cells remained low in the H56-immunized animals throughout the study. In contrast, all of the tested adjuvants enhanced TFH cell frequencies, and the most increase in frequencies of TFH cells relative to CD4 + T cells compared with H56 alone were observed 7 days after immunization. Similarly, adjuvant induced increase in total number of TFH cells was observed on day 7 compared with H56 alone, followed by a dramatic decrease on day 10 for IC31 and day 12 for GLA-SE and CAF01 (Fig. 6b). The mean frequencies of TFH cells were comparable among the different adjuvant groups on day 7, although the intra-group variations ensured that only the Alum and IC31 groups showed a statistically significant increase in TFH frequency compared to H56 alone (Fig. 6d). Higher numbers of TFH cells were only detected in the GLA-SE group (P < 0.01), but with a trend toward higher numbers in the other three adjuvant groups also when compared with the H56 alone group, reflecting a difference in total cell numbers in dLNs among the different adjuvant groups (Fig. 6e). Fully differentiated GC TFH cells have been shown to express GL7 (CXCR5 + GL7 + ), in contrast to the TFH (CXCR5 + GL7 − ) cells located outside of the GC 23 . High proportion of TFH cells in the dLNs of the different adjuvant groups expressed GL7, indicating that most of the cells falling within our TFH gate were fully differentiated with the highest frequency of GL7 TFH cells observed in CAF01 and GLA-SE groups on day 7 (Fig. 6f).
Taken together, these results show that all of the tested adjuvants enhanced the frequencies of TFH cells compared with H56 alone. Although a big variation in the TFH induction kinetics was seen among the adjuvants, no significant difference in frequencies was observed. GLA-SE induced the highest absolute numbers of TFH cells within the sampling period reflecting a significantly larger expansion of the dLN after vaccination when compared to the other adjuvants. The major proportion of adjuvant-induced TFH cells was also shown to express GL7 following immunization, indicating that these cells had completely differentiated into TFH cells.

Clinically tested adjuvants induced different GC reactions. TFH cells induce and regulate GC for-
mation as well as B cell maturation and differentiation culminating in high-affinity Ab responses. We investigated the induction of GC B cells in dLNs on days 3, 5, 7, 10, 12 and 14 after immunization, i.e., the same time points used in the TFH cell experiment described above. Although the induction of TFH cells peaked 7 days after immunization for all the adjuvant groups, the differentiation of B cells into GC B cells peaked on day 10 in terms of both frequency and numbers (Fig. 7a). GC B cells were gated as GL7 + and CD95 + among the CD19 + B cells and represented a clearly distinct population on day 10 in all of the adjuvant groups (Fig. 7b). Significantly higher frequencies of GC B cells were observed in the Alum (P < 0.01) and GLA-SE (P < 0.05) groups compared with the group receiving H56 alone, whereas more variation was observed in the CAF01 and IC31 groups (Fig. 7c). The greatest increase in the total number of GC B cells was observed in the GLA-SE (P < 0.001) group compared with the H56 group (Fig. 7d). Furthermore, GLA-SE (P < 0.05) was the only adjuvant that enhanced the number of cells with the PC phenotype (CD138 + B220 int/low ) at the early time points (Fig. 7e). Next, we examined the formation of GCs in dLNs on day 10 after immunization using immunohistochemistry with antibodies against GL7 (blue), B220 (red) and CD4 (green) (Fig. 7g). Immunohistochemical staining revealed that although the number of GCs in dLNs did not differ significantly among the immunized groups, all of the mice immunized with adjuvant plus H56 showed an increased GC size compared with the H56 alone-immunized mice (Fig. 7f). Statistically significant differences were only observed for the CAF01 and GLA-SE groups (P < 0.05) (Fig. 7f). The induction of GC B cells and PCs strongly correlated with the induction of TFH cells (Fig. 6h).
Although GC formation lagged behind TFH cell induction in mouse dLNs following immunization with the H56 antigen and the adjuvants, these results reveal that the two processes were correlated. Interestingly, all of the tested adjuvants enhanced GC formation, as demonstrated by the GC size and GC B cell frequency. However, differences were observed in the adjuvants' ability to enhance total numbers of GC B cells and early PC phenotypes in the dLNs.

Discussion
Systems biology approaches have recently been employed to define early signatures of immune responses induced by different human vaccines and adjuvants (reviewed in refs 5 and 6). However, direct comparative analyses of different human adjuvants have been hampered by a lack of uniformity in experimental design and data analysis, combined with the unavailability of proprietary adjuvants. Herein, we report the early molecular signatures in mouse WB and dLNs of three clinically tested vaccine adjuvants, namely CAF01, IC31 and GLA-SE, along with Alum.
Wide variation in transcriptional changes in both magnitude and kinetics was observed among these adjuvants, with Alum and GLA-SE inducing the smallest and largest numbers of transcriptional changes, respectively, in both WB and dLNs. These results are in accordance with a previous report showing that immunization with Alum induced only limited transcriptional changes in WB and dLNs when compared with GLA-SE 24 . CAF01 and IC31, on the other hand, induced fewer transcriptional changes in WB and dLN compared to GLA-SE. CAF01 and IC31 were previously shown to form a depot at the injection site with only slightly increased activity in the dLN over a prolonged time period 25,26 , and therefore it is plausible that a higher magnitude of gene expression could be observed at their injection site. The massive transcriptional changes induced by GLA-SE both in WB and in dLNs correlated with a much higher expansion of the total cell number in the dLNs of the GLA-SE group compared with those of the other three adjuvants. The differences between the CAF01 and IC31 groups and the GLA-SE group in terms of the overall kinetics and magnitude of in vivo transcriptional changes presumably reflect the different cell subpopulations targeted by the adjuvants and the in vivo kinetics of the adjuvants. Recently, GLA-SE was found to induce neutrophils to engulf antigen at injection site and travel to subcapsular lymphatics as well as T cell areas of the dLNs within 24 hours where they produced high amounts of IFN-γ 27 . Further, GLA-SE was shown to dramatically enhance the number of activated CD11c + GR1 + myeloid DCs in mouse LNs within 24 h after immunization 24 whereas CAF01 and IC31 were found to target a minute number of vaccine-associated DCs in vivo over a period of 2 weeks in mice 25,26,28 . These observations could, at least in part, explain the limited transcriptional changes induced by IC31 and CAF01 in both WB and dLNs as compared to GLA-SE. Furthermore, the early onset of transcriptional changes induced by GLA-SE in the dLNs (6 h) is in line with a previous report 24 . Nevertheless, the massive cell infiltration and early transcriptional changes associated with inflammatory responses induced by GLA-SE returned to the baseline levels within 7 days after administration. It is noteworthy that all of these adjuvants have previously been reported to be well tolerated with acceptable safety profiles in human trials.
Co-expression analyses have been widely used to elucidate the molecular mechanisms of various immunological properties, including human blood transcriptome modules, following vaccination 29,30 . The co-expression analysis of DEGs in WB presented herein revealed both common and unique biofunctional signatures for the three novel adjuvants. Both CAF01 and GLA-SE induced modules with biofunctions that were related to innate immune responses, such as the recruitment of various innate cells, chemotaxis and phagocytosis. IC31, on the other hand, did not exhibit a distinct module of DEGs associated with innate immune responses, and the corresponding DEGs were scattered across different modules. This is in accord with the observation that IC31 was the only adjuvant among the three without highly upregulated DEGs belonging to GO terms for innate inflammatory responses. The observed limited pro-inflammatory transcript expression is in line with a previous report that IC31 mainly targets monocyte-derived DCs by increasing the level of IFN-β and limiting the production of the pro-inflammatory cytokines TNF-α and IL-6 31 .
Gene modules in WB that were mainly associated with adaptive immune-related biofunctions were identified for all four adjuvants. An extensive immunological comparison of the same four adjuvants was recently reported by Knudsen et al. revealing that the adjuvants induced different CD4+ T cell polarization as well as IgG Ab subclass switching independent of the choice of antigen 7 . Our analyses of the main molecules involved in biofunctions related to adaptive immunity revealed intriguing differences that might explain the observed differences in B and T cell responses induced by the different adjuvants.
Interestingly, our co-expression analysis identified adjuvant-specific WB gene modules containing genes with known function in GC TFH cell-mediated B cell responses, including Nfatc1, Nfatc2 and IL21R (induced by GLA-SE) 32-34 , Batf (induced by CAF01) 22 and Pou2af1 (induced by IC31) 20,21 . All four adjuvants tested herein One representative section of a dLN showing immunohistochemical staining for GL7 (blue), B220 (red), and CD4 (green) on day 10 after immunization (g). The number of GCs per section was determined by microscopy, and the surface (μ m 2 ) of the GC was measured using ImageJ 1.49v software. Correlations between the numbers of TFH cells and GC B cells on day 10 for all groups are shown on the left, and the numbers of TFH cells and PC cells on day 10 are shown on the right (h). The results are expressed as the means ± SEM and significance was calculated using the Kruskal-Wallis test. Correlations were assessed using Spearman's rank test. *P < 0.05, **P < 0.01, ***P < 0.001. induced elevated TFH cell differentiation in dLNs 7 days after co-administration with H56, when compared to H56 alone. This was followed by induction of GC B cells peaking on day 10. It should be noted that the small percentages of TFH and GC B cells detected in the dLN of the H56 alone group support our previous findings in which H56 alone was found immunogenic, albeit its immunogenicity was dramatically enhanced when combined with CAF01, GLA-SE and IC31 7 . GLA-SE induced the highest numbers of TFH cells and GC B cells at day 10, which agrees with a recent report on the ability of GLA-SE in induction of a strong Ab response already after one immunization, whereas IC31, CAF01 and Alum induced significant Ab responses only after two or three immunizations 7 . The squalene-based adjuvant, MF59, has been reported to reach the dLNs within 6 h following co-injection with Ag where it promoted retention of intact Ag entrapped in immune complexes within macrophages translocating the vaccine to the follicular DCs pivotal in initiating the GC reaction 35,36 . The close resemblance between the SE part of GLA-SE and MF59 suggests that similar processes could be mediated by the two adjuvants. Thus, it is likely that the broad activation of dLN-resident cells in response to GLA-SE may play a role in the massive dLN expansion, the strong transcriptional changes and the early induction of GC B cells and PCs reported herein and the subsequent Ab response as previously reported 7 . Our data also support the notion that the number of TFH cells strongly correlates with the adjuvant-induced B cell responses 37 .
It has previously been reported that IC31 and CAF01 induce durable and potent characteristic CD4+ T cell responses in mice and humans 7,8,[38][39][40] . Interestingly, Batf and Tgfb, with documented roles in Th17 differentiation, were among the co-expressed molecules within adaptive immune-related biofunctions for the Th17-inducing adjuvant CAF01 41,42 . Further, genes involved in Th1 differentiation such as Pou2af1 and Btla, were among the co-expressed molecules for the Th1-tilting adjuvant IC31 43 . Nevertheless, we were unable to perform a statistical correlation analysis on the early molecular signatures and immunological readouts due to the lack of stratified immune responses in inbred mice.
In line with other publications, we observed subtle transcriptional changes induced by Alum compared with the other tested adjuvants, despite the ability of Alum in mounting a potent antibody response and the induction of TFH and GC B cell responses in the draining lymph nodes of mice 24,[44][45][46] . Mode of action of Alum is incompletely understood and suggested mechanisms of action include increased Ag uptake, recruitment of innate inflammatory cells, such as neutrophils, eosinophils, inflammatory monocytes, myeloid DCs and plasmacytoid DCs to the site of injection, and subsequent activation of innate immune responses through various danger signals 47,48 . Although our systems biology analysis of Alum induced response presented herein add to the growing body of evidence in the literature on the mode of action of Alum, further studies are needed.
In summary, the head-to-head comparative transcriptomics and systems biology analyses of clinical-grade adjuvants reported herein have revealed gene expression profiles, biological pathways, biofunctions and gene modules in mouse WB and dLNs that are shared by or exclusive to the different adjuvants. We could also identify characteristic blood gene modules that were highly enriched for adaptive immune-related biofunctions, including molecules with documented roles in TFH and GC responses. These observations were supported by immunological analyses showing that, although with different magnitude, all groups receiving adjuvanted H56 showed enhanced TFH and GC responses in dLNs, compared with the H56 antigen alone. These results provide new insights into the mechanisms of action of three clinically tested adjuvants, along with Alum, and may inform rational development of new adjuvants for human use.

Materials and Methods
Animals. Six-to eight-week-old female CB6F1/OlaHsd mice (Harlan Laboratories, The Netherlands) were housed in ventilated cages with free access to food and water. All animals were housed under standardized pathogen-free conditions at the Experimental Biomedicine Animal Facility, University of Gothenburg. Ethics Statement. The use of mice was performed in accordance with the regulations set forth by the Ethical Committee for Animal Experimentation in Gothenburg, Sweden and in accordance with European Community Directive 86/609. All the techniques and procedures were refined to provide for maximum comfort and minimal stress to mice. Experiments performed were approved by the Ethical Committee for Animal Experimentation in Gothenburg, Sweden under license 64/2015.

Adjuvants and antigen.
The effects of the adjuvants were evaluated with 5 μ g of clinical-grade Mycobacterium tuberculosis H56 fusion protein provided by Statens Serum Institute (SSI, Denmark). The following clinical-grade adjuvants were tested: Alhydrogel 2% (hereafter called Alum) provided by Brenntag Biosector (Denmark), CAF01 from SSI, IC31 from Valneva Austria GmbH (Austria), and GLA-SE from Infectious Disease Research Institute (Seattle, WA, USA). The adjuvants and antigens were combined strictly according to the manufacturers' recommendations. Alum was diluted to 3.4 mg/ml of aluminum in saline. Then, H56 in saline was added to Alum at a 1:1 ratio (v/v) and rotated at 4 °C overnight for adsorption. CAF01 was mixed with H56 in Tris buffer (10 mM, pH 7.4) at a 1:1 ratio (v/v) and gently vortexed. IC31 was mixed with H56 in phosphate-buffered saline (PBS) such that the final concentration of each dose was 100 nmol KLK/4 nmol ODN1a, and then the mixture was gently inverted several times to ensure a homogeneous suspension. GLA-SE was vortexed with H56 in PBS to obtain a final concentration of 5 μ g/dose. All vaccines were used within 2 h of final preparation.

Study design and organ collection.
Mice were lightly anesthetized with isoflurane (Baxter Healthcare) before immunization with H56 with or without adjuvant, as described above. Each mouse received a total volume of 100 μ l of the vaccine mixture injected s.c. on both sides of the base of the tail. Each group included 28 mice. A group of non-immunized mice (n = 7) were sacrificed and served as a naïve control group (time point 0 h). Seven mice from each group were sacrificed at each of the four time points (6 h, 24 h, 72 h and 168 h) after immunization. WB was collected from the axillary plexus in PAXgene RNA stabilizer (Qiagen GmbH, Germany) such that the final ratio of stabilizer to blood was maintained at 2.76, as recommended by the manufacturer. Immediately after bleeding, the mice were sacrificed by cervical dislocation and the inguinal LNs were harvested in 1 ml of RNAlater (Qiagen).
In another set of experiments, groups of mice were immunized with the vaccine formulations described above. dLNs were collected 3, 5, 7, 10, 12 and 14 days after primary immunization to evaluate effects of the different vaccine formulations and immunization schemes on TFH cells and humoral responses.

RNA extraction. RNA was extracted from dLNs and WB using the RNeasy Mini QIAcube kit (Qiagen) and
PAXgene kit (Qiagen), respectively, according to the manufacturer's instructions. Briefly, dLNs were homogenized in 350 μ l of RLT buffer with metal beads in a TissueLyser II for 2 min. Further homogenization was performed by adding the samples to QIAshredder columns, followed by centrifugation at 14,000 rpm for 2 min. The samples, DNase from an RNase-free DNase kit (Qiagen), 70% ethanol and all the necessary reagents from the appropriate RNA extraction kits were then added to the QIAcube, according to the manufacturer's protocol. The extracted RNA was kept on ice upon extraction, followed by an assessment of the RNA concentrations using a Nanodrop spectrophotometer (Thermo Scientific, USA). All the samples were examined on an Agilent 2200 TapeStation (Agilent Technologies, USA) to check the quality by measuring the RNA integrity number (RIN). High-quality RNA samples, which were defined as having a 260/280 ratio of approximately 2 and an RIN > 8, were used in downstream applications.
Whole-genome microarray analysis. High-quality RNA samples from WB and dLNs were subjected to a whole-genome microarray analysis on an Agilent platform (Agilent Technologies). The RNA was labeled with a Fluorescent Linear Amplification Kit according to manufacturer's instructions. The quantity and labeling efficiency were verified before the samples were hybridized to whole-genome 8 × 60 k mouse expression arrays, which were scanned at 5 μ m using an Agilent scanner. Image analysis was performed with Feature Extraction software (version 11.5.1.1, Agilent Technologies) to generate raw microarray data.
Microarray data acquisition and analysis. Microarray data were pre-processed and normalized using the limma package for the R suite. The raw data were first corrected for background using the "normexp" method. Background-corrected signals were normalized together using the quantile method. The normalized gene signals were further evaluated for differential gene expression using the moderate Student's t-test to compare the different groups with the control group at each time point, and the P values were further adjusted for multiple testing using the Benjamini-Hochberg method. Changes in expression induced by the adjuvants were calculated by dividing the log 2 expression value of each individual adjuvant-treated group by the log 2 expression value of the group receiving H56 alone at any given time point. A directional gene set GO enrichment analysis was performed using the PIANO package 49 . The results were plotted as a heatmap of the enrichment score (− log 10 (enrichment P value)). To identify global functional profiles, GO terms were traced back to their corresponding main GO hierarchy. GO relationships were retrieved from go-basic.obo (released November 10, 2014), and the different levels of the GOs were generated using a custom script. Differential gene expression (adj P value < 0.05) was further used to identify concerted responses using a gene co-expression analysis method called weighted gene co-expression analysis 50 . This analysis was performed in the R suite using the R package WGCNA 51 . IPA software (Qiagen) was used to evaluate biofunctions associated with gene lists derived from co-expression analyses. Visualization of the selected co-expressed modules was performed with Cytoscape software 52 . The complete set of the microarray data was deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE85339 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE85339).
Immunohistochemistry. Inguinal LNs (n = 3) were collected in Histocon (HistoLab) 10 days after immunization, embedded in Tissue-Tek OCT (Sakura), snap-frozen in isopentane cooled in liquid N 2 and stored at −80 °C until they were cut into 7 μ m sections. The sections were fixed in acetone and stored at −20 °C until staining. Three sections of different depths from each mouse were stained. Slides were blocked with 10% goat serum and stained overnight with biotinylated anti-B220 (RA3-6B2, BD Pharmingen), anti-CD4-FITC (L3T4, eBioscience) and anti-GL-7-eFluor660 (GL-7, eBioscience) at 4 °C, followed by a 1-h incubation with streptavidin-Alexa Fluor 594 (Life Technologies) at room temperature. Slides were mounted with Fluorescent Mounting Medium (Dako) and examined with a Zeiss Axio Imager.Z2 microscope. Images were acquired with Zeiss Zen Pro software. The number of GCs per section was counted under the microscope, and the GC area in each section was quantified using ImageJ 1.49 v software (NIH, USA). Assessment of cytokine and chemokines in serum. Serum samples were collected at 6 h and 24 h and stored at − 20 °C until analyzed with a Bio-Plex Pro Mouse Cytokine 23-plex Assay (Biorad). The luminex assay was performed according to the manufacturer's protocol and run on MagPix equipment (Invitrogen). Fold change (FC) calculation of proteins induced by the different adjuvants were calculated by dividing the mean values of each protein in different adjuvant groups by the mean values of those of the H56 alone group.
A Kruskal-Wallis post-test was used with a 95% confidence interval to compare multiple groups. Spearman's correlation analysis was performed, and P values < 0.05 were considered significant.