Introduction

Sickle cell disease (SCD) is a chronic, inherited hematologic disorder causing significant morbidity and mortality across the globe. SCD is characterized by pathologic red blood cell hemolysis leading to complications stemming from vaso-occlusion and vascular stress, ultimately resulting in end organ damage. In a 2018 study, it was estimated 1,950 children are born with SCD in the U.S. annually, and the life expectancy of a child with SCD is 54 years compared to 76 years for the age- and race-matched U.S. population ¹. While there has been a steady increase in life expectancy for SCD patients resulting from diagnostic screening, early intervention, and better evidence-based treatment guidelines, SCD hospitalization rates have increased ². Pain crises, the clinical manifestation of acute vaso-occlusion, were the number one reason for SCD-related hospital admissions from 2004-2012². While gene therapy will play an important role in SCD treatment, there are significant barriers to widespread utilization ³ and with SCD births on the rise globally ⁴ alternative avenues of treatment must be pursued. In addition, novel biomarkers are valuable to better stratify SCD patients for targeted early intervention and more precise treatment decisions.

The gut microbiome is a physiological system with complex interactions with SCD. It was shown in a murine model of SCD that microbial antigens transit from the gut to the bloodstream where they activate neutrophils to a disease-promoting phenotype, termed “aged”, or “aged-like”, neutrophils (hereafter, “ANs”), and that depletion of gut microbes using broad-spectrum antibiotics decreased the AN population and improved inflammation-induced organ damage ⁵. Human SCD patients were also found to have higher levels of the microbial antigen lipopolysac-charide in their blood and an increased AN percentage when compared to a control group of iron deficiency anemia patients with similar hemoglobin levels ⁶. These findings led to the proposal of gut-targeting antibiotics as a treatment for SCD, and small studies in human subjects have shown antibiotics reduce the fraction of ANs and lipopolysaccharide in SCD patients ^7–9. The clinical momentum, though, has outpaced investigation of the SCD gut microbiome itself.

Gut microbiome studies that seek to identify microbial biomarkers of disease or disease severity look for abundance differences in individual taxa between case and control cohorts, as well as differences in summary measures related to community diversity. To date, there have been four studies that probed bacterial composition of the SCD gut microbiome compared to controls to determine if there are characteristic features of the SCD gut microbiome. Lim et al. did not observe differences in overall diversity for 35 SCD patients compared to a sickle trait control cohort, though one genus of Bacteroidetes had lower relative abundance in the SCD cohort ¹⁰. Similarly, in a cohort of 32 pediatric SCD patients, Mohandas et al. did not observe a difference in alpha diversity or beta diversity compared to an immunocompromised cohort and a non-immunocompromised control cohort ¹¹. In contrast, Brim et al. observed patterns of abundance differences at all taxonomic levels below phylum when 14 SCD patients were compared to healthy controls ¹². Brim et al. also observed that SCD patients have a lower Firmicutes:Bacteroidetes (F:B) ratio, a gut microbiome community-level metric relevant to human health ¹³. Finally, in a cohort of Angolan SCD pediatric patients, Delgadinho et al. did not observe a difference in alpha diversity compared to healthy siblings but did observe abundance differences at all taxonomic ranks ¹⁴, though fewer in number and mostly different than observed by Brim et al.

These prior studies suffer from two limitations in characterizing the gut microbiome of SCD patients. First, all studies suffered from small cohort sizes, limiting their ability to detect statistical differences in microbial taxa and population metrics. Second, the studies utilized 16S sequencing for community profiling, which cannot resolve bacterial species/strain level variation, cannot identify other microbes such as viruses, and underperforms relative to metagenomic sequencing in identifying low abundance taxa ¹⁵. Additionally, these studies found different changes in the SCD microbiome, leading to differing conclusions about its relationship with the disease. A robust evaluation of SCD patients’ gut microbiomes, together with clinical data, is necessary to determine whether there exists a relationship between the microbial community and disease pathology.

Microbiome markers have been shown to stratify treatment response in other diseases ^16,17, providing di-rection for the development of clinically-meaningful gut microbiome biomarkers. Additionally, gut bacteriophage populations have been recognized for their association with chronic disease ¹⁸ and have been shown to shift life-cycle distribution at the population level in inflammatory bowel disease ¹⁹. Interrogation of the SCD microbiome, in conjunction with studying the effect of AN levels in SCD patients, could identify novel biomarkers of disease severity in SCD. Thus, we conducted a study of blood and stool samples from patients with SCD compared with age and race matched controls. Using whole community metagenomics, we investigated the gut microbiome on multiple axes, revealing complex yet informative patterns of microbiome change. By comparing gut microbiome signatures with blood AN and cytokine levels, we begin to characterize the pathological interplay between the gut microbiome and immune system activation in SCD.

Results

Sickle cell disease patient cohort

Patients and controls in the study are followed at a large hospital system in the United States. Cohorts were matched for age, sex, ethnicity, and race (Supplemental Table 1). For patients, clinical history (Supplemental Table 2) and treatment history (Supplemental Table 3) were aggregated from the clinical record. Additionally, laboratory measurements were collected from the clinical record by their most recent reading (Supplemental Table 4) or assayed as part of the study (Supplemental Table 5).

Gut microbiome changes in sickle cell disease

To study the SCD gut microbiome, we sequenced fecal samples from patients and controls. Sample sequencing reads were profiled for microbial taxa using a library of taxa-specific marker genes ²⁰. To determine whether there is a shift in the SCD gut microbiome at the community level, we used Shannon diversity to compare the species richness between patients and controls and found that SCD patients have less diverse gut microbiomes (Figure 1a). To identify microbes positively or negatively associated with SCD, we modeled taxa abundance as a function of age, sex, ethnicity, race, and SCD status using a multivariable generalized linear model (Supplemental Table 1). Overall, 25 taxa are significantly associated with SCD and 6 taxa are associated with age (Figure 1b). We noted that SCD samples had higher abundances of the phylum Bacteroidetes (B) and lower abundances of the phylum Firmicutes (F), the two taxonomic groups that make up the composite metric F:B ratio, and observe a reduction in the ratio in our study (Figure 1c), as had been observed previously ¹². Looking at human gut microbiome indicator species identified in a meta-analysis of human microbiome studies ²¹, we observed a decrease in health indicators (p=5.6e-9; Mann–Whitney–Wilcoxon test) and increase in disease indicators (p=4.4e-3; Mann–Whitney–Wilcoxon test) in SCD patients compared to controls (Figure 1d), pointing to changes in the SCD gut microbiome that are consistent with changes observed in other human diseases.

Figure 1:

We also performed a functional characterization of samples using gene-based pathway annotation ²⁰. We again modeled abundance with the same fixed effects but included an additional effect, the F:B ratio, to both capture community-level changes and identify SCD-specific changes that are independent of the high-level shift (Figure 1e). High F:B ratio had exclusively negative associations that were significant, with the most significant associations with pathways related to nucleotide biosynthesis. SCD status was positively associated with three pathways related to butyrate fermentation, fatty acid oxidation, and vitamin B6 biosynthesis, and negatively associated with a pathway related to anaerobic metabolism. The small number of pathways robustly associated with SCD in our heterogeneous patient population point to specific functional changes that shed light on conserved contributions of the gut microbiome to SCD pathology.

Gut microbiome markers and known SCD biomarkers

We investigated whether gut microbiome changes in SCD are related to the AN population by asking whether microbiome markers correlate with the fraction of ANs of the total neutrophil population in a patient’s blood (Supplemental Figure 1). For this analysis, we utilized a subset of the SCD cohort whose blood was assayed (n=57). Controls were not assayed for ANs and are excluded from the following analysis. We evaluated four groups of markers: (1) species positively associated with SCD-E. lenta, F. plautii, and C. bolteae; (2) species negatively associated with SCD-O. sp. 57 20, E. siraeum, textitA. shahii, B. intestinihominis, R. callidus, B. angulatum, C. catus, S. isoflavoniconvertens, M. smithii, O. sp CAG 241, and C. comes; (3) functions significantly different between patients and controls-succinate fermentation to butanoate, fatty acid beta oxidation VI mammalian peroxisome, and fatty acid beta oxidation II; and (4) system-level summary metrics shown to be differential between patients and controls-Shannon diversity, F:B ratio, health indicators, and disease indicators. A single marker, the bacterial species A. shahii, significantly correlated positively with the AN percent (ρ = 0.26, p=0.048).

Given the substantial changes in the gut microbiome of SCD patients, we next investigated whether gut microbiome features are associated with SCD clinical and molecular measures of disease severity. We utilized a molecular blood panel to profile SCD patients for levels of cytokines and chemokines (Supplemental Table 5). Additionally, we included a number of hemolysis-related clinical measurements (Supplemental Table 4).

We first evaluated the AN fraction given its known role in SCD pathology. Neutrophil activation assay measures had significant correlations with white blood cell measures collected from the health record (Supplemental Figure 2), providing confidence in comparing AN data with clinical data. ANs were significantly correlated with four inflammatory cytokines: interferon gamma, interleukin-1β, interleukin-10, and interleukin-17A, and two chemokines: monocyte chemoattractant protein-1 and macrophage inflammatory protein-1B (Figure 2, top). In contrast, most of the gut microbiome metrics and all of the community metrics did not correlate with inflammatory markers (Figure 2, bottom). Exceptions to this pattern are A. shahii and B. angulatum, both negatively associated with SCD and inversely correlated with multiple marker levels. Gut microbiome features did have significant correlations with clinical disease markers, including absolute reticulocyte count (ARC), total bilirubin, and indirect bilirubin.

Figure 2:

Gut microbiome viruses

Having observed limited correlation between bacterial gut microbiome markers and ANs, we interrogated the viral population in metagenomic samples. Assembled metagenomes were profiled for viral sequences using a hybrid approach that combines viral marker genes and sequence features to distinguish chromosomal, plasmid, and viral sequences ²². Overall, there was a decrease in the number of viral sequences in SCD compared to controls (Figure 3a), likely reflecting the loss in diversity observed in SCD samples given its correlation with sample Shannon diversity and health indicators (p=2.6e-5; Mann–Whitney–Wilcoxon test) (Figure 3c). Proviruses are lysogenic virus sequences integrated into host genomes and can be labeled in assembled viral genomes by identifying sequence features such as viral integrase proteins and integration sites ²³. Proviruses constitute a small fraction of the viral sequences present in samples, yet there was a significant enrichment of provirus sequences in SCD compared to controls (p=4.7e-4; Mann–Whitney–Wilcoxon test) (Figure 3b), which was robust to provirus or lysogenic viral sequence calling method (Supplemental Figure 2). While the provirus fraction was not correlated with bacterial community-level metrics Shannon diversity and F:B, it was inversely correlated with health indicators and positively correlated with disease indicators (Figure 3c), providing evidence that provirus enrichment mirrors changes in the bacterial species present in the sample.

Figure 3:

From assembled provirus sequences it is not possible to determine whether a provirus is active, meaning replicating independently of its host via a lytic lifecycle, or inactive. Using a method developed previously ²⁴, we analyzed the sequencing read distribution of provirus sequences and their flanking host sequences to determine if there was increased read coverage of proviruses, which would indicate independent replication. Across SCD samples, >90% of proviruses were inactive (Figure 3d), and the active provirus fraction distribution in SCD was not different from controls (p=3.208e-1; Mann–Whitney–Wilcoxon test). Finally, we performed the same correla-tion experiment for viral-based metrics as we performed for bacterial markers. While the provirus fraction did not significantly correlate with AN fraction (ρ = 0.02, p=0.91), the provirus fraction had a number of statistically significant positive correlations with molecular cytokines (Figure 3e), which was not true of non-integrated temperate phages (Supplemental Figure 4).

Integrated viruses

After finding that the provirus fraction in the gut microbiome correlated with the level of multiple inflammatory cytokines in SCD patients, we decided to further investigate the provirus sequences. A total of 4,962 proviruses were identified across the SCD patient gut microbiomes. Provirus lengths had a bimodal distribution with a median size of 23,775 base pairs (Supplemental Figure 5). All proviruses with predicted taxonomy (99.5%) were bacteriophages of the class Caudoviricetes (99.3%), meaning the proviral sequences are prophages.

Because the vast majority of phages were categorized as inactive in host bacterial genomes, we wondered whether these sequences would be similar between patients. We clustered all prophage sequences at 99% identity over at least 70% bi-directional coverage and found 395 clusters of sequences with median size of 2 (Figure 4a). In total, 1,249 (25%) sequences had a nearly identical sequence in another sample in the dataset. These clustered prophages were homologous to sequences in common gut commensals, however Bacteroides stood out with 113 clusters matching species in the genus (Supplemental Figure 6). The decreased F:B ratio in our patient population may also represent an increase in Bacteroides strains carrying these conserved prophages, although at the community level, the prophage fraction did not inversely correlate with the F:B ratio as would be expected if it were such strains driving the community shift. To determine whether clustered prophages differed from singletons, we compared the high-level functional content of sequences in each group (Figure 4b). Interestingly, clustered prophage sequences had lower fractions of genes with functions related to lysis, transcriptional regulation, and host takeover and a higher fraction of genes with unknown function compared to non-clustered prophages. A reduction in proteins that function in host lysis was also observed in active prophages compared to dormant prophages (Supplemental Figure 7).

Figure 4:

Discussion

The primary goal of our study was to ask whether there exist consistent, significant interactions between the gut microbiome and disease pathology in SCD. Such interactions, when observed across heterogeneous patients with varying clinical histories and burdens of disease, could lead to novel treatment modalities and diagnostics to improve patient care and monitoring in this disease of increasing importance around the world. With orthogonal computational and experimental approaches, we can begin to build a model of the interactions between gut microbiome features, immune system activity, and pathology in SCD.

Our study reveals a two-tiered dialogue between the gut ecosystem and systemic inflammation in SCD (Figure 5). Shotgun metagenomics of 98 patients showed a disease-wide contraction of bacterial diversity, a drop in the F:B ratio, and characteristic changes in pan-disease indicator taxa. These shifts parallel patterns seen across disparate chronic illnesses ^21,25, suggesting that SCD provokes a generic disease signature at the microbiome community level.

Figure 5:

Neutrophil biology and the microbiome appear to occupy distinct, complementary niches. The AN subset correlated with an overlapping cytokine panel but not with bacterial community metrics. Conversely, bacterial α-diversity and indicator taxa associated with indirect bilirubin and reticulocyte count—clinical proxies for hemolysis, yet show little relationship to the immune milieu. This dissociation implies that neutrophil activation and microbial changes represent parallel arms of SCD pathophysiology, each informative for a different disease facet.

Beyond bacteria, dormant prophages emerge as a previously unrecognized axis of SCD biology. We detect a 1.7-fold increase in the prophage fraction, yet >90% of these elements are predicted to be replicatively silent. Strikingly, the proportion of prophages, but not total phage population or lytic phage fraction, tracks with IL-1α, IL-10, IL-17A and IFN-γ levels. Immune signaling molecules are known to be elevated at baseline in SCD patients ^26–28, and 3/4 cytokines correlated with prophage level are observed to be lower in patients in vaso-occlusive crisis compared to SCD patients not in vaso-occlusive crisis, including IL-10, IL-17A, IFNγ ²⁶. Interestingly, the gut microbiome in SCD patients harbors prophage sequences that are highly conserved between unrelated individuals, with one such sequence present in over 50% of metagenomes. These conserved, shared prophages were an unexpected feature of SCD patient microbiomes, and stand in contrast to studies of lytic phages, that found non-integrating bacteriophage populations are highly individual-specific ²⁹. We found that clustered sequences had a higher fraction of genes coding for proteins involved in host takeover, gene expression manipulation, and lysis. We speculate that lysogens bearing ‘domesticated’ prophages maintain a survival advantage in the setting of a disease-associated gut microbiome. Enrichment of lysogens may modulate host immunity by altering bacterial surface antigens or metabolite flux. Experimental induction assays and longitudinal sampling will be needed to test whether prophage excision precedes cytokine spikes or merely registers host stress.

Functional signature investigation of the gut microbiome, enabled by metagenomics, can have therapeutic implications ³⁰. Among the functional module changes observed in SCD, three merit particular attention: (i) enrichment of nucleotide biosynthesis pathways, (ii) succinate to butanoate fermentation, a non-canonical route for microbial butyrate generation, and (iii) fatty-acid β-oxidation pathways, which typically require molecular oxygen and are not common in the anaerobic gut environment. The communal shift toward a low F:B ratio was tightly coupled to an expansion of nucleotide-biosynthesis modules, offering a narrow-spectrum antimicrobial target set that could re-balance the ecosystem more precisely than broad antibiotics trialed to date ^7–9. Looking at functional modules enriched in SCD independent of community shift, we see the succinate to butanoate butyrate production pathway; a departure from the primary gut butyrate acetyl-CoA pathway ³¹. Butyrate is the primary energy source for gut epithelial cells and has been shown to be anti-inflammatory and supportive of epithelial-barrier integrity ³². This shift in butyrate production pathway in the SCD gut may reflect the loss of members of taxonomic groups that produce butyrate via complex carbohydrate fermentation, such as C. catus, C. comes and E. siraeum ^31,33. Conversely, elevated β-oxidation genes suggest greater flux through short- and medium-chain fatty-acid catabolism, a shift that might deplete beneficial metabolites and warrants metabolic follow-up. Notably, the SCD enriched β-oxidation pathways II and IV are mitochondrial and peroxisomal, and therefore generally aerobic, pointing to the availability of molecular oxygen in the gut, a hallmark of gut barrier disruption ^34,35.

There are a number of recognized limitations for our study. First, while we collectively grouped SCD patients into a single cohort, the disease can vary in its presentation and clinical course, which may affect the gut microbiome. Our cross-sectional design precludes causal inference and we acknowledge that as a chronic disease, SCD presentation can be highly variable. Second, flow-cytometry assays for neutrophil activation markers were performed within 8 hours of venipuncture, and while we have correlative evidence that our results reflect accurate measures, we did not formally benchmark marker stability to this time window. Next, prophage activity was inferred from coverage ratios, as is common in the field, but with this method asynchronous induction events could evade detection. Finally, findings derive from a single U.S. center and may not generalize to global SCD populations with different diets or genotypes.

In summary, we propose that the gut microbiome bacterial shift mirrors hemolytic burden while ANs and dormant prophages mirror immunological activation tone. Longitudinal and interventional studies, particularly those integrating metatranscriptomics and targeted bacteriotherapy, are now warranted to unravel causality and translate these markers into clinical tools.

Methods

Study design and subjects

The study is a single-center outpatient cohort study of SCD patients and age- and race-matched controls. Study participants were recruited from the Children’s Hospital at Montefiore and Montefiore Medical Center (IRB NUM-BER: 2018-9080). Patients were enrolled with the following inclusion criteria: 1) diagnosis of sickle cell anemia (sickle cell disease SS genotype or sickle beta zero thalassemia, Sβ⁰-thal), or sickle cell trait in association with hereditary persistence of fetal hemoglobin (S-HPFH); 2) age ≥ 4; and 3) currently in a usual state of health. Exclusion criteria included: 1) malignancy; 2) inflammatory bowel disease or other gastrointestinal disorder; 3) immunocompromised due to additional underlying disease or immunosuppressive medication; 4) history of hematopoietic stem cell transplant; 5) history of C. difficile infection in the preceding 2 months; 6) receipt of chemotherapy in the preceding 2 months; 7) receipt of systemic antibiotic other than Pen VK in the preceding 2 months; 8) Hospitalization within the preceding 2 weeks; 9) intercurrent febrile illness or sickle cell related pain episode requiring opioids within the preceding 2 weeks. Controls were included with age ≥ 4 and excluded with the same criteria in addition to any systemic antibiotic use in the preceding 2 months.

In total, 101 SCD patients and 66 controls were enrolled in the study from 2018-2020. Of these participants, 98 SCD patients and 46 controls provided stool samples. One patient sample did not have enough material for metagenomic sequencing. Fecal samples were collected by study participants into DNA/RNA Shield Fecal Collection tube (Zymogen Research) and samples were split into 6 tubes and stored at -80°C. Blood was also drawn from study participants around the time of fecal collection.

Clinical and molecular data

Clinical data was collected for SCD patients based on chart review. Information included: (i) demographicsage (years), sex, ethnicity, race, hemoglobin genotype (SS, Sβ⁰, SS with high F, AS, AA, where ‘S’ indicates a sickle cell allele, ‘β⁰’ indicates a thalassemia allele coding for no wild type hemoglobin, ‘F’ indicates a fetal hemoglobin allele, and ‘A’ indicates a wild-type allele); (ii) patient history-hydroxyurea (on/off), folic acid (on/off), glutamine (on/off), stroke (yes/no), asthma (yes/no), allergies (yes/no), bacteremia (yes/no), silent infarct (yes/no), history of meningitis, osteomyelitis, or urinary tract infection (yes/no), number of acute chest syndrome past year, number of pain admissions past year, number of 30 day readmissions past year, number of emergency room visits past year, number of transfusions in the past year, number of exchanged transfusions in the past year; and (iii) clinical measurements-weight (kg), body mass index (BMI), urine microalbumin (MA/Creatinine; mg/gm), lactate dehy-drogenase (LDH; U/L), platelet (PLT; k/uL), white blood cells (WBC; k/uL), absolute neutrophil count (ANC; k/uL), hemoglobin (Hb; g/dL), sickle hemoglobin (HbS; %), fetal hemoglobin (HbF; %), wild type hemoglobin (HbA; %), absolute reticulocyte count (ARC; k/uL), total bilirubin (mg/dL), direct bilirubin (mg/dL), indirect bilirubin (mg/dL), alanine aminotransferase (ALT; U/L), Creatinine (mg/dL), iron studies: iron level (ug/dL), transferrin (ug/dL), saturation (%), total iron binding capacity (TIBC; ug/dL), ferritin (ng/dL). All yearly measures were annualized to 1 year if patient records are < 1 year in duration. Acute care in the past year is the sum of pain admissions and emergency department visits.

Venous blood in Acid Citrate Dextrose was collected from study participants. Inflammatory cytokine and chemokine levels were measured from serum using the MILLIPLEX Human Cytokine/Chemokine/Growth Factor Panel A kit (MilliporeSigma, HCYTA-60K). Luminex assay data were acquired on a Luminex MAGPIX and analyzed with the MILLIPLEX Analyst program (MilliporeSigma). The following cytokines and chemokines were profiled: G-CSF (pg/mL), IFNγ (pg/mL), IL-10 (pg/mL), IL-17A (pg/mL), IL-1α (pg/mL), IL-1beta (pg/mL), IL-6 (pg/mL), IL-8 (pg/mL), IP-10 (pg/mL), MCP-1 (pg/mL), MIP-1β (pg/mL), TNF-α (pg/mL). Cytokines and chemokines with < 50% of samples assayed at the lower limit of detection were excluded from further analysis.

Neutrophil activation biomarker assessment

Blood was also evaluated by flow cytometry for neutrophil adhesion and activation markers. Flow cytometry was performed using LSRII equipped with FACS Diva 8.0.1 software (BD Biosciences) and analyzed with FlowJo software (Tree Star). Neutrophils were identified by forward and side scatter characteristics combined with CD16b expression, and the AN subset evaluated by CD62L^loCXCR4^hi expression within the neutrophil population. Whole blood from the same sample used for flow cytometry assays was diluted 1:10 in PBS for complete blood count with differential counts on ADIVA 120 (Siemens Healthcare Diagnostics). Total WBC and absolute neutrophil counts along with percent AN from flow cytometry were used to calculate the absolute aged neutrophil count. Samples were excluded if analyses were not performed within 8 hours of blood draw.

DNA extraction and sequencing

Sample preparation and sequencing was done by the Molecular Microbiology Facility of the Integrated Genomics Operation at Memorial Sloan Kettering Cancer Center. DNA was extracted from samples using a custom phenol chloroform extraction optimized for fungal and bacterial isolation. Metagenomic sequencing was performed using next generation sequencing with the Illumina HiSeq platform using 2x150bp paired end reads. For some samples, multiple sequencing runs were performed to achieve targeted sequencing depth. Base calling was performed by platform software resulting in paired FASTQ files for each sample. Samples with multiple runs were concatenated to single paired files.

Metagenomic sequence profiling

For bacterial taxa and function analysis, the bioBakery 3 suite of metagenomic sequence tools was used ²⁰. First, raw reads were fist processed for quality and removal of human contamination using kneaddata (v0.10.0) with the following parameters -t 8 -p 12 –cat-final-output ²⁰. Bacterial taxa were profiled with MetaPhlAn3 (v3.0) with default parameters and –no map –nproc 12 and utilized the mpa v30 CHOCOPhlAn 201901 clade-specific marker gene database ²⁰. MetaPhlAn output for all samples were collapsed to a single table using the merge metaphlan tables.py utility script ²⁰. Function pathways were profiled with HUMAnN3 (v3.7) ²⁰. Paired read files were merged to a single table and HUMAnN was run with default parameters and –threads 96²⁰ and database versions mpa vOct22 CHOCOPhlAnSGB 202212 and uniref90 201901b full.dmnd. For comparison between samples, pathway profiles were normalized for read sequence depth using the script humann renorm table to copies per million with parameters –units cpm –update-snames; stratified profiles were split using the humann split stratified table script; and merged to a single table using the humann join tables script ₂₀.

For viral profiling, raw sequences were processed with the nextflow nf-core/mag (v2.4.0) metagenomics pipeline ^36,37for quality and contamination removal, assembly with MEGAHIT ³⁸, and gene calling with prodigal ³⁹. The pipeline was run with-profile singularity and parameters –skip binning –skip spades. Viral sequences were profiled from MEGAHIT assemblies using geNomad (v1.7.1) with the end-to-end command and default parameters ²² and using VIBRANT (v1.2.1) using the VIBRANT run.py script with -f nucl -t 8 -no plot parameters ⁴⁰. Both methods predict provirus sequences but only VIBRANT predicts lytic vs. lysogenic bacteriophage life cycle. To calculate the non-integrated lysogenic virus count in the VIBRANT data, the count of provirus sequences was subtracted from the count of lysogenic viruses per sample. PropagAte (v1.1.0) with default parameters was used to determine the activity of provirus sequences ²⁴from raw sequence reads, MEGAHIT assemblies, and provirus coordinates called by both geNomad and VIBRANT. All fraction metrics were calculated by dividing the count of interest by the total number of viruses called by the respective method (e.g. the provirus fraction for geNomad is the count of proviruses divided by the total count of viruses per sample).

SCD metagenome analysis

Multiple metrics were calculated for the bacterial population in each sample using its MetaPhlAn profile. Shannon diversity was calculated using the MetaPhlAn utility script calculate diversity.R script with parameters -d alpha –m shannon -s s. The F:B ration was calculated by dividing the abundance of ‘t Firmicutes’ by the abundance of ‘t Bacteroidetes. Health and disease associated taxa were collected from a gut microbiome meta-analysis ²¹, and the number of taxa in each set was counted based on the presence of the taxa at any abundance in the sample, as performed in ²¹.

To determine bacterial taxa and pathways differentially abundant in SCD patient samples compared to controls, we utilized a multivariate generalized linear model approach to model profiles as log-linear with MaAsLin 2 (R; v1.8) ⁴¹. For taxa, MaAslin2 was run with min prevalence = 0.1 and min abundance = 0.0001. For pathways, MaAslin2 was run with min prevalence = 0.1 and normalization = ‘NONE’. We included a number demographic variables as fixed effects, modeling

where scd is an indicator variable for the disease status of a sample and abundance refers to either taxa or pathways. Effect heat maps were generated from MaAslin2 significant result tables, and significant effects were reported for q<0.05.

Prophages investigation

Prophage sequences identified using geNomad were clustered at 99% sequence identity over >= 70% of bidirectional sequence coverage using mmseqs2 (v14.7e284) ⁴² with parameters –cluster-mode 0 –cov-mode 0 –c 0.7 –min-seq-id 0.99. Clusters were converted to fully connected networks using the python package networkx (v3.1) ⁴³and visualized using cytoscape (v3.10.0) ⁴⁴. VPF-PLM ⁴⁵was used for functional annotation of prophage sequences. Clustered prophage homology to bacterial species sequences was done using blastn (v2.16.0+) and the core nt database with parameters evalue 1e-10 and perc identity 90. The best hit bacterial species was assigned to each cluster where such a hit was identified.

Statistical analyses and visualizations

Mann–Whitney–Wilcoxon was used to compare distribution means between SCD patients and controls for Shannon diversity, F:B ratio, healthy indicator fraction, and disease indicator fraction. Values were removed in +/-3 standard deviations from mean of all study samples. Spearman’s ρ was used for all tested correlations with the exception of flow cytometry and electronic health record white blood cell measures which was tested with Pearson correlation coefficient. Significance of a correlation was measured using a permutation test with a two-sided alternative hypothesis. Mann–Whitney–Wilcoxon test was used for all tested association of categorical variables test with a two-sided alternative hypothesis. All methods were implemented with scipy stats package ⁴⁶: spearmanr, pearsonr, and mannwhitneyu with default parameters and permutation test with alternative=‘two-sided’, permutation type=‘pairings’. Significance was assigned for p < 0.05 for all comparisons. Visualizations were produced using python packages seaborn ⁴⁷and statannotations ⁴⁸. Pandas ⁴⁹and numpy ⁵⁰python packages were used for analysis.

Supplemental material

Supplemental Figure 1:

Supplemental Figure 2:

Supplemental Figure 3:

Supplemental Figure 4:

Supplemental Figure 5:

Supplemental Figure 6:

Supplemental Figure 7:

Supplemental Table 1:

Supplemental Table 2:

Supplemental Table 3:

Supplemental Table 4:

Supplemental Table 5:

Data availability

The dataset "Sickle cell disease patient gut microbiome study" is available in the repository NCBI and can be accessed via the following accession ID: PRJNA1320713. Patient level clinical, neutrophil, and immune data can be found in the Supplementary Material file.

Acknowledgements

The authors thank Marcel R. M. van den Brink and Jonathan U. Peled at the Molecular Microbiology Facility of Memorial Sloan Kettering Cancer Center for their generous assistance with microbiome sequencing.

Additional information

Ethics approval and consent to participate

This study was approved by the Institutional Review Board of the Albert Einstein College of Medicine and Montefiore Medical Center (IRB 2018-9080). All study participants were recruited from the outpatient clinics at the Children’s Hospital at Montefiore and Montefiore Medical Center. Written informed consent was obtained from all participants or their legal guardians prior to enrollment.

Funding

Support for this study was provided in part by the NIH grants R01HL069438 (L.K. and D.M), 1U01DE035632 (L.K.), P30 AI124414 (B.C.H.), and the Einstein Medical Scientist Training Program 1T32GM149364 (Z.N.F.) and the Price Family Foundation.

Authors’ contributions

Z.N.F. developed the hypotheses, analyzed the data, and prepared all figures. Z.N.F and L.K. wrote the manuscript text with D.M. and B.C.H contributing. D.M. and P.S.F. designed the study and secured IRB approval. J.D.L.S., K.I., and J.K. supported study participant enrollment, sample collection, and data aggregation. J.S.K. conducted the neutrophil assay. A.M.M. and A.C. conducted the immune profiling assays.