This study was approved by the Ethics Committee of Shanghai Ninth People's Hospital affiliated with Shanghai Jiao Tong University, School of Medicine, China (Document No. 201262), and conducted in accordance with the Declaration of Helsinki. All participants provided written informed consent before enrolment.

The participants constituted two groups: 11 pregnant women who were recruited during their regular visits to the Obstetrics & Gynecology Hospital of Fudan University and seven healthy non-pregnant volunteers from Shanghai Jiao Tong University School of Medicine who participated as the control subjects. The inclusion criteria were as follows: (i) presence of 28 natural teeth excluding third molars; (ii) maintenance of good oral hygiene; (iii) gestation <42 weeks and term delivery without complications for the pregnant group or a regular menstrual cycle for the non-pregnant group; and (iv) good systemic health (i.e., no diabetes, endocrine disorders, or hypertension). The exclusion criteria included the following: (i) antibiotic or steroid hormone use within the preceding 3 months; (ii) untreated oral lesions that could be clinically diagnosed (i.e., dental caries, chronic periodontitis, periapical abscess, or severe halitosis); (iii) previous diagnosis of Sjögren's syndrome or any disease characterized by xerostomia; and (iv) smoking or drinking habits. Demographics were acquired using self-reported questionnaires. Clinical examinations and evaluations were performed at the sites of the following teeth: #3, #8, #14, #19, #24, and #30. The clinical indices included plaque index (PLI), gingival index (GI), and sulcus bleeding index (SBI).

For each of the pregnant women, supragingival plaque samples, and salivary samples were collected during four visits: the first trimester (P1: 11–14 weeks gestation), second trimester (P2: 20–25 weeks), third trimester (P3: 33–37 weeks), and postpartum period (P4: 6 weeks after childbirth). Non-pregnant volunteers also complied with the four visitation times at the same time intervals (NP: NP1–NP4) and were sampled during the luteal phase of the menstrual cycle. One pregnant woman was lost to follow-up after delivery. Thus, 71 samples were gathered at the end of the follow-up period.

Sampling was performed 2 h after eating in the morning, according to the methods described in the Manual of Procedures for the Human Microbiome Project (http://hmpdacc.org/tools_protocols/tools_protocols.php) with minor modifications. After rinsing, 1.5 mL of unstimulated saliva was aspirated with a disposable sterile syringe (with the tip removed) into a sterile Eppendorf tube. Sex hormones including salivary oestradiol and progesterone were measured using a kit from Demeditec (Kiel Wellsee, Schleswig-Holstein, Germany). The selected teeth (#3, #8, #14, #19, #24, and #30) were isolated with cotton rolls and gently air-dried. A sterile Gracey curette was used to collect a pooled supragingival plaque sample from the buccal surfaces of the selected teeth. The collected plaque samples from each participant were separately stored in a sterile Eppendorf tube containing 300 μL of thioglycollate medium, transported on ice to the microbiology laboratory, and kept frozen at −20°C until use.

Total bacterial genomic DNA was extracted using the QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA). The DNA quality was evaluated by measuring the absorbances at 260/280 and 260/230 nm by spectrophotometry (NanoDrop® ND-1000; Thermo Scientific, Wilmington, DE, USA). Only DNA samples with A260/A280>1.7 and A260/A230>1.8 were used for further analysis (Wang et al., 2016).

PCR amplification of the 16S rDNA hypervariable V4–V5 region was carried out using the forward primer 515F (5′-GTG CCA GCM GCC GCG GTA A-3′) and reverse primer 926R (5′-CCG TCA ATT YYT TTR AGT TT-3′). A sample-unique 8-base barcode was incorporated into each forward primer to identify the sequences from the different samples. The PCR programme began with an initial denaturation for 3 min at 98°C, followed by 25 cycles of denaturation (30 s at 98°C), annealing (30 s at 50°C), and extension (30 s at 72°C) with a final extension for 5 min at 75°C. The PCR products were separated by 2% agarose gel electrophoresis, and after purifying with Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN), the final DNA concentrations were determined with the Pico-Green kit (Invitrogen, Carlsbad, CA, USA). Amplicons from different samples were mixed at equimolar concentrations for sequencing, which was performed on an Illumina MiSeq PE300 sequencing platform (San Diego, CA, USA).

Raw sequencing data were processed with the pipeline tools QIIME (version 1.7.0, http://qiime.org/) and MOTHUR (version 1.31.2, http://www.mothur.org/). To retain high-quality sequences, some sequences were eliminated; the eliminated sequences were <150 bp, contained ambiguous bases or homopolymeric stretches, had average Phred scores of <25, or were confirmed as chimeric artifacts (Bokulich et al., 2013). After trimming, the qualified sequences were clustered into operational taxonomical units (OTUs) using the Uclust programme of QIIME at 97% sequence similarity (Edgar et al., 2011) and were aligned against the Greengenes database (Release 13.8, http://greengenes.lbl.gov). After removing the singleton OTUs, a modified OTU table was generated for downstream analysis.

The rarefaction curves were constructed at 97% inter-sequence similarity using QIIME. Venn diagram and alpha diversity indexes were analyzed using MOTHUR, and the significant difference was assessed with the Dunn's multiple comparisons test using the GraphPad Prism software (La Jolla, CA, USA). To perform the beta diversity analysis, one representative sequence per OTU was selected to generate a phylogenetic tree using ARB software (Ludwig et al., 2004). The resulting tree, together with the sequence abundance data, were used for a UniFrac analysis based on a weighted matrix of distance in QIIME: the principal coordinates analysis (PCoA) (Lozupone et al., 2007). A linear discriminant analysis (LDA) effect size (LEfSe, http://huttenhower.sph.harvard.edu/galaxy) method was used to find biologically relevant features, which underlines biological consistency, statistical significance, and effect relevance (Segata et al., 2011). The Kruskal–Wallis rank-sum test was used to detect taxa with significantly different abundances with an α value < 0.05. To estimate the effect size and ranking of each taxon, LDA was applied. The threshold of the logarithmic LDA score for discriminative features was >2.0. A high LDA score indicated a significantly higher abundance of certain taxa (Guerrero-Preston et al., 2016). A partial least squares discriminant analysis (PLS-DA) was used to discriminate groups by pregnancy (Pérez-Enciso and Tenenhaus, 2003). It was performed via the mixOmics package (version 6.0.0) using R software (version 3.2.0). Variable importance in projection (VIP) >1 was used to select key variables, which were the most important contributors to generation of the model (Chen et al., 2011). Co-occurrence networks of the 50 most abundant genera were analyzed using Mothur and visualized via Cytoscape (http://www.cytoscape.org/). Spearman's correlation coefficients were also calculated. Edges were set among the genera for which the ρ-value was >0.6 and significant (P < 0.01). Redundancy analysis (RDA) was performed using CANOCO for Windows 4.5 (Microcomputer Power, NY, USA).

The subjects were Chinese and lived in Shanghai for many years. The demographic and clinical characteristics are listed in Table 1. The pregnant women had higher scores for PLI, GI, and SBI, indicating that the pregnant women had relatively more severe gingival inflammation. However, the subjects were still periodontitis-free as the scores were within the reference range. There was a similar trend of variations in progesterone and oestradiol during pregnancy and the postpartum period (Table S1). The hormone levels of the 3rd trimester were the highest in the pregnant group, whereas there were no significant differences from the 1st visit to the 4th visit in the non-pregnant group.

Illumina MiSeq sequencing produced 3,334,433 raw sequences, and after pre-processing, 3,053,338 usable high-quality sequences with an average of 43,005 ± 15,560 sequences per sample remained in the dataset (Figure S1A). The average length of the sequence was 392 bp without the primers, ranging from 388 to 426 bp (Figure S1B). After removing the low-credibility OTUs, taxonomic assignment of the sequences resulted in the identification of a total of 554 OTUs in the supragingival microbiota. The 200 most abundant OTUs represented 98.76% of all sequences.

The community richness, as estimated by the Chao 1 index, was similar among the groups (Figure 1A). The Shannon index, an alpha diversity estimator, was higher in P3 than in P1, P2, P4, and NP, and there were a significant difference between P3 and NP per the Dunn's multiple comparisons test, which indicated higher bacterial diversity of the supragingival microbiota in P3 (Figure 1B). The data from the non-pregnant group were merged as a whole due to the lack of a significant difference after the calculation (P > 0.05). The eight rarefaction curves for P1–P4 and NP1–NP4 stabilized at >150 OTUs with 20,000 sequences, showing that the depth of sequencing was close to saturation and included most species found in the samples (Figure S2). Compared with the non-pregnant groups, P2 and P3 had a higher number of OTUs. These results were corroborated by the Venn diagram (Figure S3), which also revealed a higher number of OTUs in the pregnant group (mean = 553) than in the non-pregnant group (mean = 505). Specifically, the pregnant group had a higher number of OTUs at each given time point.

Structural similarity was explored with a PCoA of the beta diversity analysis (Figure 1C), which contained the first two principal coordinate axes (factors). Factor 1 accounted for 31.94% of the total variation, whereas factor 2 explained 16.15% of the variation. The clustering and ordination were correlated with biotic and abiotic variables that control the microbial community composition (Barberán et al., 2012). The plot of the non-pregnant group indicated similar microbial structure during four visits (Figure S4). Although the non-pregnant and P1, P2, and P3 samples partly overlapped, a clear segregation of the five bacterial communities could be seen. In addition, differences between the trimesters still existed when gestation was prolonged.

To characterize the bacterial distribution, the plaque microbiota of the pregnant and non-pregnant women was analyzed in terms of relative taxonomic abundance. A total of 17 phyla, 28 classes, 51 orders, 98 families, 172 genera, and 554 species were observed. Figure S5A shows the phylum structure of the oral microbiota for each participant during the four time intervals. Of the major phyla, Proteobacteria (27.80% of the total OTUs achieved from the pregnant group), Bacteroidetes (24.93%), Firmicutes (20.48%), Fusobacteria (13.48%), and Actinobacteria (11.95%) were the most abundant in the pregnant group, together accounting for 98.63% of the total sequences. In the non-pregnant group, the top five phyla were Firmicutes (29.76%), Proteobacteria (21.44%), Bacteroidetes (21.44%), Actinobacteria (13.65%), and Fusobacteria (13.13%). These phyla were present in all participants at all four sampling times, constituting the predominant supragingival phyla (Figure S5B). The top 50 genera are displayed in the heat map (Figure S5C). Of these, Streptococcus (11.60%) was the most abundant genus in the pregnant group, followed by Neisseria (11.57%), Capnocytophaga (9.03%), Prevotella (8.23%), Fusobacterium (7.36%), Leptotrichia (6.06%), Porphyromonas (4.27%), Corynebacterium (3.76%), Veillonella (3.22%), and Actinomyces (3.04%), together comprising 68.06% of the total OTUs. In the non-pregnant group, the top 10 genera were the same as in the pregnant group, and together, they comprised 72.15% of the total OTUs.

To identify the distinguishing phylotype that most likely explained the differences between the plaque samples from the pregnant and non-pregnant groups, LEfSe was performed. A circular cladogram was generated to show the differentially abundant taxa (Figure 2A), and LDA coupled with the effect size measurement was used to assess the effect size of the taxa (Figure 2B). Of the altered families in the pregnant group, Neisseriaceae and Porphyromonadaceae and Spirochaetaceae were significantly enriched, and there were other six significantly different families. Furthermore, the following species were more predominant in the supragingival plaques of the pregnant women than the non-pregnant women: several aerobic genera, including Neisseria and Janibacter; facultatively anaerobic genera such as Cloacibacterium; and obligate anaerobic genera including Porphyromonas, Treponema, and Syntrophomonas. The microbiota of non-pregnant women was enriched with genera mainly belonging to the Firmicutes phylum, namely, Streptococcus, Veillonella, Megasphaera, and Moryella. Among the genera over-represented in the non-pregnant group, Streptococcus exhibited the highest effect size and was closely followed by Veillonella. The effect size and ranking of the LDA score indicated that differential enrichment of bacteria owing to different host environment (LDA >2, P < 0.05).

PLS-DA was applied to identify key genera that were responsible for structural differentiation due to host variations. Parameter VIP was used to select variables with the most significant contribution to discriminating the community structure. Thirty-three genera with VIP >1 were identified as key genera for separating the supragingival microbiota between the pregnant and non-pregnant groups (Table 2). Twenty-five genera were associated with the bacterial community of pregnancy, all of which were distributed across the following phyla: Bacteroidetes (5 genera), Proteobacteria (9), Actinobacteria (5), Fusobacteria (1), Spirochaetes (1), and Firmicutes (4). The table also shows eight genera that were enriched in the healthy non-pregnant group: Moryella, Megasphaera, Veillonella, Parvimonas, Eikenella, Haemophilus, Bacteroides, and SHD-231; these genera were identified as health-related bacteria (Wu et al., 2013). In particular, four genera ranked as the top two genera in the non-pregnant group (Moryella and Megasphaera) and pregnant group (Porphyromonas and Neisseria).

The top 50 genera in the relative abundance ranking were selected from the two groups. Genera that met the threshold of Spearman's ρ > 0.6 and P < 0.01 are shown in the networks (Figures 3A,B). Compared with the non-pregnant group, the positive and negative relationships of the genera in the pregnant group were more simplified. One main cluster (>5 nodes and connected with intense lines) was evident in the samples from the pregnant women, while two main clusters were evident in the samples from the non-pregnant women. From a macroscopic perspective, it was noteworthy that the genera in the clusters indicated in gray circles all exhibited a positive relationship among the genera and mainly belonged to the Proteobacteria phylum. Prevotella was the most dominant genus in the co-occurrence network of the pregnant group. Its negative relationship with the other genera was decreased, and its positive relationship was strengthened, compared with the non-pregnant group. As the most abundant genera in the non-pregnant group, Streptococcus showed mutual restriction with Fusobacterium, Porphyromonas, Tannerella, Prevotella, and Parvimonas.

A two-step redundancy analysis (RDA) was performed to explore whether a correlation between specific genera and hormonal fluctuations existed. At the first step of the RDA, non-pregnant or pregnant women were used as the constrained explanatory variables, and the relative abundances of all OTUs were used as the response variables (Figure 4A). The Monte Carlo Permutation Procedure (MCPP) showed that the constrained ordination model by the non-pregnant or pregnant group was significant (P = 0.046). One hundred ninety-five OTUs were identified as key variables associated with pregnancy, which explained 1.82% of the constrained variance. In the second step of the RDA, sex hormones were used as the constrained explanatory variables, and the relative abundances of the OTUs identified in the first step of the RDA were used as the response variables (Figure 4B), thus generating 64 OTUs that were closely relevant to the sex hormones. Blue arrows pointing in the same direction as the sex hormones were predicted to have a positive correlation. The MCPP showed that the model was also significant (P = 0.003). The heat map displays the dynamic variation in the selected OTUs and sex hormones (Figure 4C). Fifty-three OTUs, which mainly belonged to genera Treponema (12 OTUs), Prevotella (4 OTUs), Cardiobacterium (3 OTUs), Aggregatibacter (2 OTUs), Leptotrichia (2 OTUs), other bacteria of the family Lachnospiraceae (3 OTUs), or other bacteria of the order Clostridiales (3 OTUs), were more abundant in the plaques of the pregnant women. In contrast, 11 OTUs were observed to have negative correlations with sex hormones, among which 5 OTUs belonged to Streptococcus and 2 OTUs belonged to Veillonella.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.