This study was conducted at Shanghai First Maternity and Infant Hospital affiliated with Tongji University from March 2018 to July 2020. The inclusion criteria of the participants were as follows: (a) women aged between 40 and 69 years; and (b) subjects undergoing hysterectomy by any standard surgical approach for any possible benign disease and at stage I endometrial cancer (EC). The exclusion criteria were as follows: (a) pregnant women and nursing mothers; (b) women who used antibiotics within 3 months; (c) history of genital tract infection or medication within 3 months; and (d) patients receiving preoperative chemotherapy or radiotherapy. A total of approximately 120 subjects were recruited, and after filtering with the inclusion and exclusion criteria, 40 women, including 30 endometrial cancer (EC) patients and 10 healthy controls (HCs), were finally enrolled ( Table S1 ).

Endometrial tissue samples were acquired from these 40 women, who had undergone a hysterectomy. After removing the corrupted and erosive tissues of the surface, the residual tissues were macrodissected. The whole process was performed by an experienced gynecologist under strict aseptic procedures. The samples were divided into three equal parts (one for 16S rRNA sequencing, one for RNA-seq, and one for subsequent experimental verification), snap frozen and stored at −80°C.

Bacterial DNA was isolated from the endometrial samples using a DNeasy PowerSoil kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. DNA concentration and integrity were measured by a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively. PCR amplification of the V3-V4 hypervariable regions of the bacterial 16S rRNA gene was carried out in a 25 μL reaction using universal primer pairs (343F: 5′-TACGGRAGGCAGCAG-3′; 798R: 5′-AGGGTATCTAATCCT-3′) (Chen et al., 2018). The reverse primer contained a sample barcode, and both primers were connected with an Illumina sequencing adapter.

The amplicon quality was visualized using gel electrophoresis. The PCR products were purified with Agencourt AMPure XP beads (Beckman Coulter Co., USA) and quantified using a Qubit dsDNA assay kit. The concentrations were then adjusted for sequencing. Sequencing was performed on an Illumina MiSeq with two paired-end read cycles of 300 bases each (Illumina Inc., San Diego, CA; OE Biotech Company; Shanghai, China). The raw data for the Illumina reads from 40 endometrial samples were uploaded to the Sequence Read Archive (SRA) database.

Paired-end raw sequences were preprocessed using Trimmomatic software (version 0.36) (Bolger et al., 2014) to detect and cut off ambiguous bases (N), and this software was also used to cut off low-quality sequences, with an average quality score below 20 using a sliding window trimming approach. After trimming, paired-end reads were assembled using FLASH software (version 1.2.11) (Reyon et al., 2012). Sequences were further denoised using QIIME software (version 1.8.0) (Caporaso et al., 2010). Then, reads with chimeras were detected and removed using VSEARCH (version 2.4.2) (Rognes et al., 2016). Clean reads were subjected to primer sequence removal and clustering to generate operational taxonomic units (OTUs) using VSEARCH software, with a 97% similarity cutoff (Rognes et al., 2016). The representative read of each OTU was selected using the QIIME package. All representative reads were annotated and blasted against the Silva database (Version 123) using the RDP classifier (confidence threshold was 70%) (Wang et al., 2007). The microbial diversity in the endometrial samples was estimated using the alpha diversity, which includes the Chao1 index, observed species index, Simpson index, and Shannon index (Yang et al., 2019). The unweighted UniFrac distance matrix performed by QIIME software (version 1.8.0) was used for nonmetric multidimensional scaling (NMDS) plots. A Venn diagram was generated to visualize the shared and unique genera among groups regardless of their relative abundance using the R package “VennDiagram” (Zaura et al., 2009). The LEfSe (linear discriminant analysis effect size) method was used to detect differentially abundant taxa across groups using the default parameters (Segata et al., 2011). 16S rRNA gene amplicon sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China).

Total RNA was extracted from frozen endometrial tissues by using mirVana™ miRNA Isolation Kit (catalog no., AM1561, Ambion^®), and cDNA libraries were sequenced on the MGISeq2000 Platform (MGI, Shenzhen, China). Then, 150 bp paired-end reads were generated. Approximately 49.35 million raw reads for each sample were generated. Raw data (raw reads) in fastq format were first processed using Trimmomatic software (Bolger et al., 2014), and the low-quality reads were removed to obtain clean reads. Then, approximately 48.57 million clean reads for each sample were retained for subsequent analyses.

The clean reads were mapped to the human genome (GRCh38) using HISAT2 (version 2.2.1.0) (Kim et al., 2015). The FPKM (Roberts et al., 2011) value for each gene was calculated using Cufflinks (version 2.2.1) (Trapnell et al., 2010), and the read counts of each gene were obtained by HTSeq-count (version 0.9.1) (Anders et al., 2015). A differential expression analysis was performed using the DESeq (2012) R package (version 3.2.0, R Foundation for Statistical Computing, Vienna, Austria). A p value < 0.05 and fold change > 2 or fold change < 0.5 were set as the thresholds for significantly differential expression. A hierarchical cluster analysis of differentially expressed genes (DEGs) was performed to demonstrate the expression pattern of genes in different groups and samples. GO enrichment and KEGG pathway enrichment analyses of DEGs were performed using R based on the hypergeometric distribution.

The Pearson correlation coefficient was calculated to evaluate the clinical characteristics that underwent significant changes (values of serum D-dimer (DD), fibrin degradation products (FDPs), cancer antigen (CA) 125, CA199, and total protein) and endometrial microbiota at the genus level between HC and EC women. To reduce the computational load and avoid contingency, the genera that presented an abundance value of ≤ 0.1% were excluded. The significance of each genus-hematological indicator pair was determined based on a p value < 0.05 and a false discovery rate (FDR) < 0.1. Detailed scripts of the correlation calculations are provided in Table S2 .

The Pearson correlation coefficient was calculated to measure the connections between genus abundance and differential gene expression level for each genus-gene pair across all 30 patients. To reduce the computational load and avoid contingency, the genera that presented an abundance value of ≤ 0.1% were excluded. The significance of each genus-gene pair was determined based on a p value < 0.05 and FDR < 0.1. Independent Student’s t test was applied to evaluate the difference in log₂FC values calculated by DESeq between the HC and EC groups. Detailed scripts of the correlation calculations are provided in Table S3 .

Paraffin-embedded endometrial tissue sections were treated as previously described (Li et al., 2019). After retrieval of antigenic binding sites, the sections were permeabilized with 0.2% Triton X-100 for 10 min, and then blocked with 5% bovine serum albumin for up to 1.5 h at room temperature, and incubated with appropriate concentrations of primary antibodies against PRSS33 (Abmart, Shanghai, China, TD4461M, 1:200 dilution), CPB2 (Bioss Antibodies, Beijing, China, BS-7547R, 1:200 dilution), and XBP1 (Abcam, London, UK, ab37152, 20 µg/mL) overnight at 4°C. Then, the sections were incubated with Alexa Fluor^® 488/594/647 conjugated goat anti-rabbit IgG (Abcam, London, UK, 1:200 dilution) secondary antibodies for 1 h. Cell nuclei was stained with DAPI. The stained sections were observed under a fluorescence confocal microscope.

Fundamental statistical analyses were performed using SPSS software (version 19.0, IBM SPSS Inc., Chicago, IL, USA), GraphPad Prism 7 software (GraphPad software, Inc., San Diego, California, USA) and R version 3.2.0. Categorical variables were expressed as numbers (%). For normally distributed variables, the Student’s t test was used; otherwise, the Wilcoxon rank sum test was used. Pearson’s rank correlation analysis was conducted to calculate the correlation between genera, between genera and clinical characteristics, or between genera and DEGs. All statistical tests were two-tailed, and the p value was adjusted by the Benjamini-Hochberg correction. A p value < 0.05 was considered statistically significant.

All 30 endometrial cancer (EC) patients and 10 healthy control (HC) women were Han Chinese individuals from Shanghai. All diagnoses were performed based on the final surgical pathology. The demographics and other clinicopathological characteristics are shown in Table 1 . The control participants were generally matched with the case participants by age and BMI, suggesting that there were no established confounding factors that might influence group discrimination prior to sample collection. There was a significant increase in serum DD (p = 0.0122), FDP (p = 0.0003), CA125 (p = 0.0094) and CA199 (p = 0.0261) levels and a significant decrease in serum total protein in EC patients (p = 0.0441) ( Table S1 ).

We analyzed a total of 40 endometrial samples based on the pyrosequencing of barcoded 16S rRNA genes (V3–V4) and finally obtained 2,728,218 qualified sequences (median =  68,151) and 9,278 OTUs ( Tables S4 , S5 ). The rarefaction curves indicated that the estimated richness of OTUs almost reached saturation in the HC (n = 10) and EC groups (n = 30) ( Figure 1A ). Moreover, the rank abundance curves indicated decreased richness in the EC patients compared with the HCs ( Figure S1 ).

To estimate the differences in bacterial diversity between the two groups, the sequences were aligned to estimate alpha diversity and beta diversity ( Table S6 ). Significant differences between the EC and HC groups were observed in the Chao1 (2247.72 ± 151.01 versus 2476.16 ± 142.47, p < 0.001), observed species (1321.61 ± 114.95 versus 1420.98 ± 117.80, p = 0.0212) and Simpson (0.989 ± 0.007 versus 0.993 ± 0.001, p = 0.0053) indexes ( Figures 1B–D ), whereas significant differences between these groups were not observed for the Shannon index (7.99 ± 0.41 versus 8.18 ± 0.26, p = 0.0606) ( Figure S2 ). A Venn diagram showed that 100 out of 147 genera were shared between the EC and HC groups ( Figure 1E ). Notably, 20 out of 147 genera were unique to the EC group. For the beta diversity, the NMDS plots evaluated by the unweighted UniFrac distance matrix exhibited a separation of the distribution of the two groups ( Figure 1F ). The results suggest that the diversity of endometrial microbiota could be strongly influenced by the tumor burden.

The bacterial distribution was assessed according to the relative abundance of different taxa. The five dominant phyla in the EC and HC groups were Bacteroidetes, Proteobacteria, Firmicutes, Actinobacteria, and Cyanobacteria, which accounted for 92.17% and 92.67% of the total OTUs, respectively ( Figure 2A ). The compositions of the bacterial community (top 15) at the genus level are shown in Figure 2B . To provide insight into the differences in endometrial microflora between the two groups, the Wilcoxon rank sum test was performed at both the phylum and genus levels ( Tables S7 , S8 ). Compared with the HC group, Pelomonas, Prevotella, Nocardioides and Muribaculum were enriched in the EC group, whereas Oscillibacter exhibited the opposite trend (all p < 0.05, Figure 2C ).

Next, the LEfSe method was used to discover high-dimensional biomarkers. Of the above five genera, Pelomonas and Prevotella were significantly overrepresented [LDA scores (log₁₀) > 3] in the endometrium of EC patients ( Figure 2D ). The relative abundances of the two genera were further subjected to a cluster analysis, and the results suggested that Pelomonas and Prevotella were abundant in the EC group ( Figure 2E ). Sankey diagrams were also generated to visualize the major proportion of taxa (phylum and genus) between the HC and EC groups. In patients with a tumor burden, the proportion of Pelomonas and Prevotella gradually increased and became predominant, accompanied by changes in other bacteria ( Figure 2F ). Overall, these data suggest that the differentially abundant microbiota, including Pelomonas and Prevotella, were able to distinguish the microbiota between the HC volunteers and EC patients.

For the intraindividual correlation analyses, we tested the associations between clinical characteristics (abnormal values of serum DD, FDPs, CA125, CA199, and total protein) and endometrial microbiota. Pearson’s correlation-based analysis indicated that the three genera Prevotella, Acinetobacter and Brevundimonas were remarkably associated with fibrin degradation (DD and FDPs), while meaningful results based on other clinical characteristics were not observed ( Figure 3A , Table S2 ). Among the three genera, Prevotella microbes were correlated with the aforementioned different endometrial microbiota composition alterations in the HC and EC groups ( Figures 3B, C ). Therefore, Prevotella were selected for further study.

Since the EC patients demonstrated signature microbiota associated with clinical characteristics (tumor burden), we hypothesized that changes in the transcriptome of endometrial tumorigenesis may be correlated with endometrial microbiota. Thus, we performed a transcriptome analysis of the paired tumor and adjacent nontumor endometrial tissues from 30 EC patients. Based on our definition, we identified a total of 3,911 differentially expressed genes (with 1,852 upregulated and 2,059 downregulated) among the 30 paired endometrial tissues by DESeq ( Figure 4A ). The genes whose log₂FC values calculated by GFOLD were zero in > 90% of patients were excluded, thus leaving 3,410 genes for further investigation of the correlations with Prevotella.

Furthermore, Pearson’s correlation-based analysis was performed to discover Prevotella-associated genes and to examine whether the levels of DD and FDPs could be partially influenced by these genes. Among the 3410 genes, 74 genus-gene pairs were formed with FDR < 0.1 ( Table S3 ). A total of 8 genes were identified to potentially affect the degradation of fibrin (peptidase activity, proteolysis, protein digestion and absorption, and cellular secretory pathway transcription factor) ( Figure 4B ). The pathway analysis confirmed that these 8 genes converged mainly on proteolysis-related pathways ( Figure 4C ). A subsequent correlation analysis showed three examples of Prevotella-gene pairs (Prevotella-serine protease 33 (PRSS33), Prevotella-carboxypeptidase B2 (CPB2), Prevotella-x-box binding protein 1 (XBP1)) ( Figures 4D–F ). The mammalian PRSS33 is primarily responsible for degrading fibrin into DD and FDPs within blood (Izem et al., 2021). CPB2 is involved in the clotting-fibrinolytic pathway, and abnormal expression of this protein has been linked to thromboembolism and cancer (Lip et al., 2002). XBP1 is a key transcription factor of the cellular secretory system, and it is often overexpressed in cancers, such as oral squamous cell carcinoma and hepatocellular carcinoma, and correlates with the clinical outcome (Stelloo et al., 2020). Multiplex immunofluorescence confirmed that PRSS33, CPB2 and XBP1 proteins were significantly highly expressed in the endometrial tumor tissues of patients with elevated DD and FDPs ( Figure 4G ), suggesting that Prevotella could partially promote DD and FDPs by influencing host gene expression. To further validate the other 2 genes (KCNN4 and KLK11) under tumor burdens, DESeq was used to calculate the correlation between Prevotella and the KCNN4 and KLK11 genes ( Figure S3 ). The results showed that these two genes satisfied the definition of differential genes between EC and adjacent nontumor endometrial tissues.

According to the abovementioned analysis, we focused on Prevotella, DD and FDPs as potential biomarkers because they were abundant in the EC group. A receiver operating characteristic (ROC) curve analysis indicated that Prevotella (area under the curve (AUC) = 0.64; p = 0.047), DD (area under the curve (AUC) = 0.73; p = 0.031), and FDPs (AUC = 0.79; p = 0.005) were significantly associated with EC samples ( Figures 5A–C ). In addition, a ROC model in combination with the three could increase the AUC value to 0.86 (p = 0.001) ( Figure 5D ). The data indicated that these three biomarkers related to fibrin degradation have the potential to predict EC occurrence.