Introduction

Human pregnancy and parturition are complex biological processes that proceed through the coordinated action of many biochemical signaling networks. There are many pieces of evidence to suggest that normal pregnancy induces transient physiological, hormonal, and immunological changes in a highly controlled and coordinated manner1,2. The maternal immune system plays a crucial role in maintaining the balance to tolerate fetal allograft while preserving innate and adaptive immune mechanisms for protection against microbial challenges3,4. Disturbance of the delicate balance of the biological processes leads to adverse pregnancy outcomes, one being preterm birth. Immune dysregulation due to microbial pathogenesis is also increasingly appreciated in preterm birth and other pregnancy-related complications5,6.

Body fluids like blood, saliva, tears, sweat, urine, and cerebrospinal fluid are a source of putative biochemical markers that reflect various pathophysiological disorders. The saliva secreted from salivary glands and gingival crevice has some components that originate from the plasma with an overlap of almost 20–30% with plasma proteome, and the majority of these exhibit antimicrobial activity, transport, and enzymatic functions7,8. Human saliva, therefore, is a potential diagnostic fluid that can reflect many pathophysiological states of the body9. The decrease in salivary pH, calcium, glucose and an increase in phosphate levels throughout pregnancy favor gingival inflammation and oral infections, which could have some associations with adverse pregnancy outcomes10,11,12,13,14,15,16. Extensive plasma proteomics has been studied in pregnancy and pregnancy-related complications in an attempt to find potential biomarkers17,18,19. However, a few saliva specific proteomics studies have been undertaken for pregnancy. It has been shown that the level of salivary placental growth factor (PlGF) is significantly lower in preeclampsia conditions than normal pregnancy20. Similarly, the level of anti-inflammatory protein Annexin-1 in the saliva is shown to be elevated in pregnant women with gingivitis21. Recent 2D-gel based-proteomics analysis of saliva from pregnant women who deliver premature babies has identified a Metallothionein 2 (MT2A) protein as a potential marker for preterm birth prediction22.

The present study demonstrates longitudinal maternal saliva proteome changes in normal pregnancy for the first time. The SWATH-MS analysis quantifies 65 proteins that show temporal variations with gestational age. These proteins are distributed into two distinct clusters with a unique expression trajectory. Neutrophil degranulation, antimicrobial peptide function, regulation of TLR by an endogenous ligand, platelet function regulation, and glucose metabolism are the major pathways associated with these proteins. The MRM-MS analysis of 12 selected proteins confirms a similar expression pattern as observed in discovery analysis. The proteomics study of saliva in term birth pregnancy cannot be derived in adverse pregnancy outcomes. However, this information may provide useful background knowledge to target the specific molecular pathways in biomarker discovery.

Results

Clinical characteristics of the participants

Our 34 study participants who fulfilled the inclusion criteria for this study represent the semi-urban population of a northern state in India. The median age of the participants was 23 (interquartile range (IQR) 20, 23) years and 24 (IQR, 22, 25) years respectively for the discovery and the validation cohorts (Table 1). All subjects delivered at term and had a baby with appropriate body weight for age (Table 1). The median gestational age at delivery and the median birth weight for the study population is shown in Table 1. Other clinical and socioeconomic characteristics of the participants are described in Table 1.

Table 1 Clinical and demographic characteristics of the study population.

Differential proteomics of saliva with the progression of pregnancy in term delivery

The saliva-specific library was created by combining 1D SDS-PAGE followed by LC-MS/MS and bRP-C18-HPLC followed by LC-MS/MS for label-free quantitation of salivary proteins at different stages of pregnancy in SWATH-MS workflow (Fig. 1). MaxQuant analysis of combined LC-MS/MS library data yielded ~758 unique proteins in 1% FDR. The functional analysis of identified library proteins belonged to the immune system, metabolism, and signaling pathways (Supplementary Figure S1 and Table S1). The reproducibility in SWATH run between technical replicates (within triplicate) and within different time points (V1, V2, and V3) was uniform across all samples with average Pearson correlation value ~0.99 and ~0.95 respectively (Fig. 2A). The iRT peptides which were added in each sample, facilitated retention time (RT) normalization of all LC runs and ensured the specificity of fragment ion peak extraction23. However, to minimize the systematic variation all samples were normalized based on the median of total ion current (TIC) (Fig. 2B). The identification of precursor ions and protein groups within technical repeats was highly consistent in all samples, which confirms data reproducibility (Fig. 2C). The median coefficient of variance (CV) was 7, 7.1, and 6.9% over triplicate acquisition for all samples at V1, V2, and V3 time points, respectively (Fig. 2D). The data-independent acquisition (DIA) data resulted in the quantification of the average 2510 peptide precursors corresponding to ~446 protein groups with an average of 61% data set completeness across the time points in all 20 saliva samples. Sixty-five (15%; 65/446) proteins were modulated in abundance as a function of gestational age (p-value < 0.05; q-value < 0.1, Table 2).

Figure 1
figure 1

Schematic representation of the label-free quantitative discovery and MRM based targeted validation workflow. The label-free DIA approach with a saliva sample (N = 20) was taken to evaluate saliva proteome along with saliva-specific spectral library generation by DDA method. The DDA data were exported to MaxQuant for library generation and subsequently, each DIA file was analyzed with a prebuilt saliva spectral library in Spectronaut. An independent set of saliva samples (N = 14) was used to corroborate predicted central regulator proteins by MRM. More details can be found under Material and Methods.

Figure 2
figure 2

Quality assessment of label-free quantitative mass spectrometry data acquisition. (A) Scatter plot of peak intensities between the technical replicates of the salivary proteome of visit V1, V2, and V3, the Pearson correlation coefficient is indicated at the bottom of each plot. (B) The left side panel shows boxplots of precursor responses before normalization for each run at V1, V2, and V3 with technical repeats. Similarly, the right-side panel shows boxplots of the same precursor responses after normalization based on the median of total ion current (TIC). (C) Reproducibility of identification of the average precursors (green columns) and protein groups (pink columns) within the visit and technical repeats in all analyzed 20 saliva samples. Data represents mean ± SE of 20 independent data for each run at V1, V2, and V3 with technical repeats. (D) The dot plot shows the %CV distribution for all conditions, V1, V2, and V3. Data represent mean with SD of 20 independent data for each run at V1, V2, and V3.

Table 2 List of proteins modulates with the function of gestational age.

Differentially regulated proteins constitute two distinct clusters

Hierarchical with short time-series clustering of 65 modulated proteins resulted in two clusters derived from expression patterns across the gestational age (Fig. 3A). The clusters 1 and 2 were comprised of 47 and 18 proteins respectively. All proteins with their cluster possession and respective fold changes are shown in Table 2 and their expressional trajectory within a cluster is shown in Fig. 3B. The pathways enrichment by over-representation analysis of 65 proteins was performed. Overall 37 out of 65 proteins were significantly enriched in 15 biological pathways (p-value < 0.05, Table 3). The significant enrichment was observed within cluster1where 30 out of 37 proteins (81%) were from cluster 1 and rest 7 (19%) proteins were from cluster 2. Interestingly, 19 out of 37 (51%) were essentially enriched in neutrophil degranulation as the topmost pathway (Table 3). The rest of the 18 regulated proteins constituted 14 biological pathways which are listed in Table 3. The pathway analysis demonstrated that there was a significant overlap in biological processes between the two clusters. The proteins with their signature expression pattern within cluster 1 and 2 might support all the enriched pathways in a highly controlled and coordinated manner to maintain a healthy pregnancy.

Figure 3
figure 3

Hierarchical clustering and expressional trend within the cluster. A short time series was used for the construction of a hierarchical cluster. The circular plot of the dendrogram from the hierarchical clustering shows the presence of two clusters. The gene symbols in cluster 1 and cluster 2 are shown and are colored orange and blue respectively (A). The longitudinal profile of proteins across the pregnancy within two clusters 1 and 2. The line plot is showing protein abundance (log, base 2) on the y-axis and visit V1, V2, V3 respectively on the x-axis with each protein represented by a line (B). The top panel shows the proteins of cluster 1 (orange) and the bottom panel shows the proteins of cluster 2 (blue). The two clusters show the opposite trend with time.

Table 3 Enriched biological pathways linked with the proteins that are changed as a function of gestational age.

Protein-protein interaction (PPI) and network analysis identifies central regulators in pregnancy progression

Next, to investigate central regulatory proteins within saliva that hold the central network to maintain the normal pregnancy outcome, the PPI network was constructed with the expression profile of 65 proteins at V2 and V3 separately with respect to V1 (Fig. 4). The sub-network comprises of 129 nodes and 348 edges were formed with 55 seed proteins. Ten proteins were not included in this network. Interestingly, 67% (37/55) of proteins showed a downward trend (green circle) at V2, while 65% (36/55) demonstrated an upward trend (red circle) at V3 (Fig. 4A,B).

Figure 4
figure 4

Protein-protein interaction network of modulated proteins to assess central regulators in the maintenance of normal pregnancy. The upper panel shows a ratio in protein abundance between visits, V2/V1 (A), and the lower panel shows between V3/V1 (B). “Minimum network” of the modulated proteins is shown in a force atlas layout format. Nodes represent individual proteins and edges between two nodes represent interactions between them. The colors denote the expressions of nodes. The red and green indicate up- and down-regulation, respectively. The intensity of the colors signifies expression levels.

To identify central regulators within 65 proteins, network analysis of subnetwork was evaluated based on two parameters, degree centrality, and betweenness centrality. We identified 35 proteins whose degree centrality was at least 2 (Table 4). We considered these proteins as central regulators.

Table 4 Degree and betweenness centrality values of the central regulators.

Targeted validation of salivary proteins by MRM mass spectrometry

The discovery phase label-free quantitation yielded 65 altered proteins of which 35 were assigned for central regulators. Here, twelve proteins had been selected based on their pathway enrichment performance (Table 3) for targeted label-free scheduled MRM based quantification in a separate verification cohort (N = 14). The transition assessment and peak integration showed consistency in retention time and peak area across all the samples. The retention time for β-galactosidase and iRT precursors was found consistent at visit window V1, V2, and V3, and within technical repeats with average CV < 2.0 (Supplementary Figure S2). A comprehensive list of precursors/fragments is given in Supplementary Table S3. Twelve proteins were modulated in abundance as a function of gestational age (p-value< 0.05; q-value < 0.1, Supplementary Table S4). Besides the expression trajectory of 12 proteins except for SNAP23 in MRM assay was showing a similar pattern to our discovery phase results (Fig. 5 and Table 2).

Figure 5
figure 5

SWATH and MRM concurrence of temporal expression patterns of selected proteins. The x-axis of the figure represents the visits (V1, V2, and V3) and the y-axis represents a log2 of fold change with respect to V1. The mean of the proteins across samples at each visit is represented as a dot connected by lines and colored by the experimental technique.

Discussion

The current study on salivary proteomics identifies 65 proteins which change in expression with two distinctive expression patterns as a function of gestation age. The two distinct protein expression patterns deduce the interplay of relevant biological pathways. Our bioinformatics analyses have shown that the function of these processes is mediated through the 35 central regulatory proteins. Twelve of these have been validated through a separate subset of saliva samples in our cohort, in the targeted proteomics approach. The expression pattern of selected regulatory proteins in the MRM-based targeted data exhibits a similar expression pattern as observed in the SWATH-MS analysis. Our results show that the majority of differentially regulated proteins are associated with neutrophil degranulation, regulation of TLR by endogenous ligand, antimicrobial peptide function, platelet function regulation, and glucose metabolism. Here, we identified 35 seed proteins based on the degree and betweenness of interaction in the PPI network. Among them, 12 proteins have been selected in the MRM study as the majority belongs to the highest enriched pathway. Moreover, these proteins are also involved in more than one biological pathway except HSPA8 as indicated in Table 3. HSPA8 has been associated with the PPI network with the highest degree and betweenness (Table 4). We show that heat shock cognate 71 kDa protein (HSPA8/HSC70), protein S100-A8 (S100A8), protein S100-A9 (S100A9), cathepsin-D (CTSD), matrix metalloproteinase-9 (MMP9), fructose-bis-phosphate aldolase A (ALDOA), complement C3 (C3), lactotransferrin (LTF), myeloblastin (PRTN3), enolase-1 (Eno1) and galectin-3-binding protein (LGALS3BP) from cluster1, and synaptosomal-associated protein 23 (SNAP23) from cluster2 are the major 12 proteins which are associated with the above mentioned biological functions.

Heat stress cognate 70 is involved in the inflammatory signal pathways via extracellular interaction with TLR2/TLR424 besides its chaperone function. We have detected more than two folds of increased expression of HSPA8 at V3 compared to V2. Similarly, TLR4 also recognizes calcium-binding S100 family members like S100A8 and S100A925,26. In this study, we have not seen much change of S100A8/S100A9 expression in early pregnancy. However, more than three folds increased expression is observed at 26–28 weeks of gestation. It has been shown earlier that the high expression of S100A8 in serum in early pregnancy is associated with loss of pregnancy27. Protein S100A8, S100A9, and their heterodimer calprotectin (S100A8/S100A9) released by activated phagocytes have a role in inflammatory response by recruiting leukocytes and cytokine secretion28,29. A high abundance of calprotectin has been detected in amniotic fluid during intra amniotic infection, premature rupture of membranes and preterm labor30. Cathepsin-D promotes neutrophil apoptosis by activating caspase-831. The deficiency of CTSD amplifies and prolongs neutrophilic inflammation in vivo and this pathway is crucial for the resolution of innate immune responses. We observed more than a two-fold increase in the CTSD level at V3 time point compared to V2. The plasma CTSD level is significantly lower in the first trimester compared to non-pregnant women and the expression level increases in the third trimester32.

During the advancement of pregnancy, continuous uterine tissue remodeling is required to accommodate the fetus33. The enlarged uterus in pregnancy undergoes dramatic changes during labor and the post-partum period34. The degradation and reorganization of extracellular matrix components (ECM) are the ongoing processes and it has been shown earlier that cytokines are involved in the production of MMP9 in human myometrium34. MMP9 degrades the wide range of ECM components34,35. An increase in MMP9 has been implicated in vasodilation, placentation, and uterine expansion during normal pregnancy. The decreased expression of MMP9 at the later stages of pregnancy may be associated with reduced vasodilation, increased vasoconstriction, and hypertensive pregnancy36. In the present study, we observe more than 1.5-folds increased expression of MMP9 at V3 time point compared to V1.

The tightly regulated expression of complement C3 in early and mid-pregnancy confirms its activation in a coordinated manner. C3 protein shows persistent expression throughout pregnancy in our data. Normal human pregnancy is associated with complement activation in plasma to combat pathogen attack37. An antimicrobial protein, lactotransferrin (LTF) also shows lower level expression at V2 time point compared to V1. However, with the progression of pregnancy, LTF concentration increases up to 2.0 folds at V3 time point compared to V1. LTF plays an important role in cervicovaginal infection by reducing cytokines level in cervicovaginal fluid38. The neutrophil-derived myeloblastin (PRTN3) plays an important role in re-establishing vascular integrity after leukocyte transmigration. During thrombotic and inflammatory events myeloblastin also protects endothelial cells from protease-activated receptor-1 induced permeability change39. More than 1.5 folds of increased expression of PRTN3 at a later stage of pregnancy might be induced neutrophils to secrete vascular integrity modulators and helps neutrophil transmigration39. SNAP23 has been involved in the secretion of gelatinase-rich tertiary granules from neutrophils and platelet α-granule release40,41. We observed high expression (>2 folds) of SNAP23 at V2 as compared to V1. Recent mice data have also confirmed that the deletion of the SNAP23 gene results in pre-implantation embryonic lethality42. The galectin-3-binding protein (LGLS3BP) shows consistent lower expression throughout pregnancy in our data. A recent study has shown that Galectin-3 induces preterm birth in the mouse model consequent to dental infection43. However, the precise role of LGALS3BP with galectin-3 in pregnancy is not elucidated.

Our data indicate that four proteins from cluster 1 enriched primarily to gluconeogenesis and glucose metabolism. Increased expression of ALDOA, enolase (Alfa and gamma; ENO1, and ENO2), malate dehydrogenase (MDH1) supports the enhanced gluconeogenesis at late pregnancy44. Increased expression of these proteins from early to mid-stage of pregnancy has been observed in our analysis. This confirms that enhance metabolism is an essential factor for fetal growth and development.

An important limitation of this study is the absence of a pre-pregnancy sample, that may have provided more in-depth insights on how pregnancy progresses from its conception to delivery. Similarly, the presence of samples beyond 29 weeks of pregnancy would have provided useful information on the protein trajectories till delivery. The comparative study of term versus preterm saliva as well as the plasma will provide a systematic variation of proteins that may be considered as potential biomarkers in pregnancy-related complications. The biobank of GARBH-Inicohort45 hosts saliva, plasma and high vaginal fluids samples for proteomics analysis. The biomarker identification in pregnancy-related complications like spontaneous preterm birth is underway using a nested case-control design.

Salivary samples are non-invasive and more accessible than other body fluids making them a better choice for the development of biomarkers of adverse pregnancy outcomes. Several studies have described the use of saliva as a diagnostic tool for many diseases46,47,48,49,50, but the application of saliva in pregnancy study needs further research. In conclusion, the present longitudinal salivary proteomics study provides a system-wide protein expression in normal pregnancy. The results suggested that 65 proteins were altered in expression with gestational age with two distinct clusters. Further functional enrichment analysis elucidates the involvement of these proteins in pregnancy-related biological processes. To the best of our knowledge, this is the first report on the longitudinal changes in salivary proteomics across term pregnancy. We believe the findings of this study may be useful for future studies to find either biomarker or mechanistic understanding of adverse pregnancy outcomes.

Materials and Methods

Study population

A cohort of pregnant women was initiated in mid-2015 at the district hospital in Gurugram, Haryana, India with the primary mandate to generate a risk-prediction algorithm for preterm birth based on multidimensional risk factors assessed during pregnancy. The profile (detailed objectives and methodology) of this cohort, recognized as the GARBH-Ini cohort is provided in an earlier publication45. Women are enrolled within 20 weeks of gestation and are followed at 4–5 time points during the antenatal period until delivery and once within 6 months of post-partum. Various bio-specimens are collected and ultrasound scans are performed at defined intervals as per protocol. This includes maternal saliva collected at enrolment (V1), 18–20 weeks (V2), and 26–29 weeks (V3) of gestation. The collection of saliva at delivery was avoided as the participants wouldn’t fulfill the preparatory criteria for collection as detailed below. Details of the ongoing collaborative interdisciplinary program (GARBH-Ini) is published elsewhere45.

Ethical consideration

Institutional Ethics Committees (Human Research) of Regional Centre for Biotechnology (RCB), Faridabad; Institutional Human Ethics Committee of Translational Health Science and Technology Institute (THSTI); Institutional Ethics Committee of Gurugram Civil Hospital (GCH) and Institutional Ethics Committee of Safdarjung Hospital (SJH)approved this study45. The detail ethical approval and consents were described earlier45. Briefly, eligible women were enrolled in the cohort after they gave their written informed consent. A broad consent was taken for analysis of biospecimens for research use after appropriate de-identification. All methods were performed in accordance with the relevant guidelines and regulations.

Selection of participants

A small subgroup of participants who had been enrolled before 14 completed weeks of gestation during April-October of any year of enrollment, delivered at term (37–40 weeks of gestational age) with singleton vaginal delivery, and had all three time-point samples were randomly selected for our stated objectives from the ongoing GARBH-Ini cohort. Those with pregnancy-related complications (such as pregnancy-induced hypertension, preeclampsia, eclampsia), associated co-morbidity (e.g. infective or metabolic) before or at any time during pregnancy, congenital anomalies, and multiple pregnancies were excluded. This ensured a homogenous group of women who had a normal term pregnancy. Of the 34 participants selected, 20 were considered as discovery cohort for SWATH-MS discovery proteomics study, and another 14 for verification of target proteins using multiple reaction monitoring (MRM) mass spectrometry.

Sample collection, processing, and storage

The maternal saliva was collected from the enrolled participants and was processed in the research laboratory established at the hospital using standardized operating protocols that have been harmonized with other such global cohorts45,51. In brief, a non-stimulated saliva sample was collected in the morning prior to eating in a plastic container and 100 X halt protease and phosphatase inhibitors were added to a final concentration 1X. Saliva was immediately centrifuged at 10,000 g at 4 °C for 30 min. The supernatants were filtered through a 0.45-micron syringe filter and stored at −78 °C in the biorepository established for the study at the Translational Health Science and Technology Institute prior to mass spectrometry analysis. All the samples were retrieved just before the mass spectrometry experiment.

Sample preparation for mass spectrometry

All samples were buffer exchanged five times with 100 mM ammonium bicarbonate (pH 8.0) through a 3 kDa MWCO membrane (Amicon Ultra, Millipore). Protein concentration was estimated with BCA protein assay kit and in-solution trypsin digestion was performed. Briefly, saliva samples were reduced with 5 mM DTT at 56 °C for 60 min and alkylated with 20 mM iodoacetamide in dark at 25 °C for 30 min. MS-grade trypsin (Pierce Thermo Fisher Scientific) was added with 20:1 protein to enzyme ratio and the digestion mixture was incubated at 37 °C for overnight. The pH of the solution was lowered to 2.0 with 1% formic acid to quench the activity of trypsin. The digested samples were evaporated to dryness and these samples were desalted with C18 fast-flow tips (Pierce Thermo Fisher Scientific). Desalted peptides were dried in vacuum and stored at −20 °C for further LC-MS/MS data acquisition. Before the acquisition, retention time calibration iRT peptides (Biognosys, Switzerland) were spiked to all the samples at a 1:10 ratio for retention time normalization52.

Spectral library generation in Data-Dependent Acquisition (DDA)

Saliva specific peptide spectral library was created with pooled samples (N = 20) from all-time points V1, V2, and V3 in a combination of 1D-SDS-PAGE followed by in-gel digestion and in-solution digestion followed by basic reverse phase chromatography strategy (Fig. 1). SDS-PAGE was performed with 20 µg of pooled saliva protein and subsequently, in-gel digestion was performed. Basic reverse-phase chromatography was performed in the 1260 Infinity HPLC system which was equipped with an autosampler and fraction collector. A total of 600 µg of resulting peptides from pooled saliva were loaded onto C18 column (Agilent, 300 extend-C18; 3.5 µm; 2.1 × 150 mm) at oven temperature 40 °C which was previously equilibrated with solvent A (solvent A: 10 mM ammonium formate, pH 10; solvent B: 10 mM ammonium formate in 90% acetonitrile, pH 10). The peptides were eluted with the gradient from 2–50% of solvent B over 65 min total run at a flow rate of 500 µl/min. The eluted fractions were pooled, vacuum dried and desalted with C18 tips (Pierce). All fractions were resuspended in solvent C (composition: 2% (v/v) acetonitrile, 0.1% (v/v) formic acid in water) and analyzed in Sciex 5600+ Triple-TOF mass spectrometer which was coupled with ChromXP reversed-phase 3 μm C18-CL trap column (350 μm × 0.5 mm, 120 Å, Eksigent) and nanoViper C18 separation column (75 μm × 250 mm, 3 μm, 100 Å; Acclaim Pep Map, Thermo Scientific) in Eksigent nanoLC (Ultra 2D plus) system. The peptides were separated using a 90-min acetonitrile gradient from 5–50% of solvent D (composition: 98% (v/v) acetonitrile, 0.1% (v/v) formic acid) at a flow rate of 300 nl/min with column temperature of 40 °C. Data of each fraction was acquired in DDA with positive ionization mode. A maximum of 25 precursor ions with charge state 2–5 which surpass 120 cps per cycle was selected for fragmentation and each MS/MS spectrum was accumulated in high sensitivity mode. Each cycle consisted of 250 and 80–100 ms acquisition time for MS1 (m/z 350–1250 Da) and MS/MS (140–1800 m/z) scans respectively with a total cycle time of ~2.3 to 2.8 seconds. Dynamic exclusion was employed for 8–10 sec. Each fraction was analyzed majorly in duplicate. A blank of solvent C was run between each type of sample to reduce the carryover. Mass spectrometer calibration was performed at the start of each sample using a β-galactosidase digest standard.

Data Independent Acquisition (DIA) for label-free peptide quantitation

To collect quantitative data each sample was injected in 3 replicates using identical LC conditions as DDA run. MS/MS data were acquired in SWATH mode and the instrument was configured as described by the previous study53. Briefly, in SWATH-MS mode, a set of 60 overlapping variable windows were constructed covering a mass range of 350–1250 m/z. The MS2 spectra were collected from 200–1800 m/z with the accumulation time of 60 ms in high sensitivity mode. Each SWATH‐MS cycle consisted of 100 ms of survey scan resulting in a total duty cycle of ~3.7 sec.

Targeted mass spectrometry of selected proteins

Liquid chromatography-multiple reaction monitoring (LC-MRM) assay was performed on a QTRAP 6500+ mass spectrometer (Sciex) coupled with the LC system (Nexera). Saliva digested peptides (8 µg) were mixed with β- galactosidase (78 fmol) and iRT peptides (1:10) prior to sample load onto a C18 RP column (Agilent, 300 extend-C18; 3.5 µm; 2.1 × 150 mm). The gradient of acetonitrile was ranging from 5–50% of the total 60 min run at a flow rate of 0.3 ml/min. The binary mobile solvent system was used as follows; solvent E comprised 0.1% formic acid in the water, and solvent F included 100% acetonitrile with 0.1% formic acid. Spiked peptides were used for routine assessment of instrument and chromatographic performance. Skyline (v.4.1), SRM-Atlas (http://www.srmatlas.org/), and our spectral library were employed to select precursors of most intense y and b ions (transitions)54,55. Two transitions per precursors were traced for quantitation (Summarized in Supplementary Table S3). MRM data were acquired in the positive ion mode with an ion source temperature of 550 °C with the spray voltage of 5500 V. Scheduled MRM transition was performed using a detection time window of 60 sec. All MRM data were analyzed by MultiQuant (version 3.0.3, Sciex). The MQ4 integration algorithm within MultiQuant was utilized for peak integration. All MRM peaks were inspected manually to ensure correct peak detection and accurate integration. Two technical repeats were performed for each of the individual samples. The raw MRM peak area of two transitions was summed for each precursor. The average peak area values of each precursor from two technical repeats were transformed to log base 2, and demonstrated as a response of respective proteins with gestational age.

Data analysis and availability

The IDA spectra were analyzed in Andromeda within the MaxQuant analysis software (version 1.5.4.1) with default settings56,57. IDA Search parameters were shown in Supplementary Table S2. Datasets were searched against the UniProt human reference protein database58(downloaded on July 2017) appended with iRT peptides fasta database (https://biognosys.com/). A false discovery rate (FDR) of 1% was imposed for peptide-spectrum matches (PSMs) and the target-decoy approach for protein identification. The saliva-specific ‘merged’ library (Supplementary Table S2) was built from the MaxQuant analysis file using Spectronaut (version 11) with inbuilt ID Picker algorithm59,60. All DIA datasets (.wiff files) were converted to HTRMS format and calibrated with retention time dimension using the saliva-specific ‘merged’ spectral library. Data were analyzed with Spectronaut 11 by default settings61,62. Both precursor and protein level Q value was set to 0.01. Global normalization was applied to correct a systematic variance in the LC-MS performance63. Spectronaut uses all precursors for the statistical test. However, only topN values are selected for the quantification. The protein quantitation and peptide quantities were calculated as a mean intensity within the XIC peak area of the respective fragment ions at MS2. The protein CVs were calculated based on the summed intensities of their respective peptides. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier ‘PXD014800’.

Statistical analysis

Clinical characteristics of the enrolled study group were summarized as median (IQR) for continuous variables, and percentages for categorical variables using Stata 15.1 software (StataCorp LP, College Station, TX, US). Data quality assessment for DIA data between technical replicates and within visit window V1, V2, and V3 peak intensity throughout LC run were log base2 transformed and Pearson correlation coefficient was computed using Graph Pad Prism (version 7). The statistical analysis for longitudinal changes of proteins with the function of Period of gestation (POG) was performed in R statistical language and environment (www.r-project.org)64. Protein abundance from the DIA dataset of all conditions (V1, V2, and V3) with 3 replicates for 20 samples were taken from the output file of Spectronaut and transformed to log base2 to improve normality. The proteins with unique identities were retained. The proteins which shared peptides and the peptides identified as contaminants were omitted from further analyses. Linear mixed-effects models were used to model the abundance of proteins as a function of gestational age using the lme4 package65. The statistical significance of associations was estimated by comparing the models with and without gestational age by performing likelihood ratio tests using ANOVA function. Similarly, for MRM data mean peak area of each precursor was transformed to log base2 and linear mixed-effects models were applied. The changes in protein abundance across the pregnancy were considered significant if FDR adjusted p-value (q-value) was less than 0.1.

Clustering, pathway enrichment, and protein-protein interaction (PPI) network analysis

The proteins that showed a significant change as a function of POG by the linear mixed-effects models were taken forward for clustering analysis. As our longitudinal data had short time series and unevenly spaced samples we used short time-series distance, a method specifically designed to address these challenges followed by hierarchical clustering66,67. Plots to show the pattern of longitudinal changes in protein abundance were shown for visit windows 1, 2 and 3. Reactome database was queried for the pathway enrichment analysis using an R-package, ReactomePA68. The protein-protein interaction (PPI) network was built by NetworkAnalyst with the IMEX interactome database69. The node parameters degree and betweenness were considered for the identification of central regulators and were determined through NetworkAnalyst. “Minimum network” of the modulated proteins was constructed in a force atlas layout format. Subnetwork with at least 3 nodes was constructed. The expression of the nodes was symbolized by their colors.