Sample Collection
The study enrolled a total of 194 women diagnosed with infertility at the Reproductive Medicine Center of Shanghai Changzheng Hospital between November 2018 and November 2021, comprising the infertility group. Additionally, 102 healthy women who underwent routine physical examinations at the Physical Examination Center of the same hospital were recruited as the healthy control group.
Inclusion criteria for the infertility group were: regular unprotected sexual activity for at least one year without achieving pregnancy; Inclusion criteria for the healthy group were: women with a history of childbirth and no history of infertility and other gynecological diseases.
Exclusion criteria for both groups included the presence of an IUD, vaginal inflammation, any acute inflammatory condition, suspicion of cervical or endometrial neoplasia, and endocrine or autoimmune disorders. Subjects were also excluded if they had recently used hormones, antibiotics, or vaginal medications; undergone cervical treatment, endometrial biopsy, IUD removal, or hysteroscopy within the past week; performed douching within five days; or had sexual activity within 48 hours prior to sampling. Additionally, none of the participants were pregnant, lactating, or menstruating at the time of sampling.
Age, current residence, menstrual history, and fertility history were collected from all participants. For nulliparous subjects, they were asked about their pregnancy plans. All participants provided written informed consent, and the study received approval from the Ethics Committee of Shanghai Changzheng Hospital(2023SL051).Vaginal secretions were obtained from the posterior fornix region of both groups. The samples were then rapidly frozen using liquid nitrogen, stored at − 80°C, and transported in dry ice to BGI-Shenzhen for further analysis.
DNA Extraction and 16S rRNA Amplicon Sequencing
DNA was extracted from vaginal secretions and amplified by PCR using high quality DNA. More accurate distance-based clustering of reads was achieved by using V4-V6 primers and PCR premixes with targeted primers for these regions. Subsequently, these clusters were categorized as amplicon sequence variants (ASVs) at the species level[18]. Finally, PCR products were purified using Agencourt AMPure XP magnetic beads and dissolved in elution buffer. The fragment ranges in the libraries were assessed by an Agilent 2100 Bioanalyzer and quality control libraries were sequenced on the HiSeq 2000 platform).The primers for the V4–V6 regions were listed as follows: 8F-‘AGAGTTTGAT[YM]TGGCTCAG’, 518R-‘ATTACCGCGGCTGCTGG’. Y and M represent bases C/T and C/G, respectively.
Quantitative Real-Time Polymerase Chain Reaction
Real-time PCR detection was carried out by using primers for amplifying 16s rRNA gene and beta-actin. The primers were listed as follow[19] : Lactobacillus crispatus: forward primer 5′-AGCGAGCGGAACTAACAGATTTAC-3′, reverse primer 5′-AGCTGATCATGCGATCTGCTT-3′; Lactobacillus iners: forward primer 5′-AGTCTGCCTTGAAGATCGG-3′, reverse primer 5′-CTTTTAAACAGTTGATAGGCATCATC-3′; beta-actin: forward primer 5′-AAAAGCCACCCCACTTCTCT-3′, reverse primer 5′-CTCAAGTTGGGGGACAAAAA-3′. The 20-µl PCR mixture contained 1 µl of DNA sample, 1 µl of each primer, 6 µl of ultra-pure water and 12 µl of 2*SYBR Green Mix. The Eppendorf realplex system (Eppendorf, USA) was used with the thermal cycling profile of 95°C for 5 min, and 40 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s. Each sample has three technical duplicates. The bacterial abundance was calculated by dividing the average CT value of 16s rRNA gene by the average CT value of beta- actin.
Cleaning the Raw Sequencing Data
The steps to process the raw data are as follows: 1) discard readings with low base quality: set 30 bp as the window length, truncate the end sequence of readings from the window if the average quality of the window is less than 17, and delete readings whose final reading length is less than 75% of the original reading length; 2) discard readings contaminated by adapters: the default adapter sequence overlaps with the read sequence at 15 bp, sets it to 15 bp, and allows mismatch of 3. 3) eliminate readings with Ns; 4) eliminate readings with low complexity: the length of continuous occurrence of bases in the readings is ≥ 10. The resulting data in fastq format were termed as clean data.
Amplicon Sequence Variants and Taxonomy Analysis
The 16S rRNA clean data were processed by DADA2 package in R[20]. DADA2 provides a sensitive and specific workflow in amplicon sequencing. The DADA2 pipeline proceeds as follows: 1) filter and trim clean data: discarding reads at the first instance of a quality score less than or equal to 2; 2) remove duplicated sequence entries in fastq files; 3) merge paired reads; 4) learn the error rates and infer the sample composition using the error rates; 5) construct a sequence table; 6) remove chimeras; and 7) assign taxonomy using naive Bayesian classifier method and the RDP_16sRNA reference databases.
Normalizing the Relative Abundance of Bacteria
For each sample, we divide the original count of bacteria by the total reading count, thus normalizing the abundance of each bacteria at the genus or species level. The proportion of each type of bacteria is used as the normalized abundance. In addition, bacteria that account for less than 0.5% are merged into "other". Bacteria that account for more than 0.5% of at least two samples will be retained for further data analysis.
Diversity Analysis
The α-diversity analysis was performed with the R software package. Based on the OUT table, we calculated the observed species, Shannon index, chao 1 index and ACE value to estimate the α-diversity of the colony and performed principal component analysis (PCA) based on the sequencing data. The Chao index and observed species reflect the abundance of OTUs in the samples, while the Shannon index and Simpson index reflect the diversity of OTUs in the samples. PCOA belongs to a type of beta diversity, which is usually used to represent the material differences between different environmental communities.
The Functional Prediction of Microbiota
Tax4fun[21] package in R was used to estimate Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog (KO) scores for each sample, and further used the scores to predict the relative activity of the KEGG pathway.
The Differential Abundance Analysis
Wilcoxon rank-sum test was used to compare the abundance of microflora between the two groups. The p value was adjusted by Benjamini and Hochberg methods to avoid multiple tests, and Kruskal-Wallis test was used to make multiple comparisons between groups.