Single-cell analysis identifies conserved features of immune dysfunction in simulated microgravity and spaceflight

Abstract

Microgravity is associated with immunological dysfunction, though the mechanisms are poorly understood. Here, using single-cell analysis of human peripheral blood mononuclear cells (PBMCs) exposed to short term (25 hours) simulated microgravity, we characterize altered genes and pathways at basal and stimulated states with a Toll-like Receptor-7/8 agonist. We validate single-cell analysis by RNA sequencing and super-resolution microscopy, and against data from the Inspiration-4 (I4) mission, JAXA (Cell-Free Epigenome) mission, Twins study, and spleens from mice on the International Space Station. Overall, microgravity alters specific pathways for optimal immunity, including the cytoskeleton, interferon signaling, pyroptosis, temperature-shock, innate inflammation (e.g., Coronavirus pathogenesis pathway and IL-6 signaling), nuclear receptors, and sirtuin signaling. Microgravity directs monocyte inflammatory parameters, and impairs T cell and NK cell functionality. Using machine learning, we identify numerous compounds linking microgravity to immune cell transcription, and demonstrate that the flavonol, quercetin, can reverse most abnormal pathways. These results define immune cell alterations in microgravity, and provide opportunities for countermeasures to maintain normal immunity in space.

Fig. 1. Simulated microgravity without pattern recognition receptor stimulation alters the single-cell transcriptional landscape of human PBMCs.

A UMAP plot of unstimulated PBMCs single-cell transcriptomes (10X Genomics), pooled together from a male (36 years old) and a female (25 years old) donor, that underwent either 1G or simulated microgravity (uG) for 25 hours total. Cells were resolved into 28 distinct clusters. B Quantification of relative abundance of each cluster of single PBMCs by percentage, or log2Fold Change (FC) between simulated uG and 1G conditions. Source data are provided with this paper. C Volcano plot of differentially expressed genes (DEGs) across all immune cell types between uG and 1G; DEGs (including log2FC and adj. p) were calculated by the MAST method. The Benjamini–Hochberg (B-H) method was used for multiple comparison adjustments. Adjusted p value (adj. p) cutoff is 0.05, and log2FC cutoff is 0.25. D Dot plot showing the top DEGs from C and their expression levels across 22 immune cell populations. DEGs (including log2FC and adj. p) were calculated by the MAST method. P values were adjusted by the B-H method. Spot color reflects Log2FC of uG vs 1G, while spot size shows −log10 (adj. p). E UMAP of trajectory analysis of 1G and simulated uG unstimulated PBMCs. White circles represent the root nodes of the trajectory. Black circles indicate branch nodes, where cells can travel to a variety of outcomes. Light gray circles designate different trajectory outcomes. F Canonical pathway enrichment analysis obtained from Ingenuity Pathway Analysis (IPA) is shown across 19 immune cell clusters. Spot color reflects IPA z-score enrichment of simulated uG vs 1G, with red meaning predicted activation of the pathway in uG and blue meaning repression of the pathway in uG. Spot size shows the level of significance via −log10 (adj. p). Adj. p was calculated by the Fisher’s Exact Test (right-tailed) followed by B-H adjustment.

Fig. 2. Simulated microgravity induces predictive functional alterations in immune cells.

A Differences in iAge index between all cell types (left) and across 22 individual immune cell types (right) at 1G or simulated uG. Each box spans from the 25th to 75th percentiles (interquartile range, IQR), and features a median value denoted by a horizontal line. The whiskers extend to values within 1.5 times the IQR range from the 25th and 75th percentiles. 1G group (n = 11,934 cells examined over 2 independent experiments) is shown in blue and uG group (n = 16,568 cells examined over 2 independent experiments) is shown in yellow. The median for 1G group is 9.1, with min = −5.8 and max = 22.4. The median for uG group is 12.5, with min = −1.1 and max = 26.2. Two-tailed Mann–Whitney test (p value < 2.2e-16). ****p ≤ 0.0001, ***p ≤ 0.001, **p ≤ 0.01, *p ≤ 0.05. Source data are provided with this paper. B Differences in cellular senescence secretory product score, calculated from the SenMayo gene set, between all cell types at 1G or simulated uG. 1G group’s (n = 11,934 cells examined over 2 independent experiments) median is 0.001, with min = −0.084 and max = 0.302; uG group’s (n = 16,568 cells examined over 2 independent experiments) median is 0.006, with min = −0.091 and max = 0.569. Two-tailed Mann–Whitney test (p value < 2.2e-16). ****p ≤ 0.0001. Source data are provided with this paper. C Meta-transcriptome detection of mycobacteria, retrovirus, and total virus abundance in 1G (n = 11,934 cells examined over 2 independent experiments) and uG (n = 16,568 cells examined over 2 independent experiments) conditions. The bar plot shows the mean with an error bar representing the standard error of the mean (SEM). For Mycobacterium canettii (p value < 0.0001), mean of 1G = 2.5e-4 ± 9.0e-6 and uG = 5.5e-4 ± 1.4e-5. For Gammaretrovirus (p value < 0.0001), mean of 1G = 1.4e-6 ± 2.0e-7 and uG = 4.5e-6 ± 3.4e-7. For the total virus (p value = 0.0421), mean of 1G = 1.3e-5 ± 6.3e-7 and uG = 1.6e-5 ± 6.6e-7. Two-tailed Mann–Whitney test; *p ≤ 0.05, ****p ≤ 0.0001. Source data are provided with this paper. (D to F) NicheNet predicted significant ligand-receptor interaction between total T cells (Receiver) and the antigen-presenting cells (Sender) as D B cells, E DCs, and F monocytes in uG vs 1G condition (i.e., induced in uG over 1G).

Fig. 3. Simulated microgravity induces a distinct single-cell transcriptional landscape of human PBMCs following TLR7/8 stimulation.

A UMAP plot of TLR7/8 agonist stimulated (9 hours stimulation + 16 hours conditioning prior to stimulation = 25 hours total culture) PBMCs single-cell transcriptomics, pooled from a male (36 years old) and a female (25 years old) donor, that underwent either 1G or simulated uG. Cells were resolved into 23 distinct clusters. B Quantification of relative abundance of each cluster of single PBMCs by percentage, or log2FC between stimulated uG and 1G conditions. Source data are provided with this paper. C Volcano plot of DEGs across all immune cell types between TLR7/8 agonist simulated uG and 1G; DEGs (including log2FC and adj. p) were calculated by the MAST method. P values were adjusted by the B-H correction. Adj. p cutoff is 0.05, and log2FC cutoff is 0.25. D Dot plot showing the top DEGs from C and their expression levels across 19 immune cell populations. DEGs (including log2FC and adj. p) were calculated by the MAST method; p values were adjusted by the B-H correction. Spot color reflects log2FC of TLR7/8 agonist simulated uG vs 1G, while spot size shows −log10(adj. p). E UMAP of trajectory analysis of 1G and simulated uG TLR7/8 agonist stimulated PBMCs. Filled circle nomenclature was described in Fig. 1E. F Canonical pathway enrichment analysis obtained from IPA is shown across 19 immune cell clusters. Spot color reflects IPA z-score enrichment of TLR7/8 agonist-activated simulated uG vs 1G, with red meaning predicted activation of the pathway in simulated uG and blue meaning repression of the pathway in simulated uG. Spot size shows the level of significance via −log10 (adj. p). Adj. p was calculated by the Fisher’s Exact Test (right-tailed) followed by B-H adjustment.

Fig. 4. Simulated microgravity induces predictive functional alterations in immune cells following TLR7/8 stimulation.

A Differences in iAge index between all cell types (left) and across 19 individual immune cell types (right) after TLR7/8 agonist activation in 1G or simulated uG. Data are presented in the same way as described in Fig. 2A. 1G group (n = 14,916 cells examined over 2 independent experiments) is shown in blue and uG group (n = 12,230 cells examined over 2 independent experiments) is shown in yellow. The median for 1G group is 16.4, with min = 2.9 and max = 25.9. The median for uG is 13.2, with min = −1.7 and max = 23.9. Two-tailed Mann–Whitney test (p value < 2.2e-16). ****p ≤ 0.0001, ***p ≤ 0.001, **p ≤ 0.01, *p ≤ 0.05. Source data are provided with this paper. B Differences in cellular senescence secretory product score, calculated from the SenMayo gene set, between all cell types with TLR7/8 agonist activated 1G or simulated uG. 1G group’s (n = 14,916 cells examined over 2 independent experiments) median is −0.008, with min = −0.099 and max = 0.588; uG group’s (n = 12,230 cells examined over 2 independent experiments) median is 0.022, with min = −0.095 and max = 0.551. Two-tailed Mann–Whitney test (p value < 2.2e-16), ****p ≤ 0.0001. Source data are provided with this paper. C–E NicheNet predicted significant ligand-receptor interaction between total T cells (Receiver) and the antigen-presenting cells (Sender) as C B cells, D DCs, and E monocytes in TLR7/8 agonist activated simulated uG vs 1G condition (i.e., induced in uG over 1G).

Fig. 5. Validation of single-cell signatures identifies overlapping features of immune dysfunction in simulated microgravity and spaceflight.

A Volcano plot of DEGs from simulated uG vs.1G (25 hours) Bulk RNA-seq. Genes that are consistently upregulated across single-cell and bulk sequencing are labeled in red; genes that are consistently downregulated across the two datasets are labeled in blue. DEGs (including log2FC and p value) were calculated by the DESeq2 method; p value was determined by two-tailed Wald test and adjusted by the B-H method. Data were obtained from PBMCs from 3 male (ages 37, 22, 32 years old) and 3 female (age 27, 26, 40 years old) donors. B Spearman correlation of normalized counts between single-cell and bulk RNA-seq from simulated uG (R = 0.82, p < 2.2e-16) and 1G (R = 0.8, p < 2.2e-16) conditions. Two-tailed p value. C Venn diagram summarizing the overlapping DEGs between single-cell (SC; adj. p < 0.05, log2FC > |0.1|) and bulk RNA-seq (Bulk; p < 0.05) simulated uG vs. 1G. DEGs that are upregulated in both datasets are listed in the red box; DEGs that are downregulated in both datasets are listed in the blue box. The overlapping p value was calculated by the Fisher’s Exact Test, right-tailed. D Volcano plot of DEGs from Flight (ISS 33 days, n = 10) vs. Ground (n = 10) mouse spleen bulk RNA-seq (GLDS-420). Genes that are consistently upregulated across single-cell human PBMCs and bulk mouse spleen RNA-seq are labeled in red; genes that are consistently downregulated across the two are labeled in blue. DEGs (including log2FC and adj. p) were calculated by the DESeq2 method; p value was determined by two-tailed Wald test and adjusted by the B-H method. E Venn diagram summarizing the overlapping DEGs between human PBMCs single-cell (SC; adj. p < 0.05, log2FC > |0.1|) simulated uG vs. 1G and the mouse orthologous DEGs from Flight vs. Ground spleen bulk RNA-seq (GLDS-420; p < 0.05). DEGs that are upregulated in both datasets are listed in the red box; DEGs that are downregulated in both datasets are listed in the blue box. The overlapping p value (5.4e-14) was calculated by the Fisher’s Exact Test, right-tailed. F Heatmap of overlapping DEGs between human PBMCs simulated uG vs 1G and the I4 mission (n = 4) post-flight (R + 1) vs preflight (L-44) dataset. Both datasets are single-cell RNA-seq with DEGs defined by adj. p value < 0.05 and log2FC > |0.1|. Genes that are consistently upregulated across single-cell human PBMCs and I4 datasets are labeled in dark red (left). Genes that are consistently downregulated across datasets are labeled in dark blue (right). Genes that significantly overlap, but show reversal in their expression directions are labeled in gray. The overlapping p value (1e-200) was calculated by the Fisher’s Exact Test, right-tailed. G Heatmap of IPA canonical pathways enriched from DEGs between human PBMCs SC (single-cell RNA-seq uG vs 1G) and I4 (n = 4, R + 1 vs L-44). Enriched pathways have B-H adjusted p values < 0.05 (−log10(adj. p)>1.3). Red indicates a predicted activation in pathways, whereas blue indicates a predicted inhibition in pathways.

Fig. 6. Simulated microgravity induces distinct modifications to immune cell cytoskeletal morphology and cytokine production.

A Super-resolution microscopy analysis of actin in 2D for cell area (left), intensity (middle), and texture as punctate over diffuse index (PDI, variance/mean, right) between 25 hours of 1G or simulated uG. Dots represent individual PBMCs (n = 159 cells for 1G, and n = 154 cells for uG) from 4 independent donors. Donors were male (25 years old), and females (35, 38, and 46 years old). One outlier for actin intensity and actin PDI from each condition is removed based on Grubbs’ test. Two-tailed Welch’s t test was used for all comparisons. ***p ≤ 0.001. Data are plotted as mean ± standard error of the mean (SEM) and source data are provided with this paper. B Representative super-resolution microscopy images (2D left, 3D right) of PBMCs from 1G and simulated uG (25 hours) from donor 1 (35 yr F) and donor 2 (25 yr M; 2 of total 4 donors from A are shown here). 3D images better highlight changes to overall cell shape and actin protrusions in simulated uG. Scale bar = 2 μm and 1 μm, respectively. C Sixteen-channel granularity spectrum measurement of PBMCs stimulated with TLR7/8 agonist (9 hours stimulation, 16 hours conditioning prior to stimulation) from 1G (pink line) and uG (brown line) minus the corresponding unstimulated cells (25 hours total culture). The effect of simulated microgravity on unstimulated granularity spectrum is plotted in gray. Asterisks compare pink vs brown lines only. P values generated from unpaired two-tailed t test. n = 3 donors tested from A, 35-year-old female sample was not used. **p ≤ 0.01, *p ≤ 0.05. Data are plotted as mean ± SEM, and source data are provided with this paper. D Super-resolution microscopy analysis of 3D actin surface area (left, 1G n = 179 cells and uG n = 194 cells) and actin spike length (right, 1G n = 162 cells and uG n = 165 cells) between 25 hours of 1G or simulated uG. Dots represent individual PBMCs from 4 independent donors. Donors were male (25 years old), and females (35, 38, and 46 years old). Two-tailed Welch’s t test was used to calculate p values. *p ≤ 0.05, **p ≤ 0.01. Data are plotted as mean ± SEM and source data are provided with this paper. E G-LISA levels of active GTP-bound Cdc42 in PBMCs either unstimulated (25 hours) or treated with TLR7/8 agonist (9 hours + 16 hours conditioning) from 1G and simulated uG. n = 7, donors were male (25 years old), and females (38, 46, 25, 27, 26, and 40 years old). Two-tailed paired parametric t test was used to calculate p values, *p ≤ 0.05, ***p ≤ 0.001. Data are plotted as mean ± SEM and source data are provided with this paper. F ELISA levels of secreted IFNs by PBMCs treated with TLR7/8 agonist (16 hours conditioning + 9 hours stimulation) from 1G and simulated uG. n = 9, donors were male (36 years old), and females (33, 25, 38, 46, 27, 25, 26, and 40 years old). Two-tailed paired parametric t test was used, *p ≤ 0.05. Data are plotted as mean ± SEM, and source data are provided with this paper. G ELISA level of secreted ILs by PBMCs exposed to 25 hours simulated uG and 1G. n = 10 for IL-8 and n = 11 for IL-6, donors were females (32, 25, 38, 46, 25, 27, 26, 40 years old) and males (36, 33, 26 years old); 38-year-old female sample was not used for IL-8. Two-tailed paired parametric t test was used, *p ≤ 0.05. Data are plotted as mean ± SEM and source data are provided with this paper.

Fig. 7. Reversing simulated microgravity effects on the immune system with quercetin.

A Pipeline of microgravity and gene interacting compounds from discovery to validation. B Heatmap of top 50 simulated uG altered gene to compound interaction candidates. Compounds are listed on the right, and the predicted interacting genes are listed at the bottom. The color indicates the STITCH confidence score for compound-gene interaction. C Quercetin (50 μM) reverses the core gene expression signatures in simulated uG (25 hours). Log2FC levels of 106 core DEGs from simulated uG vs. 1G are plotted side-by-side to quercetin-treated uG vs. 1G in the heatmap. Red indicates positive log2FC, and blue indicates negative log2FC. 70% of the genes are reversed after quercetin treatment. The scatter plot below shows a negative association (Pearson correlation R = −0.35, p < 0.001) between the log2FC levels of the 106 core genes from simulated uG vs. 1G and quercetin-treated uG vs. 1G. D Gene set enrichment analysis (GSEA) shows the reversal effect of quercetin on the 106 core DEGs plotted in C heatmap. Quercetin treatment inverts the enrichment score (ES) in the upregulated core genes (from 0.8 to −0.64) and increases the ES of the downregulated core genes (from −0.75 to −0.55). All p values are <0.0001. E IPA Canonical pathways altered by quercetin. Heatmap plots the comparison of quercetin-treated samples against non-treated bulk RNA-seq controls, single-cell TLR7/8 agonist stimulated and unstimulated samples, and GLDS-420 mouse spleens in space. Red indicates a predicted activation in pathways, whereas blue indicates a predicted inhibition in pathways. The datasets were clustered by quercetin treatment into 2 major groups via complete linkage hierarchical clustering method. F Quercetin reduces senescence and age-associated inflammatory gene outputs. Both SenMayo scores and iAge index are reduced in the quercetin-treated group with p value of 0.0746 and 0.0268 respectively. Compared with the untreated group, quercetin downregulates more senescence-related and age-associated inflammatory genes (from ↓48.8% to ↓63.2%). n = 6, donors were 3 males (age 37, 22, 32 years old) and 3 females (age 27, 26, 40 years old). Two-tailed paired t test, *p ≤ 0.05. G Quercetin (25 hours treatment) reduces ROS levels measured by 2’,7’-dichlorofluorescin diacetate (DCFDA) assay. n = 6 for 1G vs uG, donors were males (32, 37, and 38 years old) and females (34, 32, 37 years old). 34 yrs and 32 yrs female samples were not treated with quercetin, resulting in n = 4 for comparisons between 1G vs 1G+Quercetin and uG vs uG+ Quercetin. Two-tailed paired parametric t test was used, *p ≤ 0.05. Data are plotted as mean ± SEM and source data are provided with the paper.