This study utilized pooled-level genetic data to conduct a bidirectional, two-sample MR analysis, aiming to test for the association between GM and HSCR. In order to reduce bias and obtain reliable results, we tried to satisfy the following three hypotheses when using MR analysis: (1) Correlation hypothesis: instrumental variables (IVs) are closely related to exposure; (2) Exclusivity hypothesis: IVs do not affect outcomes by other means; and (3) Independence hypothesis: IVs are not affected by confounding factors (Davies et al., 2018; Figure 1).

Genetic variation for GM comes from the largest genome-wide meta-analysis of GM composition published by Kurilshikov et al. (2021) and MiBioGen Consortium (n.d.). They collected genome-wide genotypes and 16S fecal microbiome data from 18,340 individuals across 24 cohorts based on the MiBioGen consortium. A total of 196 taxa have been included in the study, classified under 5 biological categories: phylum, class, order, family, and genus. It should be noted that participants of European origin were exclusively considered, and 15 GM taxa (unknown family or genus) without specific species names were excluded from the analysis. GWAS summary statistics for HSCR were downloaded from the NHGRI-EBI GWAS Catalog (Sollis et al., 2023) on 10 November 2023 for study GCST005289 (Fadista et al., 2018). This study used the phenotype of HSCR and combined the most extensive collection of sporadic Hirschsprung cases to date, including 586 cases and 5,620 controls of European ancestry, which combined data from Denmark, the U.S., Finland, and Sweden. In the process of data analysis, strict quality control was carried out on the samples (Fadista et al., 2018).

The following selection criteria were used to choose the IVs: (1) single-nucleotide polymorphisms (SNPs) associated with each genus at the locus-wide significance threshold (p < 1.0 × 10⁻⁵) were selected as potential IVs (Li et al., 2023); (2) 1000 Genomes project European samples data were used as the reference panel to calculate the linkage disequilibrium (LD) between the SNPs, and among those SNPs that had R² < 0.001 (clumping window size = 10,000 kb), only the SNPs with the lowest p-values were retained (Fan et al., 2023); (3) SNPs with minor allele frequency (MAF) ≤ 0.01 were removed; and (4) when palindromic SNPs existed, the forward strand alleles were inferred using allele frequency information. (5) In order to control bias caused by IVs, the strength of each IV was evaluated using the formula F = β²/SE² to calculate the F-statistic, where β represents the effect size of SNP on exposure, and SE represents the standard deviation of SNP on exposure (Burgess et al., 2011). If the corresponding F-statistic was >10, it was considered that there was no significant weak instrumental bias (Burgess et al., 2011). Power calculations were conducted with the online Power Calculator Tool (Burgess, 2014; Online sample size and power calculator for Mendelian randomization with a binary outcome, n.d.).

MR is a statistical analysis method used for causal inference. It uses genetic variation as an IV to deduce the causal relationship between exposure and outcome, which can effectively avoid the influence of confounding bias in traditional epidemiological studies (Emdin et al., 2017). In this study, we analyzed the potential causal relationship between GM and HSCR with two-sample MR and strictly reduced the effect of confounders. The MR analysis process is shown in Figure 2. In the process of analyzing causal correlation, we used a variety of methods to evaluate it, including inverse-variance weighted (IVW; Burgess et al., 2017), MR-Egger test (Bowden et al., 2015), weighted median (WM), weighted mode (WMO), simple mode, and MR-PRESSO. The IVW method uses a meta-analysis method to combine Wald estimates of each SNP to summarize GM’s overall estimate of HSCR. If there is no horizontal pleiotropy, the IVW results are highly effective. MR-Egger tests the existence of the intercept term and uses it to evaluate pleiotropy. If the intercept term is very close to zero, then the MR-Egger regression model is very close to IVW, but if the intercept term is very different from zero, it indicates that there may be horizontal pleiotropy between these IVs. The WM is a method of considering the weight of the data when calculating the median, which can be used when the weight of the data should be considered more. The WMO takes into account the weights of the data points and is able to give a suitable mode if there are duplicate values with different weights. The simple mode can be used when each data point in the data set is weighted equally. MR-PRESSO improves the accuracy of causal associations by detecting and removing outliers. In the present study, we used the IVW method as the primary analysis method, while other methods were used to evaluate the stability of our results.

Heterogeneity was examined using Cochran’s Q-test, and a p-value of <0.05 was considered heterogeneous. The horizontal pleiotropy of IVs was measured by the MR-Egger intercept term. Meanwhile, MR-PRESSO was used to analyze pleiotropy further and pick out outliers in IVs. In addition, to identify potential heterogeneous SNPs, the “leave-one-out” analysis was performed by omitting each instrumental SNP in turn, aiming to verify the reliability and stability of the causal association (Chen et al., 2022).

To further evaluate whether there is a causal relationship between intestinal flora and HSCR, we performed a reverse MR analysis. Bacteria with a significant causal association with positive MR were selected as outcome variables, and HSCR was used as exposure for reverse MR analysis. The analysis methods and SNP selection conditions were consistent with those of positive MR.

False discovery rate (FDR) correction was conducted by applying the q-value procedure, with a FDR q-value of <0.1 (Storey and Tibshirani, 2003). Genera of GM and HSCR were considered to have a suggestive association when p < 0.05 but q ≥ 0.1 and a significant association when p < q < 0.05.

To characterize the bacterial community taxonomically, we performed 16S rDNA amplicon sequencing. First, genomic DNA was extracted from feces using the QIAamp Fast DNA Stool Mini Kit. Illumina-compatible primers (338F, 5′-ACTCCTACGGGAGGCAGCAG-30; 806R, 5′-GGACTACHVGGGTWTCTAAT-30) were then used to amplify the 16S rDNA V3–V4 regions by conventional PCR (Tian et al., 2023). The amplification products were subjected to gel purification prior to quantification by QuantiFluor-ST. After reverse transcription, the cDNA library was processed by fragmentation, end repair, and A-tailing using the TruSeq DNA PCR-free sample preparation kit. Then, the pooled library was run and sequenced on the NovaSeq600. According to QIIM, the effective Tags were obtained after chimera picking and quality filtering. Clustering by Uparse v7.0.1001 was performed to generate a list of open reference operational taxonomic units (OTUs) with an identity of 97%. The taxonomic assignment was subsequently achieved with the SSU rRNA database (Tian et al., 2023).

A targeted mass spectrometry assay was developed to detect the short-chain fatty acid (SCFA) levels in feces. Briefly, acetic, propionic, butyric, isobutyric, valeric, isovaleric, and hexanoic acids were purchased from Sigma Aldrich, as standard solutions with 10 concentration gradients (0.02, 0.1, 0.5, 1, 2, 5, 10, 25, 50, and 100 μg/mL). The sample was mixed with 50 μL of 15% phosphoric acid, 10 μL of 75 μg/mL of isohexanoic acid (internal standard), and 140 μL of diethyl ether, homogenated for 1 min, then centrifuged at 15,000 g at 4°C for 10 min. Next, the supernatant was collected for the test. The Agilent HP-INNOWax column was used for split injection, and the sample volume was 1 μL, the split injection ratio was 10:1, the inlet temperature was 250°C, the ion source temperature was 230°C, the transmission line temperature was 250°C, and the quadrupole temperature was 150°C. The programmed temperature increased from 90°C to 120°C at 10°C/min, next to 150°C at 5°C/min, and finally, to 250°C at 25°C/min for 2 min. The helium carrier gas flow rate was 1.0 mL/min. Mass spectrometry conditions included an electron bombardment ionization (EI) source, SIM scanning mode, and electron energy of 70 eV. The samples and standard solutions of SCFAs were measured by MS. The area under the curve for standard SCFA solution normalized to the internal standard was used as a single point on that SCFA standard curve. A standard curve for each SCFA was established based on the concentration gradient and served to calculate the concentration of this SCFA in the samples. All the samples were measured repeatedly, at least six times (Tian et al., 2023).

R software (version 4.3.1) and GraphPad Prism (version 8.0.1) were used to conduct the statistical analysis. We performed MR by the TwoSampleMR (version 0.5.7), MR-PRESSO (version 1.0; Verbanck et al., 2018), q-value (version 2.34.0; Storey and Tibshirani, 2003), and ggplot2 (version 3.3.3) R packages. Normally distributed data were presented as means standard error of the mean (SEM), and non-normally distributed data were presented as medians (interquartile range). Depending on the distribution and variances, a two-tailed Student’s t-test was applied to compare differences between groups. Significance was set at a p-value of <0.05. All figure panels represent data independently repeated at least three times, yielding similar results.

According to the selection criteria for IVs, 2,699 SNPs (p < 1.0 × 10⁻⁵) were selected for 196 intestinal flora in our MR analysis process. Details of all SNPs are provided in Supplementary Table S1.

IVW estimates were employed as the primary MR assay, and seven intestinal flora significantly associated with HSCR were identified. As shown in Figure 3, Eggerthella (OR: 2.66, 95%CI: 1.23–5.74, p = 0.01) was associated with an increased risk of HSCR, while Peptococcus (OR: 0.37, 95%CI: 0.18–0.73, p = 0.004), Ruminococcus2 (OR: 0.32, 95%CI: 0.11–0.91, p = 0.03), Clostridiaceae1 (OR: 0.22, 95%CI: 0.06–0.78, p = 0.02), Mollicutes RF9 (OR: 0.27, 95%CI: 0.09–0.8, p = 0.02), Ruminococcaceae (OR: 0.16, 95%CI: 0.04–0.66, p = 0.01) and Paraprevotella (OR: 0.45, 95%CI: 0.21–0.98, p = 0.04) were associated with reduced risks of HSCR. We corrected all p-values with FDR, and the Q-values were all greater than 0.05 (Supplementary Table S2). The other four methods are consistent with the direction of the IVW beta value (Figure 4).

Among these seven causal associations, the F-statistics of the IVs ranged from 17.42 to 32.26, eliminating the bias of weak IVs (Supplementary Table S3). As shown in Table 1, Cochran’s Q-test showed no significant heterogeneity for these IVs (p > 0.05; Supplementary Table S4), the MR-Egger regression intercept analysis found no horizontal pleiotropy (p > 0.05; Supplementary Table S5), and the MR-PRESSO test showed no outliers (global test p > 0.05; Supplementary Table S6). Moreover, after correcting the outliers, there is still no pleiotropy, which is supplemented by the MR-Egger regression intercept analysis. Subsequently, by drawing scatter plots (Figure 4) and leave-one-out plots (Figure 5), we did not find significant outliers for all IVs.

Reverse MR analysis showed that HSCR was associated with an increased abundance of Eggerthella (Supplementary Table S7; Figure 4), while no significant causal association was found with the other six bacteria (Supplementary Table S7). Cochran’s Q-test showed no significant heterogeneity in the IVs of HSCR (p > 0.05; Supplementary Table S8). Meanwhile, the analysis results of MR-Egger and MR-PRESSO (p > 0.05; Supplementary Table S8) show that horizontal pleiotropy is insignificant.

To investigate the alteration of the significant intestinal flora in HSCR, 16S rDNA sequencing was employed. Compared with the control, the relative abundance of Clostridiaceae (P: 0.002), Ruminococcus (P: 0.0003), and Paraprevotella (P: 0.00003) was significantly decreased (p < 0.05), while Eggerthellaceae was increased without statistical significance (p < 0.1) in HSCR. There were no significant changes in Ruminococcaceae (P: 0.671) or Peptococcus (P: 0.185; Figure 6A; Supplementary Tables S9, S10).

According to MR and the difference analysis, Clostridiaceae and Ruminococcus, as protective flora, showed a significant decrease in HSCR, and were reported as the producers of SCFAs (Vojinovic et al., 2019; Wu et al., 2023). Thus, we performed metabolically targeted mass spectrometry for SCFAs in fecal samples. Compared with the control, acetic, propionic, and butyric acid levels were lower in the HSCR group, while caproic acid levels were higher (Figure 6B; Supplementary Tables S11, S12). We further analyzed the correlation between significant flora and SCFAs. In the HSCR group, Ruminococcaceae was negatively correlated with acetic acid (R: −0.998, P: 0.041), while Clostridiaceae was positively correlated with isobutyric acid (R: 1.000, P: 0.002), and Ruminococcus was positively correlated with acetic acid (R: 1.000, P: 0.001; Figure 6C; Supplementary Table S13). However, in the control group, only Paraprevotella was negatively associated with acetic acid (R: −1.000, P: 0.012; Figure 6D; Supplementary Table S14).

These results suggested that the protective flora alternation was closely associated with SCFA levels in the intestinal microenvironment of HSCR patients.