Metagenome and metabolome study on inhaled corticosteroids in asthma patients with side effects

2.1 Clinical

The clinical arm of this study involved a cohort of 24 participants (equal gender), stratified into three distinct groups: individuals diagnosed with asthma experiencing inhaled corticosteroids (ICS) related side effects (Asthma.YES), and patients who do not experience any side effects from treatment. These patients tolerate the therapy well, without any negative reactions or complications (Asthma.NO), and a control group of healthy (HEALTHY) individuals. There are eight subjects per group.

The participants in the study were diagnosed more than 1 year ago with moderate, partially controlled Non-Allergic Asthma (Not linked to allergies). They were prescribed inhaled corticosteroids, which were used for more than three weeks. The participants did not have any other allergic conditions, and they were not taking any other medications. The patients were using inhalers 25 mcg/125 mcg, each single actuation provides salmeterol xinafoate equivalent to 25 μg of salmeterol and 125 μg of fluticasone propionate, 2 puffs 2 times daily.

Inclusion Criteria:

1. Adults aged 18–65 years with a confirmed diagnosis of asthma.

2. Regular use of inhaled corticosteroids for a minimum of 3 months.

3. Ability to provide informed consent and comply with study procedures.

Exclusion Criteria:

1. Use of systemic corticosteroids within the past 3 months.

2. Recent respiratory infections (within the past 4 weeks).

3. Presence of other chronic respiratory diseases (e.g., COPD, bronchiectasis).

4. Pregnancy or lactation.

The side effects include:

1. Hoarseness of voice (dysphonia): This occurs due to the deposition of corticosteroids on the vocal cords, leading to local muscle atrophy and mucosal irritation.

2. Dry mouth and throat irritation: ICS can reduce salivary secretion and cause mucosal dryness, resulting in discomfort and increased thirst.

3. Tachycardia and palpitations: These are attributed to the systemic absorption of ICS and the β2-agonist components in combination inhalers, which can stimulate cardiac β2-adrenergic receptors.

4. Severe cough: This may result from local irritation of the airway mucosa due to the inhaled particles.

5. Rashes on the mucous membranes of the mouth: Prolonged ICS use can lead to local immunosuppression, making the oral mucosa susceptible to infections and inflammatory reactions, presenting as rashes.

Each participant in the asthma group taking ICS was examined by professional pulmonologists to control asthma symptoms, as well as for the correct technique of using inhaled medications. For control symptoms, there were clinical examinations by a pulmonologist with the collection of all complaints, and patients also used the asthma control test questionnaire [13], 14], as well as the introduction of self-monitoring of peak expiratory flows [15]. To control the correct use of ICS, patients demonstrated their skills to the pulmonologist according to the rules of inhalation using an aero chamber [16].

Alongside stool sample collection, baseline health metrics including height, body mass index (BMI), and blood pressure were systematically recorded for each participant to provide a comprehensive clinical profile. The collection protocol was meticulously followed, ensuring the uniformity and integrity of the collected samples. Post-collection, all biological materials were immediately flash-frozen at −82 °C to preserve microbial and molecular fidelity, as recommended for microbiome and metabolomic analyses [17].

The study design, including the sampling framework and participant recruitment strategy, was reviewed and received ethical approval from the Kyrgyz National Surgical Center named after Mamakeev Ethics Committee on 07 June 2022, protocol number N16. [Supplementary Material 4]. This endorsement ensured adherence to ethical standards and the safeguarding of participant welfare throughout the study. The informed consent was obtained from all subjects and/or their legal guardian(s). Further details are available in Supplementary Material 4.

2.2 Sample preparation

The integrity of biological samples is critical for the reliability of downstream meta-omics analyses. Consequently, all collected materials were frozen at −82 °C immediately upon collection. This stringent protocol for sample preservation and storage was integral to minimizing microbial and enzymatic activity, thereby preserving the native microbial and metabolomic profiles for subsequent analysis. Furthermore, the samples were transported under cryogenic conditions for overseas analysis, ensuring the stability and integrity of the samples were maintained throughout the logistic process.

Specifically, stool samples were collected following the standardized procedures outlined [17], 18], ensuring the preservation of microbial DNA integrity for subsequent analysis.

2.3 DNA. Whole genome sequencing (WGS)

The genomic library was prepared to adhere to the DNBseq technology guidelines, employing a paired end read length of 150 nucleotides (PE150).

The WGS data for the stool microbiome were initially processed using the nf-core/tax-profiler pipeline, incorporating a series of quality control steps with FastQC and adapter sequence handling with Fastp. BBDuk was utilized for filtering out low-complexity sequences and those of low quality. Raw data with adapter sequences or low-quality sequences was filtered. We went through a series of data processing to remove contamination and obtain valid data. This step was completed by SOAPnuke [19] with the following filter parameters: “ −n 0.01, −l 20, −q 0.4, −adaMis 3, −outQualSys 1, −minReadLen 150”.

Other steps include:

1. Filter adapter: if the sequencing read matches 50.0 % or more of the adapter sequence (maximum 3 base mismatches are allowed), remove the entire read.

2. Filter read length: if the length of the sequencing read is less than 150 bp, discard the entire read.

3. Remove N: if the N content in the sequencing read accounts for 1.0 % or more of the entire read, discard the entire read.

4. Filter low-quality data: if the bases with a quality value of less than 20 in the sequencing read account for 40.0 % or more of the entire read, discard the entire read.

5. Obtain clean reads: the output read quality value system is set to Phred+64, which is integral for ensuring accurate base calling in sequencing data [20], 21].

6. Subsequent processing involved the removal of adapter sequences, potential contaminants, and low-quality reads, aligning with the best practices for next-generation sequencing data analysis [21]. See Supplementary Material 1–3 for more technical details.

2.4 Metabolomics

The analytical setup encompassed a Low-Speed Cryogenic Universal Centrifuge (Centrifuge 5430, Eppendorf), a Vortex Mixer (QL-901, Kylin-Bell Lab Instruments Co., Ltd.), and an Ultra-pure Water Meter (Milli-Q Integral, Millipore Corporation, USA). Sample preparation was executed using a Refrigerated Vacuum Concentrator (Maxi Vacbeta, GENE COMPANY). Sample disaggregation utilized a Tissue Grinder (JXFSTPRP, ShanghaiXinNing, China). Reagents of analytical calibration, including Methanol (A454-4) and Acetonitrile (A998-4) obtained from Thermo Fisher Scientific, Ammonium Formate (17,843–250 G, Honeywell Fluka), and Formic Acid (50,144–50 ml, DIMKA) were employed in the protocol.

Metabolite extraction commenced with a gradual defrosting at 4 °C, followed by the weighing of 25 mg of sample. The sample was combined with 800 µL of a pre-chilled extraction solvent (methanol:acetonitrile:water, 2:2:1,v/v/v) and 10 µL of an internal standard within a 1.5 ml Eppendorf tube. Tissue disaggregation was achieved via grinding with steel balls (50 Hz, 5 min), succeeded by ultrasonication in a 4 °C water bath for 10 min and a subsequent 1-h stabilization at −20 °C. Centrifugation at 25,000 g for 15 min at 4 °C facilitated the collection of the supernatant, which, following vacuum concentration, was re-dissolved in a 600 µL mixture of acetonitrile:water (7:3, v/v). Post-vortex and centrifugation, the supernatant was allocated into analysis vials. A composite quality control (QC) sample was formulated by amalgamating 50 µL from each sample to ensure analytical reproducibility in the LC-MS procedure [22].

Metabolite separation and detection were facilitated by a Waters UPLC I-Class Plus system coupled with a Q Exactive high-resolution mass spectrometer (Thermo Fisher Scientific, USA). Chromatographic separation was conducted on a BEH Amide column (1.7 µm, 2.1 × 100 mm, Waters, USA) employing a dual-solvent system: Mobile phase A (95 % acetonitrile with 10 mM ammonium formate and 0.1 % formic acid) and Mobile phase B (50 % acetonitrile with 10 mM ammonium formate and 0.1 % formic acid), both pH-adjusted to 9 using ammonia. The column temperature was maintained at 30 °C, utilizing a gradient elution program. The flow rate was set at 0.35 mL/min, and the injection volume was 2 µL.

Mass spectrometric analysis was conducted using the Q Exactive system, configured for both primary and secondary mass spectrometry data acquisition. Operational parameters included a full scan range of 70–1,050 m/z, with resolution settings of 70,000 for MS and 17,500 for MS/MS. The automatic gain control (AGC) target for MS acquisitions was set to 3e6, with stepped normalized collision energy (15, 30, and 45 eV) applied for precursor fragmentation. Sheath and auxiliary gas flow rates, as well as electrospray ionization voltages, were optimized for both positive and negative ion modes, and temperatures for the capillary and auxiliary gas heater were regulated at 320 °C and 350 °C, respectively.

Offline mass spectrometry data were imported into Compound Discoverer 3.3 (Thermo Fisher Scientific, USA) for metabolite identification and quantification against established databases such as the BGI Metabolome Database (bmdb), mzCloud Database, and ChemSpider. The generated data matrix, comprising metabolite peak areas and identification results, was further processed using the Compound Discoverer 3.3 platform (https://mycompounddiscoverer.com/).

The metaX software platform was pivotal in preprocessing this data, enabling the creation of detailed compound profiles and intensity metrics required for downstream analysis. This step was critical in discriminating between human and bacterial metabolites.

The software’s algorithms employed Probabilistic Quotient Normalization (PQN) for enhancing sample comparability: for every feature the mean response is calculated across all QC samples. A reference vector is then generated. The median between the reference vector and every sample is computed obtaining a vector of coefficients related to each sample. Each sample is then divided by the median value of the vector of coefficients; this median value is different for each sample.

Quality Control-based Robust LOESS (QC-RLSC) was used for correcting batch effects. Furthermore, a cutoff was established to exclude metabolites with a Coefficient of Variation greater than 30 % in the quality control (QC) samples, ensuring data reliability [23], 24].

2.5 Bioinformatics

Taxonomic categorization was conducted using Kraken2, supplemented by Bracken for species-level read reassignment. The differential abundance of microbial species across the cohorts was assessed with the DESeq2 R package, providing a quantitative landscape of microbial presence, FDR cutoff α = 0.1 was used to adjust p-value by the Benjamini–Hochberg procedure. To further probe into the meta-omics data, Multi-Omics Factor Analysis (MOFA) and Advance Software (ASAR) were employed to perform to facilitate advanced integrative analysis [25], [26], [27], [28], [29], [30], [31]. R package MOFA2 was used to calculate MOFA factors, source code is available from the GitHub repository.

The innovative GSA FBA-PRCC [32] methodology was leveraged to dissect the influence of exo-metabolite production and consumption on bacterial growth dynamics. By integrating metagenomic and metabolomic datasets, the platform provided insights into the stability of microbial communities and the roles of specific microbial taxa and metabolites within these complex ecosystems [32].

3.1 Taxonomical analysis

Our analysis revealed alterations in the abundance of various bacterial strains, underlining the complex influence of corticosteroids on the microbiota that extends beyond specific strains [33]. Notably, we identified an increased abundance of the following pathogens in asthma patients with and without corticosteroid side effects, designated in Table 1 as ASTHMA.YES and ASTHMA.NO respectively. The differential abundance of species was quantified using the DeSEQ2 R package, which demonstrated that specific species are overrepresented in asthma patients on corticosteroids [34].

Table 1 shows bacterial species that are present in higher quantities in the gut microbiome of asthma patients when compared to healthy individuals. Megamonas funiformis reduces complement C3 level, glucocorticoids and aldosterone (in alveolar epithelium) as well. So, the bacteria may increase the side effects of these steroids via complement C3 [35].

Bacteroides stercorischanges in abundance. have been correlated with asthma [36]. Phascolarctobacterium faecium is noted for its SCFA production. SCFAs have been implicated in the modulation of systemic immune responses, which are relevant to asthma [37]. Limosilactobacillus mucosae is traditionally linked with healthy mucosal surfaces, the implications of its abundance in asthma patients are not fully understood and represent a novel area of research [37]. Bifidobacterium pseudocatenulatum and Bifidobacterium longum are known for their health-promoting effects, including the potential suppression of allergic reactions, their unexpected prevalence in the gut microbiome of asthma patients may indicate a nuanced interaction with the host immune system [38].

Details of microbial species (Table 2) that are found in greater abundance in healthy individuals compared to those with asthma, potentially contributing to gut homeostasis and overall health.

Research on this specific bacteriophage like Blohavirus americanus: is minimal; however, the role of phages in regulating the bacterial populations and the immune response is an area of active research. Phages are increasingly being studied for their therapeutic potential and their impact on the microbiome [39]. Catenibacterium sp. co_0103 and Catenibacterium mitsuokai are part of the gut microbiota [40]. Known as a sulfate-reducing bacterium, the genus Desulfovibrio has been implicated in various gut-related issues. However, its role in asthma and respiratory health has not been well established and is being explored [41]. As mentioned earlier, variations in Bacteroides spp. have been correlated with several health conditions, including asthma [42].

Table 3 highlights a range of microbial species whose numbers are notably higher in the gut microbiomes of asthma patients with no side effects in comparison to the patients with side effects. Bacteriophages may impact gut bacterial populations and influence the side effects [43]. Typically, environmental Chryseobacterium rhizoplanae and Polaribacter sejongensis might indicate a state of dysbiosis. Such changes in microbial composition could potentially be linked to the immunomodulatory effects of corticosteroids [44]. While Bifidobacterium spp. is generally recognized for their health-promoting roles, an increase may reflect a compensatory reaction to inflammation or indicate a disrupted microbiome due to corticosteroid treatment [45].

3.2 MOFA analysis

Advancements in statistical software now permit the simultaneous analysis of diverse omics data derived from the same sample set. For this study, we employed the MOFA2 [46] package to uncover latent factors that simultaneously describe metagenomics and metabolomics data.

Three distinct analyses were conducted: the first encompassed all 23,970 features detected by LC/MS pipeline, the second analysis was limited to only the 5,081 identified metabolites, excluding the 18,889 features that lacked molecular annotation, and the third analysis was limited to as little as 1,214 metabolites that have KEGG reference [47].

Table 4 summarizes the variance explained in MOFA model by both metabolomics and metagenomics data. It is interesting that despite a decrease in the metabolomics dataset by order of magnitude the proportion of variance explained by metabolites did not change much.

These results indicate that metabolomics data account for a larger proportion of the variance explained by the MOFA model compared to metagenomics data, suggesting that metabolites may have a more pronounced role in the observed differences between patients. However, the significant variance in metagenomics data, particularly in the identified metabolites, emphasizes the need for further research to understand the complex interactions between the microbiome and metabolome in patients experiencing side effects from asthma treatment. As anticipated, a reduction in the number of features leads to a decrease in the portion of variance explained by metabolomics and an increase in the importance of metagenomics. In both models, six latent factors were identified, as shown in Table 5 below. The number of factors was limited by the number of samples in the dataset.

We observed that Factor2 (Figure 1) is predominantly influenced by metagenomics in all three models, while the remaining factors are more representative of metabolomics data. It is also evident that Factor2 better discriminates between asthma and healthy data compared to other factors.

Figure 1:

MOFA Beeswarm factor plots show the distribution of samples along their factor values. The samples are grouped and colored according to the diagnosis (A), diagnosis and side effects (B), patient gender(C) and diagnosis with patient gender (D).

The analysis of the factor plot in Figure 1 shows that Factor 2, which contrary to other factors influenced more by metagenomics rather than metabolomics, better differentiates healthy and asthma patients. Moreover, that factor does not separate samples either on side effects (Figure 1B) or on patient gender (Figure 1C and D). On the other hand, Factor 1 assigns values according to a combination of diagnosis and gender (Figure 1D). It is interesting that the healthy female samples groups with asthma males, while asthma females samples groups separately. It is also clear that factors 5 and 6 are of little use because they are dominated by single healthy female and male samples correspondently. Factors 3 and 4, despite the explanation of a meaningful part of the metabolomic data variation in the model, do not show meaningful data separation along any of the annotation factors. That is why we won’t analyze factors 3 to 6 any further.

Table 6 shows how many features, which have high weight absolute value, mapped to KEGG pathways. Factor 2, despite explaining 5 % of the variability in metabolomics data, has no metabolic features with a weight above 0.1 in absolute value. However, if we select 20 features with the most negative weights two of them: N-acetylserotonin (KEGG ID C00978) Indole-3-acetaldehyde (KEGG ID C00637) will be mapped to the Trytophan pathway (KEGG ID map00380).

Factor 1 mainly deals with metabolomics data and demonstrates a good separation of samples along with diagnosis and gender annotations. Seven out of 13 metabolic features with weight below −0.5 are mapped to the Biosynthesis of unsaturated fatty acids (KEGG ID map01040): Nervonic acid (KEGG ID C08323), Arachidic acid (KEGG ID C06425), Erucic acid (KEGG ID C08316), Palmitic acid (KEGG ID C00249), Docosanoic acid (KEGG ID C08281), Eicosapentaenoate (KEGG ID C06428), Lignoceric acid (KEGG ID C08320). We discuss the implications of fatty acid biosynthesis and tryptophan metabolism together with metabolome analysis.

Like what we saw in Factor 2 for metabolomic features, we could observe in metagenomics features for Factors 1 and 3: that there are no metagenomic features with an absolute value of weight above 0.1. That corroborates the finding that Factor 2 describes the variation of metagenomics data while Factors 1, 3, and 4 describe the variation of metabolomics data with little influence across views.

3.3 Metabolome analysis

The figures below illustrate major findings in metabolomics. Figure 2 shows overall changes in metabolites between all three groups: healthy people, asthma patients with and without corticosteroid side effects side effects. Figure 3 shows which pathways are affected in comparison to data from patients with and without side effects. Figure 4 shows which pathways are affected in patients with side effects compared to healthy people. Figure 5 shows which pathways are affected in patients with no side effects compared to healthy people.

Figure 2:

Overall changes in metabolites. Differences between three groups of patients represented: Healthy people (healthy), asthma patients with corticosteroid side effects (asthma.YES) and asthma patients without side effects (asthma.NO).

Figure 3:

Major changes in pathways. Pairwise analysis of data from asthma patients with corticosteroid side effects vs. asthma patients without side effects.

Figure 4:

Major changes in pathways. Pairwise analysis of data from asthma patients with corticosteroid side effects vs. healthy people.

Figure 5:

Major changes in pathways. Pairwise analysis of data from asthma patients without corticosteroid side effects vs. healthy people.

3.4 FBA-PRCC analysis

The structure of the microbial community is controlled by the interactions between community members and the environment. Most such interactions are mediated by small molecules called exo-metabolites, which are produced and consumed by members of the community. We applied FBA-PRCC [32] to investigate which bacteria and metabolites belong to which category. For each species in AGORA2 [48] collections of whole-genome models of human symbiotic bacteria we have calculated sensitivity coefficients linking metabolite influx/outflux and bacterial growth. By filtering out only coefficients significant at the level of 10 % we build a sensitivity graph in which 7295 bacterial strains are linked to 843 small molecules. Each edge of the graph has a weight corresponding to the partial rank-correlation coefficient between bacterial growth and metabolite availability. A negative sign of weight indicates that growth is suppressed by an increase of metabolite concentration. The direction of the edge defines either the production or the consumption of the metabolite influenced the growth.

Similar types of graphs are well studied about the protein-protein interaction (PPI) networks, where they are used to identify putative multi-subunit complexes and protein-disease relationships by property propagation from the seed set of known proteins [49]. However, PPI is usually undirected, and their edge weights are positive. Santolini and Barabasi [50] extend that type of analysis to directed graphs and signed weights in their DYNAMO algorithm. We adopted DYNAMO algorithm as it was implemented in BioNAR package [51] to evaluate the stable distribution of species and metabolite abundances in the community. The propagation algorithm simulates the mutual directional influence of species and metabolites on each other. For this we used seed abundances, for example, setting 10 most abundant species to their sample counts keep the rest as zero, and apply the propagation algorithm until convergence. The presence of negative weights could mean their extinction from the community. To avoid negative weights, we modified the propagation algorithm in a way that once a node acquires a negative weight its weight is set to zero. That modification allows the simulation of full consumption of metabolites and extinction of species.

We did have neither a human model in our graph nor diet information available, so we could use metabolomics data as the predefined environment and checked the structure of the steady-state distribution of the species dictated by metabolomic conditions. To compare species distribution, we used the Top-Down Concordance Coefficient (TDCC) [52] because it considers the fact that variability in high-abundance species is smaller than in low-abundant species, and so the contribution of the high-abundance species to the similarity measure should be higher. The top-down concordance coefficient (TDCC) was used to compare species weights with metagenomic abundances. Relative abundance of mapped 179 metabolites used as a “seed” set for propagation. In the FBA-PRCC graph, four Bacteroides species are significantly dependent on tryptophan. Differential abundance data shows Bacteroides intestinalis underabundant in asthma patients (Figure 6) supports findings. Linoleic acid is for one species Roseburia nov and one of the Roseburia genus is overabundant in asthma (Figure 6).

Figure 6:

FBA-PRCC tryptophan confirmation.

The central node, which becomes red is L-glutamate (Figure 7). Its high value is associated with Asthma.NO patients, but the number of samples is not sufficient for the rigorous statistical analysis.

Figure 7:

Steady-state propagation values are distributed within Roseburia strains. Blue is low value; yellow is medium, red is high value.