Introduction

Caesarean section (CS) is a prevalent surgery worldwide, and its rate has increased in recent decades [1]. Although CS can significantly reduce dystocia and stillbirths [2], Caesarean section scar diverticulum (CSD) affect about 19.4-88% of women receiving this operation [3]. CSD emerges from poor healing of the local uterine incision, forming a depression or cavity that connects with the uterine cavity [4]. Recently, CSD has drawn widespread attention because of its potential damage to subsequent fertility. For example, Gurol and colleagues found that the possibility of subsequent pregnancy decreases by an average of 10% after CS relative to a previous vaginal delivery [5]. A niche can reduce the chances of embryo implantation and increase the likelihood of spontaneous miscarriages if the implantation occurs close to or in the CSD [6]. Our previous studies have shown that the persistent effusion of CSD is a key cause of failed embryo implantation [7].

The unique microbial community composition in the female reproductive tract is essential in maintaining female reproductive health [8-10]. Our previous studies have demonstrated that the abnormal alterations in the local microbiota of Cesarean scar defects (CSD) cause continuous leakage through local inflammation and immune imbalance [11]. Further investigations of ours revealed the underlying mechanism demonstrating that abnormal bacteria in CSD deplete protective fatty acids and generate N-(3-hydroxy-eicosanoyl)-homoserine lactone, thus leading to promoting apoptosis of vascular endothelial cells and endometrial epithelial cells, ultimately impairing women’s reproductive capacity [12]. Nonetheless, the causative factors behind bacterial dysbiosis in complex microbial communities remain poorly understood.

Fungi and bacteria, as integral components of the human microbiome, establish complex interactions with each other and the host. The collection of genomes and genes carried by fungal species coexisting within a specific environmental or biological niche is commonly called as the “mycobiome” [13]. However, in CSD, the mycobiome remains poorly understood, particularly in comparison to the microbiome. The interaction between microbial communities in the reproductive tract, particularly bacteria and fungi, is critical in maintaining female reproductive health [14].

Therefore, this study aims to examine the composition and function of fungi in CSD, as well as their interactions with bacteria, to provide a more comprehensive understanding of the role of microbial communities in CSD. This study also aims to identify potential therapeutic approaches to enhance fertility by building on current research.

Results

Recruitment of participants and metagenomic sequencing

Forty-eight participants were included in this study, including 24 in the CSD group and 24 in the CON group. The clinical characteristics of the participants are shown in Table 1. To better characterize the composition of the cervical microbiota, we performed ultra-high-depth metagenomic sequencing. The effective data volume of each sample ranged from 14.25 to 18.63 G, the N50 statistical distribution of contigs ranged from 230 to 245 bp, and the number of ORFs in the gene catalogue (non-redundant gene set) constructed after redundancy removal was 281,107.

Table 1.

Composition and characteristics of bacterial communities

Alpha diversity based on Pielou index (Fig.1 A) and Shannon index (Fig.1 B) calculations showed that species richness and evenness in the CSD group were significantly higher than those in the control group (T-test). Beta diversity based on the Bray-Curtis distance (Fig.1 C) indicated that the distances between samples within the CSD group were greater than those in the CON group (PERMANOVA test).

The composition of CST community types in the CSD group and CON group differed. The proportion of CST III decreased in the CSD group, while the more unstable CST IV-B and CST IV-C increased (Fig.1 D). The CST community types in the CON group were dominated by Lactobacillus spp. (Fig.1 D). CST subtype analysis revealed more detailed community types. In CST I (dominated by L. crispatus) and CST III (dominated by L. iners), the CSD group was mainly distributed in the B subtypes with a lower abundance of L. crispatus and L. iners (Fig.1 E). The stacked bar chart of species composition shows the species composition characteristics of each sample (Fig.1 F). One hundred and sixty-two differential species were identified between the two groups. The differential species analysis indicated that Bifidobacteriaceae bacterium spp. and Gardnerella spp. were significantly higher in the CSD group than in the CON group (Fig.1 G).

Fig. 1.

Composition and characteristics of fungal communities

A total of 431 fungal species were annotated (Supplementary File 1). There was no difference in α diversity based on the Pielou index (Fig.2 A) and Shannon index (Fig.2 B) between the two fungal communities. β diversity suggested similar distances between samples within each group (Fig.2 C). Regarding species composition, the composition of the two groups was similar at the phylum (Fig.2 D) and species (Fig.2 E) levels. The tree map revealed the relationships among fungal phylum (Fig.2 F). Forty-two differential species were identified between the two groups, and the top 10 species are shown in Fig.2 G.

Fig. 2.

The co-occurrence network among microbial community

We constructed co-occurrence networks for fungal species in the CSD and CON groups (Fig.3 A). The network in the CSD group had fewer connections and weaker random robustness than the CON group (Fig.3 B), indicating greater fragility of the CSD fungal co-occurrence network. Subsequently, a cross-domain network between bacteria and fungi was constructed (Fig.3 C). The connection between bacteria and fungi in the CSD group was reduced compared with that in the CON group. At the same time, the CSD group had weaker random robustness than the CON group (Fig.3 D). We selected species with a sample coverage rate greater than or equal to 60% for Spearman correlation analysis and visualized the data using heatmaps (Fig.3 E) and networks (Fig.3 F). The results showed that Lachnellula suecica, Arthrobotrys oligospora, and Piptocephalis cylindrospora have close relationships with bacteria (Fig.3 F). Piptocephalis cylindrospora has a close symbiotic relationship with Lactobacillus jensenii (Fig.3 F).

Fig. 3.

Functional gene analysis

Functional gene analysis indicates significant differences in the functional gene composition of microbial communities between the CSD and CON groups (Figure 2-figure supplement 1 A). The CSD group has two significantly higher KEGG modules than the CON group, including M00142 (NADH: ubiquinone oxidoreductase, mitochondrion) and M00151 (cytochrome bc1 complex respiratory unit) (Figure 2-figure supplement 1 B). The CSD group has significantly increased activity in the inflammatory mediator regulation of TRP channels and platelet activation, while the activity of the galactose metabolic process has decreased (Figure 2-figure supplement 1 C). The CAZy annotation results indicate that the activity of GH73 (glycoside hydrolase family 73) in the CSD group has decreased (Figure 2-figure supplement 1 D).

Figure 2-figure supplement 1

Untargeted metabolomics revealed unique metabolic characteristics

LC/MS untargeted metabolomics was performed on 40 of the subjects. The heatmap shows differences in the distribution of metabolites between the two groups (Fig.4 A). Enrichment analysis of signaling pathways suggests that the CSD group has significantly increased activity in the Thermogenesis and Pentose phosphate signaling pathways, while the activity in Steroid biosynthesis is significantly decreased (Fig.4 B). Spearman correlation analysis suggests a close correlation between changes in the metabolite contents of Lachnellula suecica and Arthrobotrys oligospora (Fig.4 C). At the same time, we also explored the relationship between bacterial and metabolite content changes (Figure 4-figure supplement 1).

Fig. 4.

Figure 4-figure supplement 1

Integrating bacterial, fungal, and metabolite data, we found that the fungus Lachnellula suecica has a mutual regulation relationship with Lactobacillus spp. through GoyaglycosideA and Janthitrem E (Fig. 5). At the same time, Lachnellula suecica occupies a central position in the entire regulatory network and is closely related to changes in metabolites and bacterial composition. Clavispora lusitaniae and Ophiocordyceps australis are also fungal species that are closely related to changes in the abundance of Lactobacillus spp. (Fig. 5).

Fig. 5.

Discussion

Caesarean section scar diverticulum (CSD) significantly affects female fertility, posing a challenge to women who wish to conceive. Our research team conducted a preliminary investigation on bacteria-host interactions [12]. Current findings suggest that fungi play crucial roles in maintaining bacterial community and microbiome stability [15, 16]. This study explores the intricate fungal-bacterial interactions within the microbiome of CSD, building on our team’s previous research findings.

We achieved more precise bacterial community compositions by implementing higher-resolution metagenomics than our previous study [12]. Lactobacillus jensenii of the Lactobacillus genus exhibited the greatest discrepancy between the two groups. This species has been linked to integral aspects of women’s reproductive health. It has been demonstrated that a reduction in L. jensenii abundance has a close correlation with early-stage embryo arrest [17]. Furthermore, research suggests that L. jensenii may hold potential in facilitating vertical transmission from mother to infant [18].

Fungi and bacteria, as integral components of the human microbiome, establish complex interactions with each other and the host. The characteristics of fungal communities in CSD an CON group are similar, but there are differences in species composition. Clabispora lusitaniae is the species with the largest difference between the two groups. In addition, the co-occurrence network of fungal communities and the stability of the bacterial-fungal cross-domain network in the CSD group are poorer than those in the CON group, indicating that the cervical microbiota stability is disrupted in CSD, resulting in local microbial-immune imbalance. A previous study confirmed that persistent exudation in CSD is a crucial factor affecting embryo implantation [7]. Several studies have shown that microbial dysbiosis can affect local immune balance and cause inflammation [19, 20].

Fungi play an important role in maintaining the stability of bacteria and the entire microbiome. Studies have shown that adding specific fungi can promote the growth of certain bacteria [16], while fungi also play an important role in improving the adaptability of bacteria [21]. To further understand how abnormal fungal community status in CSD affects bacteria, leading to the disruption of local stability, we used LC/MS technique for metabolite detection. The metabolite analysis results showed differences in the metabolite spectra between the CSD and CON groups, and the regulatory mode of several different fungal species on metabolites was mainly negative regulation. By integrating fungal, bacterial, and metabolite information, we found that Lachnellula suecica played an important role in regulating the abundance of several Lactobacillus species. Lachnellula suecica reduces the abundance of several Lactobacillus species, including Lactobacillus jensenii, by decreasing the metabolites Goyaglycoside A and Janthitrem E. Meanwhile, Clavispora lusitaniae and Ophiocordyceps australis can also synergistically affect the abundance of Lactobacillus by influencing metabolite abundance.

So far, we have conducted a comprehensive analysis of bacterial and fungal communities and changes in environmental metabolites in Caesarean section scar diverticulum (CSD). While our study provides significant insights, it is important to note that our methodology did not include the use of negative controls. This is a notable limitation, as negative controls are critical in distinguishing true microbial signals from potential contamination, particularly for rare species. The absence of negative controls could potentially affect the discrimination and identification of such rare species within our samples. Consequently, while our results provide valuable information on the prevalent microbial interactions within CSD, they should be interpreted with caution, especially when considering the implications for clinical treatment plans.

Furthermore, our analysis did not extend to examining the antibiotic sensitivity of the key fungi and bacteria identified, which limits the direct application of our results to clinical treatments. Future studies should consider including negative controls and antibiotic sensitivity analysis to validate our findings and extend the implications for clinical practice.

Despite these limitations, the study’s results highlight the significant role of fungal composition and activity in driving bacterial dysbiosis in CSD. This insight is crucial for informing clinical treatment strategies and suggests that addressing fungal populations could be an important aspect of managing CSD. Future research is needed to further elucidate the complex interactions within the CSD microbiome and to develop more comprehensive treatment approaches that address both bacterial and fungal components.

Materials and methods

Participant recruitment and sample collection

The subjects were recruited at the Reproductive Medical Center of the Sixth Affiliated Hospital of Sun Yat-Sen University. The inclusion and exclusion criteria were previously described in a previous study [12]. All participants were fully informed, volunteered to participate in this study, and signed informed consent forms. All study processes were reviewed and approved by the Sixth Affiliated Hospital Ethics Committee of Sun Yat-Sen University (IRB no. 2019ZSLYEC-005S).

The sample collection procedure is as described above [12]. In short, after cleaning the external genitalia and vagina, a sterile disposable swab is inserted into the cervical canal, and the swab is rotated five times to collect the sample fully. The sample is rapidly transferred to liquid nitrogen for quenching and then stored at -80°C.

DNA extraction and metagenomic sequencing

According to the manufacturer’s instructions, the total DNA was isolated from the sample using QIAamp® Fast DNA Mini Kit (Qiagen, Hilden, Germany). The DNA concentration and integrity were evaluated by NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis. The DNA was fragmented with S220 Focused-ultrasonicators (Covaris, USA) and purified with Agencourt AMPure XP magnetic beads (Beckman Coulter Co., USA). Then, the library was constructed using the TruSeq Nano DNA LT Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. The sequencing was performed using Illumina NovaSeq 6000.

Raw Data Preprocessing

Parallel computing in research is implemented with GNU Parallel [22].

The sequencing raw data was stored in FASTQ files. After trimming and filtering with Trimmomatic (v0.36) [23], the clean paired-end reads were aligned to the host genome (hg38) with bowtie2 (Presets: --very-sensitive-local) [24], and the matched reads were discarded. The clean reads were assembled into contigs using MEGAHIT software (--k-list 21,33,55,77) [25] with a filtering threshold of 500 bp. The ORFs were predicted using Prodigal software [26] and were translated into amino acid sequences. CD-HIT software [27] was used to remove redundancy from the ORF predictions of each sample and mixed assembly with default parameyers, and to obtain non-redundant initial Unigenes. The clustering was performed with default parameters of 95% identity and 90% coverage, and the longest sequence was selected as the representative sequence. Bowtie2 software was used to align the clean reads of each sample to the non-redundant gene set (95% identity), and the abundance information of genes in each corresponding sample was calculated from the number of reads and the length of the genes.

Taxa and functional annotation

The representative sequences of the redundant gene set were aligned with NCBI’s NR database using DIAMOND software [28], and annotations with e <1e-5 were selected to obtain proteins with the highest sequence similarity, thereby obtaining functional annotations and species annotations. Species present in only <10% of samples were excluded. The Valencia software [29] was used to perform community state type (CST) assignment for each sample. The abundance of species is calculated using the sum of the gene abundance corresponding to the species. LEfSe (Linear discriminant analysis Effect Size) [30] was used to screen differentially abundant species by using MicrobiotaProcess R package [31]. The first comparison uses the Kluskal-Wallis test, the second comparison uses the Wilcoxon test, and finally LDA is used to determine the difference, and the p value is corrected using fdr. The calculation and visualization of species diversity and abundance are executed by the EasyMicroPlot R package [32] and EasyAmplicon pipline [33]. Annotations were performed using the KO, COG, and CAZy databases to explore changes in enzyme family related to the microbiome. Differential functional analysis was conducted using the STAMP software [34].

Construction of co-occurrence network

In network construction, we used more stringent species filtering criteria and only included species detected in 50% of samples in the analysis. The construction and random robustness comparison of fungal and bacterial-fungal co-occurrence networks were executed by the ggClusterNet R package [35]. Spearman correlation analysis was used to calculate correlations, the cutoff value of the correlation coefficient was set to 0.6, and p values were corrected using FDR. The stability analysis of the network is implemented by the Robustness.Random.removal function of the ggClusterNet R package with the default parameters.

Metabolite detection using liquid chromatography (LC) and mass spectrometry (MS)

A separation buffer (methanol/acetonitrile/water [2:2:1]) was used to extract metabolites from cervical swabs. Ten ul was taken from each sample and mixed to form a quality control (QC) sample [36], which was used to evaluate stability during the experiment. The metabolite signals in each sample were detected using liquid chromatography (LC) and mass spectrometry (MS) (ACQUITY UPLC I-Class plus, Waters). Progenesis QI v2.3 software was used for qualitative analysis of metabolites, with parameters of 5ppm precursor tolerance, 10 ppm product tolerance, and 5% product ion threshold (Nonlinear Dynamics, Newcastle, UK). All output data were normalized using internal standards, and the results were displayed as peak values (test sample peak area/internal standard sample peak area). Compound identification was performed using human metabolome database (HMDB), Lipidmaps (V2.3), Metlin, EMDB, PMDB, and a self-built database based on mass-to-charge ratio (M/z), secondary fragmentation, and isotope distribution.

The screening of differential metabolites was performed using the OPLS-DA (orthogonal partial least squares discriminant analysis) method, with a VIP value>1 for the first principal component and a P value<0.05 for T-test as the threshold. MetPA was used for differential metabolic pathway analysis [37]. We used Spearman correlation analysis to evaluate the correlation between the microbial community and metabolites. The network was visualized using Cytoscape software [38].

Declarations

Ethics approval and consent to participate

All study procedures were reviewed and approved by the ethics review board of the Sixth Affiliated Hospital of Sun Yat-Sen University (IRB no. 2019ZSLYEC-005S).

Authors’ contributions

Xing Yang and Peigen Chen coordinated the study, participated in the design. Xing Yang carried out the study. Peigen Chen and Haicheng Chen analyzed and interpreted the data and drafted the manuscript. Xinyi Pan, Qianru Liu and Ziyu Liu recruit subjects. Xinyi Pan and Ziyu Liu participated in the follow-up work and data collection. Qianru Liu were also responsible for preparing ethical review materials. Xing Yang did the review of the final paper. All authors read and approved the final manuscript.

Competing interests

The authors declare that there are no conflicts of interest to disclose.

Data Availability

The metagenome sequencing has been deposited in China National Center for Bioinformation (https://ngdc.cncb.ac.cn/) under reference number PRJCA016850.

Funding

This work was supported by the [National Key R&D Program of China] under Grant [number 2022YFC2703000]; [National Natural Science Foundation of China] under Grant [number 82271651]; [Nationally Funded Postdoctoral Researcher Program] under Grant [number GZC20233216].

Acknowledgements

The authors thank to OE Biotech Co., Ltd. (Shanghai, China) for assistance with metagenomic sequencing.