The participants with UFs (n = 42) and the control participants (n = 43) were recruited at The Third Xiangya Hospital of Central South University from December 2020 to May 2021. The UF patients were diagnosed by the Gynecology Department of The Third Xiangya Hospital according to the clinical practices (Stewart, 2015). The exclusion criteria included severe chronic diseases (e.g., metabolic disorders, heart failure, cirrhosis, and gastrointestinal, neurological, and/or autoimmune diseases). Moreover, individuals with a history of probiotic interventions, diarrhea, and taking antibiotics or NSAIDs within 3 months before prior to collection were also excluded. None of the female subjects was premenarchal or postmenopausal.

This study was approved by the Ethics Committee of The Third Xiangya Hospital of Central South University and was conducted under the relevant guidelines and regulations (IRB number 22003). Written informed consent had been obtained from the participants before the research, and all samples and questionnaires were voluntary. All UF subjects retained stools at the time of the diagnosis of the disease without initiating any treatment. Moreover, all fecal samples were collected after menstruation. Fresh stool samples were collected in sampling tubes with the preservative solution and stored at -80°C until further processing.

Bacterial DNA was extracted from stool samples using the Stool DNA Extraction Kit (Halgen, Ltd.) according to the procedures described in the manufacturer’s instruction. Primers F44 (RGTTYGATYMTGGCTCAG) and R1543 (GGNTACCTTKTTACGACTT) were used to amplify the DNA of the V1–V9 regions of the 16S rRNA gene. Approximately 20–30 ng of the extracted DNA was used in a 50-µl PCR reaction under the following cycling conditions: 94°C for 3 min for an initial denaturing step; then, followed by 94°C for 30 s, 55°C for 30 s, 72°C for 60s, altogether for a total of 30 cycles; and followed by a final extension step of 72°C for 3 min. Agarose gel electrophoresis was then applied to check if the PCR amplification was successful. Finally, the PCR products were directly labeled using a DNA labeling kit (Halgen Ltd., Zhong Shan, China) and further processed for microarray hybridization.

The human gut bacterial microarrays used were designed and manufactured by Halgen Ltd. The arrays use its proprietary oligo-array technology and cover more than 95% of the culturable gut microbial species found in different populations. Probes were selected from all the variable regions of bacterial 16S rRNA. The length of each probe was designed to be approximately 40 bp. The hybridization mixture consisted of 500 ng of Cy5-labeled test sample DNA and 50 ng of Cy3-labeled reference pool. Then, the hybridization buffer and the Cy3- and Cy5-labeled samples were added to a final volume of 150 µl, heated to 100°C for 5 min, and cooled on ice for 5 min. All hybridization mixtures were placed in a hybridization box and then hybridized at 37°C for 3.5 h in a hybridization oven. Finally, the slides were washed in 2× saline sodium citrate, 0.25% Triton X-100, 0.25% sodium dodecyl sulfate, and 1X Dye Protector for 15 min at 63°C. Then, the slides were rinsed in 1X Dye Protector until they were clear of water droplets after immediate withdrawal from the solution. The slides were immediately scanned using a dual-channel scanner.

The Cy5/Cy3 ratio measured by the respective channels was used to calculate the percentage of each microbial species, which is presented as the relative abundance value. R (v.4.1.2) was used in this study. Alpha diversity (α-diversity) and beta diversity (β-diversity) indices were analyzed to characterize species diversity within and among habitats respectively, to evaluate their overall diversity in an integrated manner. α-Diversity includes richness, Shannon–Wiener diversity, Gini–Simpson diversity (obtained by subtracting the value of the classical Simpson index from 1), and Pielou’s evenness. They were measured using the function diversity in the package “Vegan” based on unpumped flat OTU table; the beta diversity indices of the microeukaryotic communities were calculated using the function vegdist after data were pumped. ANOSIM was chosen to test for significance between groups (Wang et al., 2008). For non-metric multidimensional scaling (NMDS), stress less than 0.2 indicates that the results of the NMDS analysis were scientifically credible (Legendre and Legendre, 1998). The species composition at the phylum level between the two groups was visualized using the ggplot2 package. The DEseq2 package was used to analyze species difference and marker species (Topper et al., 2017). The differential expression matrix and the P-value matrix of species composition were obtained through the function DESeqDataSetFromMatrix. The significance level was P <0.05, and the absolute FoldChange value was greater than 2; the volcano map was drawn using the ggplot2 package. The coexistence network of the two groups was established based on Spearman correlation matrix and corrected by P-value matrix using the igraph package; the Benjamini and Hochberg false discovery rate was used to correct the P-value. Modules are divided according to the high intra-module connectivity and the low inter-module connectivity; the Spearman correlation coefficient and corrected P-values were 0.4 and 0.05, respectively (Yuan et al., 2021). Then, the bacteria coexistence network was constructed in Gephi software (https://gephi.org/). A random network with the same number of nodes and edges as the real network was constructed by using the erdos.renyi.game function (Erdos and Renyi, 1960). The classified information of species in network modules was presented using ggplot2. The redundancy analysis (RDA) of the effect of age, body mass index (BMI), and other body indices on the distribution of samples and the distribution of species was performed by ggplot2’s built-in vegan package.

SPSS (version 26.0) was also used in this study. Continuous data were reported by median with range (minimum–maximum) or mean ± standard deviation (SD) and were appropriately analyzed with Wilcoxon test or t-test. The categorical data were described with the number and percentages and were analyzed with χ² test or Fisher’s exact test as appropriate. P < 0.05 was considered to be statistically significant.

The demographic characteristics of the UF group (case group) and non-fibroids (control group) are shown in Table 1 and Supplementary Figure S1. There is no difference between the case and control groups in terms of BMI, menstrual history, previous history, and childbearing history, and a difference in average age was detected (P = 0.02) between the two groups.

All samples were sequenced to sufficient depth and dilution curves, which were calculated and recorded after 5 replicate random samples (Figure 1A). The control group had higher indices than the case group in terms of all α-diversity indices: richness (Figure 1B), Shannon–Wiener (Figure 1C), Gini–Simpson (Figure 1D), and Pielou (Figure 1E). β-Diversity was compared using both the algorithm of NMDS and principal coordinate analysis (PCoA), which demonstrate significant differences between the two groups. The PCoA results showed that the distribution of cases and controls was scattered between groups and clustered within groups with a smaller area in the case group and the ANOSIM test (P = 0.001). The difference between groups was greater than that within groups, implying a significant difference in diversity between the case and control groups (Figure 2A). The NMDS analysis results were similar to those of PCoA and stress = 0.159 (<0.2) (Figure 2B). The α-diversity index of gut microbes was further investigated in patients with different number and locations of tumorigenesis. Interestingly, our results revealed that the gut microbial diversity of patients decreased with the increasing number of tumors (P < 0.01). Some differences were also observed in gut microbial α-diversity indices depending on the location of tumorigenesis (Figures 1F, G and Supplementary Figure S2).

There were commonalities and differences in species composition between the case and the control groups. At the phylum level, 16 phyla were in the case group, while 17 phyla were in the control group (Figure 2C). At the genus level, 337 genera were in the case group, 370 were in the control group, and 321 were common to both groups (Figure 2D). In total, 866 species in the case group and 959 in the control group were at the species level. Among them, 764 were common, while 195 were unique to the control group and 102 to the case group (Figure 2E). In order to demonstrate the differences in taxonomy composition between the case and the control groups, we compared the differentially expressed species between the two groups at the phylum and the species level, respectively. The results suggest that the composition of the fractions was similar, but the abundance percentage of the components varied at the phylum level (Figures 3A, B). Specifically, the relative abundance of Firmicutes, Proteobacteria, Actinobacteria, Cyanobacteria, Dictyoglomi, and Spirochaetes, respectively, was significantly lower in the case group than in the control group (P < 0.05) (Figure 3C). Besides this, among all components, only Verrucomicrobia showed the opposite trend (Figure 3C). Furthermore, we analyzed the differential abundance of species based on statistical differences (Metastats), and the multiplicity of differences characterized the biomarkers (DESeq2). For Metastats analysis, we presented the top 20 differentially expressed species in relative abundance (P < 0.05) (Figure 4A). In total, 17 biomarkers were found, 3 of which (marked in red) were upregulated, while 14 (marked in blue) were downregulated in the case group relative to the controls (Figure 4B).

We performed a network co-occurrence analysis to unravel the relationships among microorganisms. The resulting case network consisted of 863 nodes linked by 17,786 edges, with a much higher number of strong positive correlations (17,311, 97.33%) than negative ones (475, 2.67%), and the average number of edges per node was 2.978. The control network consisted of 958 nodes linked by 23,105 edges, also with a much higher number of strong positive correlations (21,665, 93.77%) than negative ones (1,440, 6.23%), and the average number of edges per node was 2.798 (Supplementary Table S1). The results suggest that microbial networks were made up of closely connected nodes and formed a kind of “small world” topology (Supplementary Table S1). Compared with the topological properties of the random network with the same number of nodes and edges (Supplementary Figure S3 and Supplementary Table S2), the network of the case group exhibited a scale-free characteristic (P < 0.001, Supplementary Figure S4), and the control group also exhibited a scale-free characteristic (P < 0.001, Supplementary Figure S5), indicating that the network structure was non-random. Both the gut microbial interaction networks of the case group and the control group were divided into seven modules (Figures 5A, B). The average degree of the case group is 41.219, which is lower than that of the control group (48.236, P < 0.001, Figure 5C), and the number of sides forming triangles is also lower (P < 0.05, Figure 5D). This suggests that the total connectivity and complexity between gut microbes was higher in the control group than in the case group. The average path length is 2.978 in the case group and 2.798 in the control group (Supplementary Table S1); the average clustering coefficient of the case group is 0.593, which is significantly higher than that of the control group (0.525, P < 0.001, Figure 5E). These results manifest that the average “clustering property” of the whole network between gut microorganisms in the case group was higher than that in the control group. We drew a doughnut to show the taxonomy composition of each module of the case and the control group networks at the phylum level (Figure 5F). The results show that a difference in the components and the abundance percentage exists in each module. Given the above-mentioned findings, it can be concluded that there are differences in the gut microorganism interaction network between the case group and the control group. Compared with the control group, the case group network has lower connectivity and complexity and higher clustering property.