In the March 2020 issue of Nature Immunology, Weisel, Shlomchik, and co-workers presented interesting data pioneering the use of flow-purified B cells from BCR knock-in mice to explore substrate utilization and metabolic features of B lymphocytes (1). These included naïve B cells - in some but not all comparisons - and a population of germinal center (GC)-phenotype B cells recovered from recipients with a monomorphic B cell population designed to avoid inclusion of other B cells into GC (1–3). Comparisons also involved B cells after T-independent activation in vivo (1). In light of the limits to using bio-informatic data to reach conclusions about biological systems, the new evidence about fatty acid oxidation (1) advances insights beyond gene expression profiles comparing naïve and GC B cells (4). However, the paper evoked a need to evaluate the conclusive statement that the “GCBC transcriptome is not commensurate with [….] hypoxia” and similar broad conclusions of the text. This claim seems connected to a view of the authors that RNA-Seq data with GCBC do not contain evidence of enrichment for genes encoding glycolytic enzymes, or that such increases relative to naïve B cells would necessarily show up in a metabolomics analysis with ¹³C-labeled glucose. These issues prompted examination of these and other data sets in GEO. The results of the analyses point to limits to the conclusions as stated in (1); they also raise a broader question about the system used for this work.

Several papers (4–6) have documented results from intravital labeling with imidazole compounds that covalently modify cellular constituents when the mitochondria of viable cells operate under reductive conditions due to intracellular hypoxia (7–9). Work under controlled conditions has shown that a meaningful signal above background is obtained only when the ambient pO₂ is below about 1-1.5%, levels sufficient to yield HIF stabilization (7–9). Indeed, direct evidence of increased HIF-1α has been presented for both GC B cells (4, 5) and their Tfh counterparts (10), suggesting that many GC light zones are hypoxic. Of note, any issue of HIF function requires understanding that BCR engagement and TLR stimulation cause sustained HIF-1α and HIF-2α stabilization, which presents a drawback to comparing activated versus GC B cells. It might formally be possible that duration of the hypoxia, BCR signaling, and HIF stabilization failed to yield changes in mRNA concentrations large enough for enough gene products to yield a “statistically significant” result in a gene set enrichment analysis (GSEA) algorithm. Indeed, using a gene signature derived with the human breast cancer-like cell line MCF7 (11), the authors’ analysis of their purified GC B cells suggested that neither hypoxia-related nor HIF-1 target genes were enriched under their conditions of experimentation (1). Our previous published work had used GSEA with a gene signature indicative of biologically significant hypoxia (12) and scored the result of 2-fold normalized enrichment as statistically significant (4). In contrast, application of this gene signature to the data sets of (1) yielded a balanced mix of increased versus decreased mRNA and was not “statistically significant”. This difference along with other disparities of the data prompted us to compare the informatic findings while also using additional benchmarks that could test each report for independent replication.

Datasets comparing RNA expression data from “naïve” (IgD+ or Follicular) versus immunization-induced GC B cells were mined from the GEO depositions of raw sequencing data, all of which were generated with Illumina Hi-Seq 2000 or 2500 instruments. For a bespoke pipeline, sequences were trimmed using the fastp FASTQ preprocessor for overrepresented sequences and to remove any sequence with a quality score <10 (13). Trimmed sequences were aligned to the mm10 mouse genome using the STAR sequence aligner of Dobin et al. (14). Aligned sequences were then quality-tested using Qualimap (15) software, and counted using featurecounts in Rsubread (16). To cross-check results and to use a commercial platform readily accessible to anyone seeking to re-analyze the data herein, the RNA-Seq platform within the suites of Basepair Technologies were used. Processing of the gene expression and PCA were cross-checked by using the commercial pipeline of Basepair Technologies as applied to the primary data for naïve and GC B cells of the papers cited as (1, 4, 17, 18).

For heatmap generation ( Figures 2A, B ) via the Basepair Technology platform, first the raw read counts generated using STAR aligner and featurecounts were normalized using the DESeq2 package. DESeq2 (19) performs an internal normalization where geometric mean is calculated for each gene across all samples. The counts for a gene in each sample are then divided by this mean. The median of these ratios in a sample is the size factor for that sample. After this, a Z-score normalization is performed on the normalized read counts across samples for each gene, so that Z-scores are computed on a gene-by-gene (row-by-row) basis by subtracting the mean and then dividing by the standard deviation. Computed Z scores were then used to plot heatmaps in which each row represents one gene. There were no substantive differences between different data sets in the quality [for naïve and GC B cell RNA respectively, 92% and 91-92%, 95% and 94-95%, 94% and 93%, and 87% and 91% for the papers referenced in the main text as (1, 4, 13, and 14, respectively) via the in-house pipeline uniformly applied to all samples. Comparable values, all above 90%, were generated by the Basepair Technologies pipeline and STAR alignment algorithms. Length scoring also was indistinguishable among the different datasets, two of which [the papers cited as (1) and (18) in the main text] were generated by paired-end sequencing. These steps were followed by multi-factor differential expression analysis in which deposited datasets as well as “naïve” vs. “GCB” expression data were compared using DESeq2 (20). Effect size shrinkage of differential expression data was normalized using the apeglm method published by Zhu et al. (21). Euclidean distances and Spearman correlation coefficients were derived by applying a standard “dist” function in R to the data described above. Gene set enrichment analyses were carried out using the GSEA program of the Broad Institute (22, 23) and gene sets derived from published literature, the Rat Genome database (“RGD”), and the Broad KEGG database. As a technical and data note, the informatics pipeline reported in (1) used voom instead of DESeq2, followed by rankSumTestWithCorrelation instead of GSEA. As expected from the technical literature (19, 24–27), the results from each pipeline did not materially differ when applied to the GEO-deposited data of (1). For Figure 2D display of count data in a manner that reduces dependence on the variance of the mean for low count data, blinded dispersion estimation using the variance stabilizing transform, VST (28), was conducted on the DESeq2 count data followed by generation of a heat map in which the count data are displayed rather than relative expression.

GC hypoxia was observed with several types of immunization, including with NP-carrier (ovalbumin) (6), but only one paper (4), by Cho, Boothby, et al., had RNA-Seq data for comparisons. However, contemporaneous (2016, 2017) papers with data deposited in GEO had replicate data on naïve and bulk GC B cells (i.e., unfractionated mixtures of all of the diverse types of GC B cell, Figure 1A ) in a similar time frame (7-10 d) after immunization with the same immunogen (SRBC) (17, 18). We analyzed the sequencer output data for all four papers (1, 4, 17, 18) using the same two parallel pipelines (one assembled in-house, detailed in a technical log appended to this Perspective; a second via Basepair Technologies, Inc for an independent framework).

Several salient observations emerged from these comparisons. Unsupervised clustering with Spearman correlation analyses ( Figure 1B ), Principal Components analyses ( Figure 1C ), and the Euclidean distances among the different types of samples, revealed that the results from (1) differ substantially from the data of (3), (17), and (18). First, the mRNA expression pattern of GC B cells generated after transferring large numbers of B cells [apparently, 10⁶ - (2) as cited in (1)] biased toward a single specificity BCR (B1-8ⁱ Vκ -/-) into recipient mice with a monomorphic B cell population specific for Ig (AM14-Tg × V_κ8R-KI BALB/c mice, i.e., specific for allotype-disparate IgG2a^b) differed substantially from the other samples. A measure of overall differences between naïve and GC B cells was less in (1) (mean Euclidean distances of 126.3 ± 0.69) than the difference between the GC B cells of the transferred B1-8i, Vκ-/- cells (1) and the polyclonal GCBC (4, 17, 18) (mean Euclidean distances of 170.3 ± 10.3; p<0.01) as well as those of non-transgenic naïve versus GC B cells (147.6 ± SEM of 3.3; p<0.05). Second, the GC B cells from mice with a normal pre-immune repertoire (4, 17, 18) were substantially more similar to one another in comparing among three independent analyses (correlation coefficients in the range of +0.41 to +0.92) as opposed to anti-correlation of the B1-8i-derived GCBC, (correlation coefficients of -0.05 to -0.32). Although each of the independent data sets with the non-transgenic C57Bl/6 (B6) system differed somewhat – internally and from each other - the RNA-Seq “fingerprints” of naïve B cells in the polyclonal system were quite similar. Thus, the data were replicable when comparing independent SRBC immunizations of mice with no restrictions of the repertoire or BCR. In contrast, even the naïve B1-8ⁱ Vκ-/- B cells were quite distinct from the clusters of polyclonal naïve B cells ( Figures 1B, C ). Relative to the data from (18) as well as (4), the changes were more modest in (1). The similarities and differences among independent data sets can also be parsed by heat maps of differentially expressed genes among the data sets, analyzing either the total set thereof or highlighting specific genes that are known to encode major determinants of B cell positioning, signal initiation during activation, of differentiation into GC B or plasma cells ( Figures 2A, B ). In practice, the collective data ( Figures 1B, C ; 2B ) indicate that the GCBC generated after adoptive transfers of B cells with a restricted repertoire into mice whose endogenous B cells will contribute little to the GC are qualitatively distinct from a polyclonal response that evolved from a polyclonal repertoire.

Turning to the question, are there hypoxia-related gene signatures in GC B cells when taking into account other work with a polyclonal repertoire, GSEA algorithms were applied using several different gene sets for hypoxia ( Figure 2C ). A further facet of what the analyses revealed is that with some hypoxia modules, even the RNA-Seq data of (1) show significant enrichment in GSEA ( Figure 2C ). These increases can also be appreciated in data with the actual counts when displaying all naïve and GC B cell data ( Figure 2D ). Moreover, RNA-Seq data with GC from NP-CGG-immunized mice with a normal pre-immune landscape (29) also show enrichment for a functional hypoxia signature ( Figure 2C , data from the Pernis lab), and divergence from the B1-8i data. Metabolism rather than hypoxia was the central point of (1), a part of which was indirect evidence about glycolytic rates - perhaps on the view that if HIF were stabilized, expression of glycolytic genes should be increased. Prior work had provided evidence that the “glycolysis”, “mitochondrial respiration/oxidation” and “FAO” transcriptomes were increased in GCBC when compared to naïve B cells (4). Although this analysis was not evident in the Figure presented in (1) using the B1-8i system, a GSEA comparing naïve and GCBC with the RNA-Seq data in (1) found significant increases in glycolysis-related gene expression ( Figure 2C ), which also applied to each of the other data sets. Indeed, expression of genes encoding enzymes along the glycolysis pathway increased in GC B cell data sets of (1) as well as the polyclonal B6 GC B cells when compared to naïve B cells ( Figure 2D ).

The data presented here underscore that what is measured for an overall population of GC B cells is influenced by the structure of the experiment and the inputs analyzed. Key findings are that the cells analyzed in (1) were distinct from and changed less from naïve to GC phenotype than what is found in reproducible data with B cells that derived from a normal polyclonal repertoire (1, 17, 18). B1-8i, Vκ-/- B cells started out with substantial differences in activation markers such as Nr4a1 and differentiation-related genes, i.e., Prdm1 and Irf4. Of note, the differences do not relate to the particular informatic tools but instead to the selection of inputs. Thus, the pipelines used here yield similar results to those of (1) if, and only if, the analyses are restricted to those data sets and only use the selections in that work. Several factors may be involved in the differences between GC B cells in (1) from those in analyses published several years previously (4, 17, 18). Akin to findings with T cells (30), the experiments in (1) may involve the intra-clonal competition shown to result from using hundreds of thousands of transferred antigen receptor-transgenic cells (31, 32) - as already documented for B1-8i (31). In addition, the metabolic environments of highly diverse GC may differ from the secondary follicles that ensued with the transfer system used in (1). These factors warrant investigation, but the differences emphasize why it is vital that stated conclusions be restricted by limitations of the approach – not least when the approach is quite artificial.

GC are quite heterogeneous (33, 34). This point suggests that caution is warranted from the outset, militating against framing conclusions as blanket generalizations to apply across all GC. In terms of hypoxia or HIF gene signatures, prior work reported (4) that as many as 20% of splenic GC had no signal of intravital hypoxia, a variegation observed in parallel by others (J. Jellusova, personal communication). An integrative possibility is that hypoxia and/or its influence are reduced in GC designed to minimize bystander B cell involvement and dominated by a single specificity, as in (1–3). One technical issue is that there is no validated universal ‘hypoxia’ or ‘HIF’ module for activated lymphocytes, let alone a unitary module that includes GC B cells. Hypoxia and HIF responses are protean, but three decades of papers address cell type-specific functionality of a transcription factor. Thus, some restraint may be warranted before concluding that a GC transcriptome is not commensurate with hypoxia or conflating HIF with a particular magnitude of increased expression of genes encoding the enzymes of glycolysis or glycolytic flux. The issues connected to these questions take on added currency in considering antibody diversification during the persistent hypoxemia of many patients with severe CoVID-2019 infection. Thus, a great challenge will be to assess if human GC (in which questions about hypoxia, HIF, and in situ metabolism) exhibit effects such as those uncovered in mice (4–6) - either in normoxemic people or during concurrent hypoxemia.

Of course, such studies may rely heavily on scRNA-Seq and informatics analyses – conflicting or not – cannot settle secondary issues such as glucose uptake and utilization. Translation efficiency, post-translational modifications, and complex but unknown aspects of nutrient supplies, substrate concentrations, and substrate competition along webs of interconnected pathways all are downstream from the mRNA in question. That simply means that direct rigorous biochemical assays of glucose uptake per cell [e.g., ³H-2-deoxyglucose (DG), since 2-NBDG is known not necessarily to correlate with glucose entry rates (35)] and of glucose oxidation are needed before stating or accepting broad conclusions about germinal center B cells based on a single or special case as in (1). In this respect, normalizing a measure such as uptake to cell size – which is distinct from suitable compensation for autofluorescence in flow cytometry - lacks rigor. Normalization to size is analogous to concluding that a 160 kg person eating 4400 kCal/d has the same energy intake as one who is 80 kg and eats 2200 kCal/d, or claiming that a truck does not consume more fuel than a compact sedan if the size-normalized fuel consumption is the same. With respect to conclusions based solely on informatics, it is proposed that a reasonable standard is to (a) compare new RNA-Seq data sets to independent entries in GEO, and (b) quantify the correlation to - or difference from - those that are in GEO and readily comparable (i.e., not micro-array to RNA-Seq). Of essential importance, some restraint in statement of conclusions, along with openness about limitations, is vital at a time when the societal landscape is roiled by consequences of how scientists frame their work or evidence.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

MB wrote the text and guided analyses along with getting presubmission input from three independent faculty/PIs. AR led the primary informatic analyses, prepared figures, drafted portions of text, and edited the overall text. SC, KS, and JL performed analyses and provided insights into the RNA-Seq data and gene sets in the papers meta-analyzed. VHH provided input and guidance to SHC and MB; SHC provided input and guidance to KS and MB. All authors contributed to the article and approved the submitted version.

Departmental Funds - Department of Pathology, Microbiology, & Immunology, Vanderbilt University Medical Center.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.