Informed consent and ethics approval

This study was completely in accordance with appropriate ethics guidelines. All participants consented at an informed consent briefing at SpaceX (Hawthorne, California), and samples were collected and processed under the approval of the institutional review board at Weill Cornell Medicine, under Protocol 21-05023569. All crew members provided written informed consent for data and sample sharing.

We sequenced analysed samples from human skin, oral and nasal environmental swabs before, during and after a 3-day mission to space. This dataset comprised paired metagenomic and metatranscriptomic sequencing for each swab. A total of 750 samples were collected in this study by the four crew members of the SpaceX Inspiration4 mission. The samples were taken from 10 body sites (Fig. 1a) across 8 collection points (3 pre-launch, 2 mid-flight and 3 post-flight) between June 2021 and December 2021. The crew additionally collected 20 samples from multiple Dragon capsules from 10 different locations. We note that some crew members (two adult male, two adult female) were using wet wipes (UPC, 036000317985) to bathe themselves in-flight in between swabbing; however, not every crew member did so, and SpaceX did not require this to be a consistent protocol among the crew. Wet wipes used by the crew were neither reused nor shared, which should limit any influence of this confounding variable. No statistical methods were used to predetermine sample sizes but our sample sizes are greater than any previous publication in this field.

The crew were each provided sterile Isohelix Buccal Mini Swabs (Isohelix, MS-03) and 1.0ml dual-barcoded screw-top tubes (Thermo Scientific, 3741-WP1D-BR/1.0ml) prefilled with 400l of DNA/RNA Shield storage preservative (Zymo Research, R1100). Following sample collection, swabs were immediately transferred to the barcoded screw-top tubes and kept at room temperature for less than 4days before being stored at 4C until processing. Additional descriptions of the sample collection and sequencing methods are available in companion publications37

DNA, RNA and proteins were isolated from each sample using the QIAGEN AllPrep DNA/RNA/Protein kit (QIAGEN, 47054) according to manufacturer protocol, yet omitting steps one and two. To lyse biological material from each sample, 350l of each sample was transferred to a QIAGEN PowerBead tube with 0.1mm glass beads and secured to a Vortex-Genie 2 using an adapter (1300-V1-24) before being homogenized for 10min. Of the subsequent lysate, 350l was transferred to a spin-column before proceeding with the protocol. Concentrations of the isolated DNA, RNA and protein for each sample were measured by fluorometric quantitation using the Qubit 4 fluorometer (Thermo Fisher, Q33238) and a corresponding assay kit. The Qubit 1Xds DNA HS Assay kit was used for DNA concentration (Q33231) and the RNA HS Assay kit (Q32855) was used for RNA concentration.

For shotgun metagenomic sequencing, library preparation for Illumina NGS platforms was performed using the Illumina DNA FLEX Library Prep kit (20018705) with IDT for Illumina DNA/RNA US indexes (20060059). Following library preparation, quality control was assessed using a BioAnalyzer 2100 (Agilent, G2939BA) and the High Sensitivity DNA assay. All libraries were pooled and sequenced on an S4 flow cell of the Illumina NovaSeq 6000 Sequencing System with 2150-bp paired-end reads.

For metatranscriptomic sequencing, library preparation and sequencing were performed at Discovery Life Sciences (Huntsville, Alabama). The extracted RNA went through an initial purification and cleanup with DNase digestion using the Zymo Research RNA Clean & Concentrator Magbead kit (R1082) following the manufacturer-recommended protocol on the Beckman Coulter Biomek i5 liquid handler (B87583). Following cleanup, ribosomal RNA reduction for RNA-seq library reactions was performed using the New England Bioscience NEBnext rRNA Depletion kit (Human/Mouse/Rat) (E6310X), and libraries were prepared using the NEBnext Ultra II Directional RNA Library Prep kit (E7760X) with GSL 8.8 IDT Plate Set B indexes. Following library preparation, quality control was assessed using the Roche KAPA Library Quantification kit (KK4824). All libraries were pooled and sequenced on an S4 flow cell of the Illumina NovaSeq 6000 Sequencing System with 2150-bp paired-end reads.

For faecal collection, all participants were provided with DNA Genotek OMNIgene-GUT (OM-200) kits for gut microbiome DNA collection. Each participant was instructed to empty their bladder and collect a faecal sample free of urine and toilet water. From the faecal specimen, each participant used a sterile single-use spatula, provided by the OMNIgene-GUT kit, to collect the faeces and deposit it into the OMIgene-GUT tube. Once deposited and sealed, the user was instructed to shake the sealed tube for 30s to homogenize the sample and release the storage buffer. All samples from each timepoint were stored at room temperature for less than 3days before storing at 80C long term. Faecal samples collected using the OMNIgene-GUT kit are stable at room temperature (1525C) for up to 60days.

DNA was isolated from each sample using the QIAGEN PowerFecal Pro DNA kit (51804). OMNIgene-GUT tubes were thawed on ice (4C) and vortexed for 10s. Then, 400l of homogenized faeces was transferred into the QIAGEN PowerBead Pro tube with 0.1mm glass beads and secured to a Vortex-Genie 2 using an adapter (1300-V1-24) before being homogenized at maximum speed for 10min. The remainder of the protocol was completed as instructed by the manufacturer. The concentration of the isolated DNA was measured by fluorometric quantitation using the Qubit 4 fluorometer (Thermo Fisher, Q33238), and the Qubit 1Xds DNA Broad Range Assay kit was used for DNA concentration (Q33265).

For shotgun metagenomic sequencing, library preparation for Illumina NGS platforms was performed using the Illumina DNA FLEX Library Prep kit (20018705) with IDT for Illumina DNA/RNA US indexes (20060059). Following library preparation, quality control was assessed using a BioAnalyzer 2100 (Agilent, G2939BA) and the High Sensitivity DNA assay. All libraries were pooled and sequenced on the Illumina NextSeq 2000 Sequencing System with 2150-bp paired-end reads.

All metagenomic and metatranscriptomic samples underwent the same quality control pipeline before downstream analysis. Software used was run with the default settings unless otherwise specified. The majority of our quality control pipeline makes use of bbtools (v.38.92), starting with clumpify (parameters: optical=f, dupesubs=2,dedupe=t) to group reads, bbduk (parameters: qout=33 trd=t hdist=1 k=27 ktrim=r mink=8 overwrite=true trimq=10 qtrim=rl threads=10 minlength=51 maxns=1 minbasefrequency=0.05 ecco=f) to remove adapter contamination, and tadpole (parameters: mode=correct, ecc=t, ecco=t) to remove sequencing error38. Unmatched reads were removed using bbtools repair function. Alignment to the human genome with Bowtie2 v.2.2.3 (parameters: very-sensitive-local) was done to remove potentially human-contaminating reads39.

We assembled all samples with MetaSPAdes v.3.14.3 (assembler-only)40. Assembly quality was gauged using MetaQUAST (v.5.0.2)41. We binned contigs into bacterial metagenome-assembled genomes on a sample-by-sample basis using MetaBAT2 v.2.12.1 (parameters: minContig 1500)42. Depth files were generated with MetaBAT2s built-in jgi_summarize_bam_contig_depths function. Alignments used in the binning process were created with Bowtie2 v.2.2.3 (parameters: very-sensitive-local) and formatted into index bamfiles with samtools v.1.0.

Genome bin quality was checked using the lineage workflow of CheckM (v.1.2)43. Medium and high-quality bins were dereplicated using deRep v.3.2.2 (parameters: -p 15 -comp 50 -pa 0.9 -sa 0.95 -nc 0.30 -cm larger). The resulting database of non-redundant bins was formatted as an xtree database (parameters: xtree BUILD k 29 comp 2), and sample-by-sample alignments and relative abundances were completed with the same approach as before. Bins were assigned taxonomic annotations with GTDB-tK (v.2.1.1)44.

To identify putative viral contigs, we used CheckV (v.0.8.1)45. For downstream viral abundance quantification, we filtered for contigs annotated as medium quality, high quality or complete. This contig database was dereplicated using BLAST and clustered at the 99% identity threshold as described above using established and published approaches (https://github.com/snayfach/MGV/tree/master/ani_cluster)46. The non-redundant viral contigs were formatted as an xtree database (parameters: xtree BUILD k 29 comp 0), and sample-by-sample alignments and relative abundances were computed with the same approach as before, the only difference being the coverage cut-off used to filter out viral genomes, which was lowered to 1% total and 0.05% unique due to the fact that those in question came directly from the samples analysed.

We generated gene catalogues using an approach piloted in previous studies47,48,49. Bakta v.1.5.1 was used to call putative open reading frames (ORFs)50. The annotations reported in this study (for example, Fig. 5) derive directly from Bakta. We clustered predicted and translated ORFs (at 90% requisite overlap and 90% identity) into homology-based sequence clusters using MMseqs2 v.13.4511 (parameters: easy-cluster min-seq-id 0.9 -c 0.9)51. The resulting non-redundant gene catalogue and its annotations were used in the functional analysis. We computed the abundance of the representative consensus sequences selected by MMseqs2 by alignment of quality-controlled reads with Diamond (v.2.0.14)52. We computed the total number of hits and computed gene relative abundance by dividing the number of aligned reads to a given gene by its length and then by the total number of aligned reads across all genes in a sample.

To identify viral taxonomic abundance via short-read alignment, we mapped reads to a database of all complete, dereplicated (by BLAST at 99% sequence identity) GenBank viral genomes. We used the Xtree aligner for this method (see below); however, given the difficulty of assigning taxonomic ranks to viral species on the basis of alignment alone, we first benchmarked this process. We used Art53 to generate synthetic viral communities at random abundances from 100 random viruses from the GenBank database. We then aligned (with Xtree) back to these genomes, filtered for 1% total coverage and/or 0.5% unique coverage, and compared expected read mapping vs observed read mapping. We additionally computed true/false positive rates on the basis of the proportion of taxa identified that were present in the mock community (true positive) versus those that were not (false positive) versus those that were present but not identified (false negative). Overall, we identified optimal classification at the genus level, with >98% true positive rate (that is, 98/100 taxa identified) and low false positive/negative rates (for example, <10 taxa not present in the sample identified) (Extended Data Fig. 10a,b). Species-level classification had higher false negative rates (generally arising from multimapping reads to highly similar species) and a 6070% true positive rate. Genus-level classification also yielded a nearly perfect correlation (>0.99 on average) between expected and observed read mappings (Extended Data Fig. 10c). As a result, while we report analyses for every taxonomic rank in the supplement, in the main text we describe only genus-level viral analysis.

In total, we used and compared seven different short-read mapping methods (MetaPhlAn4/StrainPhlAn, Xtree, Kraken2/Bracken run with four different settings, and Phanta), which together utilize five different databases that span bacterial, viral and fungal life. In addition, we identified and computed the relative abundance of non-redundant genes as well as bacterial and viral metagenome-assembled genomes. Subsequent downstream regression analyses were run on each resultant abundance table at each taxonomic rank.

Unless otherwise stated, for the figures involving taxonomic data used in the main text of this paper, we used XTree (https://github.com/GabeAl/UTree) (parameters: redistribute). XTree is a recent update to Utree54 containing an optimized alignment approach and increased ease of use. In brief, it is a k-mer-based aligner (akin to Kraken2 (ref. 55) but faster and designed for larger databases) that uses capitalist read redistribution56 to pick the highest-likelihood mapping between a read and a given reference based on the overall support of all reads in a sample for said reference. It reports the total coverage of a given query genome, as well as total unique coverage, which refers to coverage of regions found in only one genome of an entire genome database. We computed beta diversity (BrayCurtis) metrics for taxonomic abundances using the vegan package in R57.

For bacterial alignments, we generated an Xtree k-mer database (parameters: BUILD k 29 comp 0) from the Genome Taxonomy Database representative species dataset (Release 207) and aligned both metagenomic and metatranscriptomic samples. We filtered bacterial genomes for those that had at least 0.5% coverage and/or 0.25% unique coverage. Relative abundance was calculated by dividing the total reads assigned to a given genome by the total number of reads assigned to all genomes in a given sample. We additionally ran MetaPhlAn4 (ref. 58) (default settings) as an alternative approach to bacterial taxonomic classification.

For viral GenBank alignments, we generated an Xtree database (parameters: BUILD k 17 comp 0) from all complete GenBank viral genomes. We first dereplicated these sequences with BLAST 99% identity threshold via published approaches (https://github.com/snayfach/MGV/tree/master/ani_cluster)46,59. We filtered for genomes with 1%/0.5% total/unique coverage. Relative abundance was calculated identically as with the bacterial samples. We additionally ran Phanta (default settings) as an alternative to this approach for viral classification60.

As another set of methods for measuring taxonomic sample composition, we used Kraken2 and bracken, both with the default settings, to call taxa and quantify their abundances, respectively55,61. We used the default kraken2 reference databases, which include all NCBI listed taxa (bacteria, fungal and viral genomes) in RefSeq as of September 2022. We ran Kraken2 with four different settings: default (confidence=0) and unmasked reads, confidence=0 and masked reads, confidence=0.2 and unmasked reads, and confidence=0.2 and masked reads. In the cases where we masked reads before alignment (to filter repeats and determine whether fungal and other eukaryotic alignments were probably false positives), we used bbmask running default settings.

To evaluate our taxonomic profiling approach, we first compared the top ten genus-level classifications by body site before and after decontamination for each classifier in metagenomic and metatranscriptomic data. We observed general concordance among the various classification methods; for instance, the predominant skin genera consistently identified included Staphylococcus, Cutibacterium and Corynebacterium. The oral microbiome included Streptococcus, Rothia and Fusobacterium. Kraken2, which uses a database comprising both eukaryotic and prokaryotic organisms, identified fungi in the skin microbiome, as expected. The swabs from the Dragon capsule predominantly contained a diverse array of environmental microbes.

We compared these results at additional taxonomic ranks and with other taxonomic classifiers. For example, to discern higher specificity of the viral changes, we additionally fit species-level virus associations. While species-level viral taxonomic classification can be difficult due to high read misalignments (Extended Data Fig. 10), we wanted to determine whether we could observe a higher-resolution picture of viral activity due to spaceflight, as this effect is known to be space-associated (as opposed to bacterial skin to skin transmission, which could be a result of sharing tight quarters and not a space-specific effect).

We observed that many of the swabs collected, especially those from the skin sites, comprised low-biomass microbial communities; there are many documented challenges in analysing these data62,63. To filter environmental contamination and the kitome64 influencing our findings, we collected and sequenced negative controls of both (1) the water that sterile swabs were dipped in before use, as well as (2) the ambient air around the sites of sample collection and processing for sequencing.

Following taxonomic classification and identification of de novo assembled microbial genes, we removed potential contaminants from samples by comparison to our negative controls. We ran the same classification approaches for each negative control sample as described in the above paragraphs. This yielded, for every taxonomy classification approach and accompanying database, a dataframe of negative controls alongside a companion dataframe of experimental data. On each of these dataframe pairs, we then used the isContaminant function (parameters: method=prevalence, threshold=0.5) of the decontam package65 to mutually high-prevalence taxa between the negative controls and experimental samples. The guidance for implementation of the decontam package, including the parameter used, was derived from the following R vignette: https://benjjneb.github.io/decontam/vignettes/decontam_intro.html. Note that we used both metagenomic and metatranscriptomic negative control samples to decontaminate all data, regardless of whether those data were themselves metagenomic or metatranscriptomic. This decision was made to increase the overall conservatism of our approach.

Four mixed-model specifications were used for identifying microbial feature relationships with flight. Time is a variable encoded with three levels corresponding to the time of sampling relative to flight: pre-flight, mid-flight and post-flight. The reference group was the mid-flight timepoint, indicating that any regression coefficients had to be interpreted relative to flight (that is, a negative coefficient on the pre-launch timepoint implies that a feature was increased in-flight). We fit these models for all genes, viruses, and bacteria identified in our dataset by assembly, XTree (GTDB/GenBank), MetaPhlAn4, Kraken2 (all four algorithmic specifications), Phanta and gene catalogue construction. Each variable encoding a body site is binary, encoding whether a sample did or did not come from a particular region.

To search for features that were changed across the entire body, we fit overall associations, oral associations, skin associations and nasal associations:

$$begin{array}{l}{rm{ln}}left(rm{{microbial}{rm{_}}{feature}}{rm{_}}{abundance}+{minval}right)\sim {beta }_{0}+{beta }_{1}{rm{Time}}+left(1{rm{|}}rm{{Crew}.{ID}}right)+{epsilon }_{i}end{array}$$

(1)

For associations with oral changes, we used:

$$begin{array}{l}{ln}left(rm{{microbial}{rm{_}}{feature}{rm{_}}{abundance}+{minval}}right)\sim {beta }_{0}+{beta }_{1}{rm{Time}}times {rm{Oral}}+left(1{rm{|}}rm{{Crew}.{ID}}right)+{epsilon }_{i}end{array}$$

(2)

For associations with nasal changes, we used:

$$begin{array}{l}{rm{ln}}left(rm{{microbial}{rm{_}}{feature}{rm{_}}{abundance}+{minval}}right)\sim {beta }_{0}+{beta }_{1}{rm{Time}}times {rm{Nasal}}+left(1{rm{|}}rm{{Crew}.{ID}}right)+{epsilon }_{i}end{array}$$

(3)

For identifying associations with skin swabs, we fit the following model:

$$begin{array}{l}{rm{ln}}left(rm{{microbial}{rm{_}}{feature}{rm{_}}{abundance}+{minval}}right)\sim {beta }_{0}+{beta }_{1}{rm{Time}}times {rm{Armpit}}+{beta }_{2}{rm{Time}}times {rm{ToeWeb}}+{beta }_{3}{rm{Time}}times {rm{NapeOfNeck}}\+{beta }_{4}{rm{Time}}times {rm{Postauricular}}+{beta }_{5}{rm{Time}}times {rm{Forehead}}+{beta }_{6}{rm{Time}}times {rm{BellyButton}}\+{beta }_{7}{rm{Time}}times {rm{GlutealCrease}}+{beta }_{8}{rm{Time}}times {rm{TZone}}+left(1{rm{|}}rm{{Crew}.{ID}}right)+{epsilon }_{i}end{array}$$

(4)

The characters in each of the above equations refer to the beta coefficients on a given variable in that given regression. The characters refer to the regression residuals. Note that in the final equation (4), the reference groups are samples deriving from the nasal and oral microbiomes; this means that highlighted taxa will be those associated with time and skin sites as compared to the oral and nasal sites. We additionally fit these same model specifications without the random effect and compared the results in Extended Data Fig. 2. Data distributions were assumed to be normal but not tested for every single microbial feature. Individual data points for each feature are present in the online data stored at figshare66 and with NASA GeneLab (see Data availability).

We used the lme4 (ref. 67) package to compute associations between microbial features (that is, taxa or genes) abundance and time as a function of spaceflight and body site. For all data types, we aimed to remove potential contamination before running any associations. We estimated P values on all models with the ImerTest package using its default settings67,68. We adjusted for false positives using BenjaminiHochberg adjustment and used a q-value cut-off point of 0.05 to gauge significance.

We grouped microbial features associated with flight into six different categories. These were determined since our model contained a categorical variable encoding a samples timing relative to flight: whether it was taken before, during or afterwards. Since the modelling reference group was mid-flight, the interpretation of any coefficients would be directionally oriented relative to mid-flight microbial feature abundances. As a result, we were able to categorize features on the basis of the jointly considered direction of association and significance for the pre-flight and post-flight levels of this variable. The below listed categories are all included in the association summaries provided on figshare66 (see Data availability).

Transient increase in-flightnegative coefficient on the pre-flight variable level, negative coefficient on the post-flight variable, statistically significant for both

Transient increase in-flight (low priority)negative coefficient on the pre-flight variable level, negative coefficient on the post-flight variable, statistically significant for at least one of the two

Transient decrease in-flightpositive coefficient on the pre-flight variable level, positive coefficient on the post-flight variable level, statistically significant for both

Transient decrease in-flight (low priority)positive coefficient on the pre-flight variable level, positive coefficient on the post-flight variable level, statistically significant for at least one of the two

Potential persistent increasenegative coefficient on the pre-flight variable level, positive coefficient on the post-flight variable level, statistically significant for at least one of the two

Potential persistent decreasepositive coefficient on the pre-flight variable level, negative coefficient on the post-flight variable level, statistically significant for at least one of the two

We used these groups to surmise the time trends reported in the figures. It would be intractable to visualize every association of interest, so we prioritized within each category on the basis of the absolute value of beta-coefficients and adjusted P values. In Fig. 1c, we removed the low priority categories (two and four above) and only looked at the top 100 most increased and decreased significant genes, by group, relative to flight. We did so to make fitting splines feasible (especially in the case of genes, which had so many associations) and filter out additional noise due to low association-size findings.

We took a similar approach for the barplots in Figs. 24 and Extended Data Figs. 79. We again filtered out the low priority associations and selected, for each body site represented in the figure (for example, oral, skin, nasal), the top N with the greatest difference in absolute value of average L2FC relative to the mid-flight timepoints. In other words, we selected for microbial features with dramatic overall L2FCs. We maximized N on the basis of the available space in the figure in question. We note that the complete, categorized association results are available in the online data resource (see Data availability).

We modelled our species/strain-sharing analysis on the basis of ref. 30. Briefly, we used the s flag in MetaPhlAn4 to generate sam files that could be fed into StrainPhlAn. We used the sample2markers.py script to generate consensus markers and extracted markers for each identified strain using extract_markers.py. We ran StrainPhlAn with the settings recommended in ref. 30 (markers_in_n_samples 1, samples_with_n_markers 10 mutation_rates phylophlan_mode accurate). We then used the tree distance files generated by StrainPhlAn to identify strain-sharing cut-offs on the basis of the prevalence of different strains (detailed tutorial: https://github.com/biobakery/MetaPhlAn/wiki/Strain-Sharing-Inference).

The single-cell sequencing approach and averaging of host genes to identify expression levels are documented in refs. 33,69. The resultant averaged expression levels across cell types were associated with microbial feature abundance/expression using lasso regression. We used the same log transformation approach as in the mixed effects modelling for the microbial features, and we centred and rescaled the immune expression data. In total, we computed one regression per immune cell type (N=8) per relevant microbial feature, with the independent variables being all human genes (N=30,601). We selected features on the basis of their grouping described above, picking only those that were increased transiently or persistently increased after flight. Due to the volume of gene-catalogue associations, we only analysed persistently increased genes. We report outcomes with non-zero coefficients in the text.

The GNU parallel package was used for multiprocessing on the Linux command line70. We additionally used a series of separate R packages for analysis and visualization67,68,71,72,73,74,75,76. Figures were compiled in Adobe Illustrator.

No statistical method was used to predetermine sample size; all possible samples from all crew members (N=4) were taken. No sequenced data were excluded from the analyses; however, samples were quality controlled before bioinformatic and statistical analysis to remove duplicated reads, trim adapters and low-quality bases, remove human contamination and remove potential microbial contamination (using negative controls). The experiments were not randomized. Data collection and analysis were not performed blind to the conditions of the experiments.

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Read more:

Longitudinal multi-omics analysis of host microbiome architecture and immune responses during short-term spaceflight - Nature.com