Analysis of the oral microbiome during hormonal cycle and its alterations in menopausal women: the “AMICA” project

The maintenance of human health is dependent on a symbiotic relationship between humans and associated bacteria. The diversity and abundance of each habitat’s signature microbes vary widely among body areas and among them the oral microbiome plays a key role. Significant changes in the oral cavity, predominantly at salivary and periodontal level, have been associated with changes in estrogen levels. However, whether the oral microbiome is affected by hormonal level alterations is understudied. Hence the main objective pursued by AMICA project was to characterize the oral microbiome (saliva) in healthy women through: profiling studies using "omics" technologies (NMR-based metabolomics, targeted lipidomics by LC–MS, metagenomics by NGS); SinglePlex ELISA assays; glycosidase activity analyses and bioinformatic analysis. For this purpose, thirty-nine medically healthy women aged 26–77 years (19 with menstrual cycle and 20 in menopause) were recruited. Participants completed questionnaires assessing detailed medical and medication history and demographic characteristics. Plasmatic and salivary levels of sexual hormones were assessed (FSH, estradiol, LH and progesteron) at day 3 and 14 for women with menstrual cycle and only once for women in menopause. Salivary microbiome composition was assessed through meta-taxonomic 16S sequencing and overall, the salivary microbiome of most women remained relatively stable throughout the menstrual cycle and in menopause. Targeted lipidomics and untargeted metabolomics profiling were assessed through the use of LC–MS and NMR spectroscopy technologies, respectively and significant changes in terms of metabolites were identified in saliva of post-menopausal women in comparison to cycle. Moreover, glycosyl hydrolase activities were screened and showed that the β-D-hexosaminidase activity was the most present among those analyzed. Although this study has not identified significant alterations in the composition of the oral microbiome, multiomics analysis have revealed a strong correlation between 2-AG and α-mannosidase. In conclusion, the use of a multidisciplinary approach to investigate the oral microbiome of healthy women provided some indication about microbiome-derived predictive biomarkers that could be used in the future for developing new strategies to help to re-establish the correct hormonal balance in post-menopausal women.


Results
Plasmatic and salivary hormones. Plasmatic FSH, LH, E2 and PGN were measured for women with menstrual cycle at 3rd (M-3rd day) and 14th day (M-14th day) and at one time point for post-menopausal women (MP) (Fig. 1A). As expected, E2 increased significantly in the M-14th day group as compared with M-3rd day, due to the ovulation stage, whereas the other hormones remained stable between these two groups. Moreover, FSH levels increased significantly in the MP group as compared to both M-3rd and M-14th days, whereas E2 decreased significantly in comparison with the same groups. Levels of salivary FSH is reported. Very interestingly, FSH was detectable only in MP samples and not in saliva coming from women with menstrual cycle (Fig. 1B). On the other hand, levels of other hormones did not change between the three groups (Supplementary Table S1).
Microbiota analysis. The 16S meta-taxonomic analysis identified 72 families and 76 genera, with relative abundances varying from less than 0.1 to 85 percent; discrepancy between number of families and number of genera identified reflects on one side the belonging of several genera to the same family, on the other side the impossibility to assign some reads to a specific genus.
In particular, the genera most abundant (present in % > 1.0) were the following: Streptococcus, Prevotella, Neisseria, Haemophilus, Veillonella, Granulicatella, Gemella, Porphyromonas, Fusobacterium, Actinomyces, and Rothia ( Fig. 2A). Genera identified with relative abundance below 1% were not equally distributed in all samples. Salivary FSH levels (not detectable, nd). * and ° indicate significant differences compared to MP and M 14th day, respectively. P < 0.05 was considered statistically significant. Two-way ANOVA, followed by Tukey post-hoc.  Table S2 (families) and Table S3 (genera) as mean ± SEM of the relative abundance. Raw counts for meta-taxonomic 16S microbiome profiling at the genera rank is reported as Supplementary excel file (Table S4). In Fig. 2B-C families and species are shown. Alpha (Shannon) and beta-diversity (Bray-Curtis distance) are reported in Fig. 2D-I. An OPLS-DA was applied to meta-taxonomic analysis and reported in Figure S1. The supervised regression yielded 0 components that is no latent components to drive class separation according to microbiota content expressed in the experimental conditions (the 3rd, the 14th fertile subjects and the menopause class). On the other side, unsupervised PCA produced a 5 components statistical model (R2 = 0.92, Q2 = 0.4) with no distinct cluster for the 3 groups. As we can observe from the scores plot in Figure S1A , there are some outliers and other samples which fall very far from the origin of the model plane; those microbiota content often belongs to the same subjects, regardless to the menstrual cycle phase (as for example M15, M16, M11, M17). Moreover, by exploring the associate loadings plot in Figure S1B we can deduce that the dominant families responsible for those data distribution are Prevotella and Streptococcus. The composition of salivary microbiome in terms of both families and genera remained relatively stable throughout the menstrual cycle and menopause. Statistical analysis at the rank of species did not show significance of any taxa after FDR correction. However, two taxa could be considered to be significantly different between experimental groups before correction ( Figure S2). Indeed, Prevotella copri was significantly higher in the MP group compared to M-14th group, and Veillonella tobetsuensis was significantly lower in MP group compared to the two other groups.
As shown in Fig. 3 almost all glycosyl hydrolase activities screened were higher in MP group (M-3rd or M-14th day vs. MP p < 0,0001). However, α-glucosidase and β-N-acetyl-glucosamminidase were higher in M-14th day group as compared to M-3rd day and MP samples (p < 0,0001). The association between GH activities and hormone levels were evaluated through Pearson's correlation analysis but no significance was observed ( Figure S3).
Untargeted NMR-based metabolomic analysis. To reveal different metabolic patterns characterizing saliva alteration during the physiological feminine life cycle from the fertile to the non-fertile period, NMRbased metabolomic analysis was conducted on saliva samples. The OPLS-DA regression resulted in a non-overfitted model (R2 = 0.51, Q2 = 0.28) with one parallel and 1 orthogonal component. The scores plot in Fig. 4A showed a clear class discrimination along the predictive component of all fertile saliva (at t [1] negative values) collected from both women at their 3rd day (M-3rd, purple squares) or the 14th day (M-14th,red squares) of Figure 3. Screening of Glycosyde Hydrolases Activities in saliva. * and ° indicate significant differences compared to MP and M 14th day, respectively. P < 0.05 was considered statistically significant. Two-way ANOVA, followed by Tukey post-hoc. The correlation analysis between salivary metabolites and plasmatic hormones ( Figure S4) showed no significant association.
Bin integrations of statistically significant metabolites are reported as bars plot in Fig. 5A-F. To further associate those metabolic alterations with the glycosidase activities measured from each sample, we integrate metabolites and enzymes expressions through hierarchical analysis. The associated heatmap (Fig. 6) revealed weak positive correlations between putrescine with α-D-glucopyranoside and α-D-mannopyranoside activities, while  www.nature.com/scientificreports/ negative correlations were found between both the short-chain fatty acids 3-hydroxyisovalerate and n-butyrate with α-and β-D-mannosidase. Moreover, to investigate metabolites variation among the same subjects when analyzing salivary samples collected at different time points of the menstrual cycle (3rd and 14th day), we adopted a multilevel approach able to reveal metabolomic changes excluding individual factors, so to emphasize only those molecular variations depending on the hormonal status ( Figure S5). Multilevel PLS-DA was applied to the reduced spectral dataset using mixOmics routine in R 19 . A 2 component regression resulted in a total of 18% of data variation explained without a discrimination of samples according to the menstrual cycle condition. In particular, the projection of all samples in the scores plot ( Figure S5A) is contained within the confident ellipses delimiting the two classes (confidence level set to 95%). The variables mainly characterizing the total saliva distribution at high values along the component1 in the correlation circle plot ( Figure S1B) are glucose (5.23, 3.75 and 3.77), acetate (1.93) and valine (2.25). On the other side, variables characterizing the distribution at negative coordinates are mostly related to polyamine (1.83-1.87), aspartate (2.85-2.87). None of those metabolites were found discriminant for the whole dataset analysis.
Moreover, to investigate the impact of hormonal alteration on female metabolic pathways and networks, we applied enrichment metabolic analysis, which provides mechanistic insight into the underlying biology of differentially expressed metabolites between fertile and non-fertile women. The algorithm found a network with 150 nodes with a threshold of p < 0.01, which are reported in Figure S6 and Table S5 NMR of the Supporting Information. In particular, the pathways enriched are the following: • Citrate cycle (TCA cycle); • Galactose metabolism; • Valine, leucine and isoleucine degradation; • Arginine and proline metabolism; • Taurine and hypotaurine metabolism; • Glycosphingolipid biosynthesis-ganglio series; • Lysosome; • Taste transduction; • Glucagon signaling pathway; • Type I diabetes mellitus; • Carbohydrate digestion and absorption; • Protein digestion and absorption; • Central carbon metabolism in cancer.

Targeted lipidomics analysis. Since increasing evidences in the last decade have shown a connection
between the endocannabinoid (eCB) system and the microbiome, especially in the gut and in the gut-brain axis, as reviewed elsewhere [20][21][22][23] , and it is well known the important role of the eCBs in the reproduction and in the interplay with estrogen system 24 , a targeted lipidomics approach was applied to identify and quantify the main eCBs and N-acylethanolamine related in saliva samples. As shown in Fig  www.nature.com/scientificreports/ Moreover, we analyzed the correlation between hormone levels and eCBs ( Figure S7) but Pearson's correlation analysis did not show any significant association.

Multi-omics analysis.
To better understand the variables associated with menopause, we integrated the results from meta-taxonomic analysis, glycoside hydrolases activities, NMR metabolomics and targeted lipidomics through mass spectrometry with estradiol levels, and anthropometric data using multiple factor analysis (MFA). We observed separation of the MP group from the two other groups (Fig. 8A). As expected, the main variables associated with menopause were age and anthropometrics (Fig. 8B). Aside from four microbiome genera, namely Bulleidia, Slackia, Scardovia and Alloscardovia, the omic and enzymatic approaches investigated in this study did not permit to separate the experimental groups. Indeed, when performing the analysis using only meta-taxonomic, NMR metabolomics, targeted lipidomics and glycoside hydrolases activities, it was not possible to separate the experimental groups, suggesting that these variables measured in saliva are not key players in these processes ( Fig. 8C-D). Nonetheless, there seemed to be potential associations between the variables we measured. We elaborated a correlation network between all the variables generated in this study (Fig. 9). This network points to many potentially central players in the saliva ecosystem. For instance, several metabolites gravitate glycoside hydrolases, including positive and negative correlations. The same is true for bacterial genera Abiotrophica, Serratia, Peptostreptococcus and Alloscardovia. We also observed two major clusters of bacteria. In addition to these general observations, the endocannabinoid 2-AG was strongly correlated with α-mannosidase. Since we observed an uneven distribution of smokers between our experimental groups, we validated if our conclusions could be affected by this lifestyle variable. Multiple factor analysis did not show differential clustering of smokers and non-smokers ( Figure S8). In addition, analysis of variance do not show significant different microbiome taxa or metabolites either for smoking status or interaction between smoking and experimental groups. Therefore, we consider that smoking is not a confounding factor in this study.

Discussion
In our study, we have observed no significant change in the oral microbiota composition between women in menopause and women with menstrual cycle. However, we report several significant findings using a multidisciplinary approach. In particular, even if no major changes in the microbiota composition of women's saliva were observed, the genera most abundant are those ones known to be predominant in healthy individuals, such as Streptococcus, Neisseria, Porphyromonas, Prevotella, and Veillonella 25 . This is in line with the medical history of our volunteers that reported no predominant pathological issues in the general and periodontal health. Moreover, this result is comparable to a recent study in which salivary microbiome from a larger cohort of women with menstrual cycle was analysed 26 . In fact, they observed no significant differences in alpha-diversity or phase-specific clustering of the overall microbiome, even though they measure the sugar intake by weekly dietary records and it appeared to influence the composition of the salivary microbiome during the menstrual cycle 26 . Blood analysis (Fig. 1A) confirmed hormonal status typical of menopause and different stages during menstrual cycle. However, in saliva only FSH levels were detectable solely in MP group (Fig. 1B), assessing their postmenopausal status. In agreement with previous data reported in literature, the activity of β-N-Acetyl-glucosamminidase was found to be the highest in saliva samples, among the other glycoside hydrolases (GHs) screened 27 . In particular, we observed that α -glucosidase and β -N-acetyl-glucosamminidase were higher in M-14th day group as compared to M-3 rd day and MP samples, but no significant correlation with hormone levels were observed ( Figure S3). A possible explanation of this data is the correlation between hormone levels and salivary flow, that in our project has not been investigated. In fact, recently, it was suggested a correlation between the changes in oral health as well as systemic health due to menopause and hormonal changes 28 . In particular, E2 level positively correlated with unstimulated salivary flow indicating that a decrease of E2 may result in decrease of salivary flow, which could www.nature.com/scientificreports/ in turn cause various oral problems associated with menopause and a possible alteration of oral microbioma. Moreover, β-N-Acetyl-glucosamminidase and α-glucosidase activities have been reported as the predominant GH not only in secreted saliva but also in healthy human dental plaque collected overnight 29 . Previous studies showed that Prevotella intermedia, Treponema denticola, Porphyromonas gingivalis, and Streptococcus gordonii, as constituent of oral microbiome, have genes coding for the β-N-acetyl-hexosaminidase and expressed this activity, which in fact indicated the presence into the oral microbiome of these organisms. In fact, the dental plaque environment favors organisms that are able to metabolize complex glycoprotein and mucins and the detection of some GH activities in human saliva have been reported and have been considered as indicator of several parameters correlated to human health [30][31][32] . In addition, a very recent study has also analyzed different physicochemical parameters of saliva in menopausal women, in which they did not report GH activities, but they stated that differences in saliva properties observed in menopause can potentially affect the oral environment, possibly increasing the risk of some pathological changes in the oral cavity and consequently indicating the need to take special care of this group of female patients in order to help them maintain proper oral health 33 . Untargeted metabolomics analysis revealed a class differentiation between the menopausal group and menstrual cycle (with no difference between 3rd and 14th day). Interestingly, among metabolites that contributed to the discrimination, we found the citrate pathway and taste transduction module in which citrate plays a key role 34 . In particular, sour taste is mainly correlated with citric, malic and other acids (acetic, adipic, fumaric, lactic, succinic and tartaric) as well as the production of citrate ions 34 . The increased citrate levels in saliva of menopausal women could be responsible for the enhanced perception of sour taste. However, Saluja and co-workers have reported that both pregnant and postmenopausal women appeared to have a reduced perception of sucrose but no significant difference in the taste perception of citric acid 35 . Moreover, we also found the involvement of the galactose metabolism, where sugars contribute in taste functions, being stimulators for sweetness 36 . Fucose and galactose showed lower concentration in menopause women's saliva which can further contribute with alteration in taste. Indeed, gustatory function and taste perception may alter eating habits and in turn modify oral microbiome affecting oral health. In addition, a galactose consumption together with low activity of galactose-1-phosphate uridyl transferase have been related to ovarian senescence 37 and dysfunction 38 among women non suffering with galactosemia, which is a known condition for premature ovarian infertility (POI) 39,40 . Another interesting pathway involved is the one responsible for polyamines, such as putrescine and spermidine, that are known to be produced by the intestinal microbiota and regulate multiple biological processes 41 . Interestingly, putrescine deficiency has been reported as one of the causes of poor egg quality in an aged mouse model 42 . Putrescine is produced in peri-ovulatory ovaries and its supplementation reduced egg aneuploidy, improved embryo quality, and reduced miscarriage rates in aged mice 42 . Moreover, we have observed low butyrate levels in MP group (Fig. 5F) that might be in line with the estrogen deficiency and related to specific microbial species (Firmicutes).
In addition, we applied a targeted lipidomics approach on saliva samples in order to investigate the possible alteration of endocannabinoid tone in association with menstrual cycle and/or menopause and oral microbiome composition. In fact, in the last decade increasing evidence have reported the cross-talk between the eCB system and the associated endocannabinoidome (eCBome), an expanded pleiotropic system that play a key role in several physio-pathological conditions, with the gut microbiome, mainly focused on metabolic and obesity-related disorders suggesting that modulation of the eCBome is related with changes in the gut bacterial community and, on www.nature.com/scientificreports/ the other hand, the modification of the gut microbiota by using probiotics, antibiotics or germ-free mice affected eCB signaling 20,43 . Even though the presence of the eCBs in saliva have been already reported and correlated with obesity and feeding status 44 , the potential interplay with oral microbiome and menstrual cycle and/or menopause has never been investigated. Interestingly, our data showed significant increase in 2-AG levels (Fig. 7B) in saliva of menopausal women, and no significant changes in the N-acylethanolamines measured. This result is in agreement with previous studies in which it is shown that plasma 2-AG is associated with menopause 45 . On the other hand, another study found an association between the Fatty Acid Amide Hydrolase (FAAH) expression in adipose tissue and anandamide circulating levels in postmenopausal women in association with obesity 46 . However, the body mass index (BMI) mean in MP group in our study does not indicate an obesity status and this may reflect in no change in AEA and the other two N-acylethanolamines analyzed. Recent studies further implicate the role of gut microbiota in eCBome signaling, as commensal microbe Bacteroides produces eCB-like N-acyl amides 47,48 that are able to bind receptors activated in the eCBome 47 .
Multiple factor analysis integrating microbiome, lipidomics, glycoside hydrolases activities and metabolomics data from saliva showed some, but very little, difference between the three groups. These results suggest that saliva may not harbor key biomarkers associated with menopause.
Some studies reviewed by Segovia-Mendoza and Morales-Montor reported a direct correlation between estradiol and the release of β-hexosaminidase in mast cell line 49 . Moreover, children with type 1 diabetes showed an increased salivary β-hexosaminidase concentration which may be useful in the diagnosis 50 . In addition, β-D-hexosaminidase seems to play a key role in the biofilm formation in saliva of some oral pathogen 51 . It is possible that the β-hexosaminidase activity has effects not only on interactions with saliva-coated surfaces, but also interactions with host epithelial cells. These enzymes may also influence interactions with other bacteria in terms of release of glucose or galactose that may act as growth substrates for cohabiting bacteria.
Our data shows that although no significant alterations in the composition of the salivary microbiome was detected, potential associations between the variables measured exist, such as the enzymatic activity of some glycoside hydrolases and the presence of specific metabolites. Multi-omics analysis did not show any significant discrimination combining simultaneously all experimental data, but revealed some metabolites and enzymes as the most responsible for the distribution.
In conclusion, the use of a multidisciplinary approach allowed an exhaustive investigation of the oral microbiome of healthy women and provided some indication about microbiome-derived predictive biomarkers that will potentially help in finding novel strategies to establish the correct hormonal balance in post-menopausal women.

Methods
Volunteers. The study was conducted as per the guidelines of Helsinki declaration and received the ethical clearance from CNR Ethical Commission (Prot. n. 3419). All participants provided written informed consent before enrolment in the study sample. Thirty-nine medically healthy women aged 26-77 years (19 with menstrual cycle and 20 in menopause) were recruited. Participants completed questionnaires assessing detailed medical and medication history and demographic characteristics (Table 1). No financial incentive was offered to any of the survey participants.
In brief a total of 39 healthy women were enrolled in the study: 19 with menstrual cycle and 20 post-menopausal women.
Inclusion criteria • Females (18-80 years), • To be able to understand and communicate in Italian, • To be able to give informed consent

Exclusion criteria
• Patients who underwent to dental care within the last 30 days • Patients who showed an irregular menstrual cycle or a cycle with a length more than 30 days, • Current infective disease or fever, nausea, diarrhea within the last 30 days, Next generation sequencing for salivary microbiota composition. Bacterial DNA was purified from salivary samples with the use of "QIAamp® DNA Microbiome Kit" which allows depletion of human host DNA from biological samples and yields enriched bacterial DNA. Purification protocol consists of a degradation of host DNA followed by a disruption of bacterial cells through a combination of mechanical and chemical lysis. Target DNA is then purified through adsorption to a silica membrane, washing steps and elution with appropriate buffer. Purified DNA was accurately quantified using a fluorescence-based quantification method, such as Qubit dsDNA HS and 5 ng DNA for each sample were used for amplification. Library production was carried out by the Ion 16S™ Metagenomics Kit, this kit contains primer pools for PCR amplification of 7 hypervariable regions of the 16S rDNA gene of bacteria. Amplicon libraries were then sequenced with the Ion PGM™ sequenc- In addition, OPLS-DA approach was also applied to taxa data matrix with observations in row and microbiota families UV (unit variance) normalized in columns. Such analysis was carried out to explore possible presence of dominant families according to subject's specific menstrual cycle time point (3rd day, 14th day of fertile women's menstrual cycle, menopause women).
NMR-based metabolomics. Saliva samples were thawed, centrifugated and 400 µL of salivary supernatant were mixed with 300 µL of Phosphate Buffer Saline (PBS, pH 7.4) with 10% D2O for lock procedure, and then transferred into an NMR tube. Final pH is on average 7.2-7.6. One-dimensional (1D-1H) spectra were recorded at 600.13 MHz on a Bruker Avance III-600 spectrometer equipped with a TCI CryoProbeTM fitted with a gradient along the Z-axis, at a probe temperature of 27 °C, using the excitation sculpting sequence for solvent H 2 O suppression 52 . Spectra were referred to internal 0.1 mM sodium trimethylsilylpropionate (TSP), assumed to resonate at δ = 0.00 ppm. For each sample, the metabolic profile was attained through the acquisition of proton spectra (1D-NMR) and 2D NMR mono and hetero nuclear experiments (TOCSY, HSQC), useful for the recognition and assignment of metabolites. Subsequently, all the one-dimensional spectra undergone to multivariate statistical data analysis with unsupervised and supervised regression methods in order to discriminate the classes of samples (fertile subjects / subjects in menopause) on the basis of the changes in their corresponding metabolic profile induced by the different hormone levels.
Multivariate data analysis. Saliva proton spectra ranging from 9.50 ppm to 0.50 ppm were automatically binned into 450 integrals of 0.02 ppm each using the AMIX 3.9.15 software package (Bruker Biospin GmbH, Rheinstetten, Germany). The residual water resonance region (5.05-4.55 ppm) was excluded, and each integrated region was normalized to the total spectrum area to avoid possible dilution effects on the signals. The obtained NMR data format, expressed by a matrix (X matrix), was firstly imported into mixOmics routine in R 19 , after UV scaling, to analyze with a multilevel PLS-DA approach the salivary metabolic variations of the reduced dataset only consisting of fertile subjects at the 3rd and 14th day of their menstrual cycle. Then, the whole X matrix was imported into SIMCA-P + 14 package (Umetrics, Umeå, Sweden) where Principal Components Analysis (PCA) and Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) were performed, after UV scaling. We first applied PCA-class to explore data trend and eventually exclude outliers (data not shown). Once homogeneity was assessed, we applied supervised OPLS-DA to investigate the metabolic differences in saliva collected in three diverse feminine hormonal conditions. Initially, PCA and PCA-class was used to reduce data dimensionality and to explore possible trends as well as the existence of outliers. Once class homogeneity was assessed for each group, supervised OPLS-DA was applied to emphasize categories discrimination, where dummy variables were assigned to define class belonging (Y matrix). Supervised regressions were conducted comparing saliva of fertile groups at different time of their cycle (3rd and 14th day) and menopausal class, in order to generate predictive models that better relate metabolites variation to the different hormonal assessments. Each model quality was evaluated by using the goodness-of-fit parameter (R2) and the goodnessof-prediction parameter (Q2) 53 together with an internal iterative 7-round cross-validation and permutation test (800 repeats). The discriminatory metabolites were selected based on the bins containing NMR signals which were not overlapping with other peaks for the relative quantification, performed using OriginPro 9.1 software package (OriginLab Corporation, Northampton, USA). Statistical significance for selected metabolites was determined by parametric (ANOVA with Tukey correction) or non-parametric (Mann-Whitney U) tests according to the results of normality test performed on data to evaluate each distribution (Shapiro-Wilk, Kolgomorov-Smirnov test). P values < 0.05 were considered as statistically significant. Moreover, spectroscopic data were integrated with the enzyme expressions evaluated for the glycosidase activity. For this purpose, a correlation map with hierarchical clustering was also generated by combining clinical test values and selected bin integrals of significant metabolites with R software 54 . The Euclidean distance was considered for the metrics and the centroid method for clustering criterion.
Network analysis. Enrichment analysis on selected and more representative metabolites found in fertile/menopause class separation was applied using the 'pagerank' method computed by the FELLA package in R 55,56 . Starting from the set of altered compounds, such analysis suggests affected reactions, enzymes, modules and pathways using label propagation in a knowledge model network based on Homo sapiens database in KEGG [57][58][59] . The resulting network and subnetwork are visualized and exported in related plot and table with a threshold of p < 0.001. Targeted lipidomics. Saliva samples collected according to the protocol reported above were subjected to a series of biochemical steps for the extraction, purification and quantification of endocannabinoids, as reported in the literature 18 . The extracts or fractions enriched in endocannabinoids were then analyzed by liquid chromatography coupled with mass spectrometry by chemical ionization at atmospheric pressure (LC-APCI-MS) as reported above 19 Screening of glycosyde hydrolases activities in human female saliva. Samples were thawed, centrifuged and the supernatants were diluted tenfold with Tris buffer 100 mM, pH 7,5 27 . These salivary supernatants were analyzed for the presence of glycosidase activities: GH activities were measured separately for each saliva sample using black 96-well microtiter plates (BD Bioscience). Aliquots (50 µL) of diluted saliva samples were added to 20 µL of 4-nitrophenyl (PNP) substrates solution in Tris buffer 100 mM, pH 7,5 to a final concentration of 4 µM. Plates were incubated at 37 °C for 1 h, before the reaction was quenched with 130 µL of 500 mM Na 2 CO 3 -NaHCO 3 buffer (pH 10.3).
Statistical analysis. Data were represented as mean ± SEM. Plasmatic hormones and glycoside hydrolases were analyzed using the Two-way ANOVA followed by Tukey post hoc test. eCB and related N-acylethanolamine levels were analyzed using one-way ANOVA, followed by Tukey post hoc comparisons.
Statistical analysis of 16S microbiome profiling was performed at the rank of genera using linear model followed by one-way ANOVA on rank-transformed relative abundances. Analyses were performed with and without 20,000 reads rarefaction with similar results. Taxa were considered significantly modulated if FDR-corrected p-values were lower than 0.05.

Multi-omics analysis.
Integrative analysis of meta-taxonomic analysis, glycoside hydrolases activities, NMR metabolomics and targeted lipidomics through mass spectrometry with E2 levels, and anthropometric data was performed using multiple factor analysis (MFA) from the FactoMiner package 60 . Variables were scaled to unit variance for MFA analysis. Prior to the construction of the correlation network, we scaled, centered and applied Yeo-Johnson transformation to reduce biases in the correlations. Spearman correlations were performed in Rstudio 1.3.1093. P-values were FDR corrected and only correlations with FDR-corrected p < 0.001 were included in the network. The network was build using Cytoscape 3.8.0. Node position was manually adjusted to improve visualization.
Ethics approval and consent to participate. The study was conducted as per the guidelines of Helsinki declaration and received the ethical clearance from CNR Ethical Commission (Prot. n. 3419).

Data availability
All data generated or analysed during this study are included in this published article and its supplementary information files.