Fenugreek Stimulates the Expression of Genes Involved in Milk Synthesis and Milk Flow through Modulation of Insulin/GH/IGF-1 Axis and Oxytocin Secretion

We previously demonstrated galactagogue effect of fenugreek in a rat model of lactation challenge, foreshadowing its use in women’s breastfeeding management. To assess longitudinal molecular mechanisms involved in milk synthesis/secretion in dams submitted to fenugreek supplementation, inguinal mammary, pituitary glands and plasma were isolated in forty-three rats nursing large 12 pups-litters and assigned to either a control (CTL) or a fenugreek-supplemented (FEN) diet during lactation. RT-PCR were performed at days 12 and 18 of lactation (L12 and L18) and the first day of involution (Inv1) to measure the relative expression of genes related to both milk synthesis and its regulation in the mammary gland and lactogenic hormones in the pituitary gland. Plasma hormone concentrations were measured by ELISA. FEN diet induced 2- to 3-times higher fold change in relative expression of several genes related to macronutrient synthesis (Fasn, Acaca, Fabp3, B4galt1, Lalba and Csn2) and energy metabolism (Cpt1a, Acads) and in IGF-1 receptor in mammary gland, mainly at L12. Pituitary oxytocin expression and plasma insulin concentration (+77.1%) were also significantly increased. Altogether, these findings suggest fenugreek might extend duration of peak milk synthesis through modulation of the insulin/GH/IGF-1 axis and increase milk ejection by activation of oxytocin secretion.


Introduction
The World Health Organization (WHO) recommends exclusive breastfeeding for infants up to 6 months of age, based on evidence of its clear health benefits on mother-infant dyad [1][2][3]. Yet, 6 months after delivery, less than 40% of mothers are still breastfeeding in several high-income countries of North America and Europe [4]. Although the early cessation of breastfeeding is multifactorial [5], the perception that their milk secretion is insufficient to support adequate infant growth is reported by about 35% of lactating women [6,7]. Though perceived milk insufficiency due to psychological placed in a room with a fixed 12 h light-dark cycle (light from 7:00 a.m. to 7:00 p.m.). Pregnant rats had access to water and food ad libitum.
During gestation, every dam received a control diet (CTL) based on AIN-93G diet with 20 g protein per 100 g of food. At parturition (G21 = L0), dams received ad libitum either the control diet (CTL) alone or the CTL diet supplemented with dry water extract of fenugreek seeds (Plantex, Sainte-Geneviève-des-Bois, France) (FEN) at 1 g/kg body weight/day as previously described [15]. The diets were manufactured by the "Experimental Food Preparation Unit" (INRAE-UPAE, Jouy-en-Josas, France).

Experimental Design
The experimental protocol was approved by the Animal Ethics committee of Pays de La Loire and registered under reference (protocol APAFIS 2018121018129789) (Angers, France, 10 December 2018). All animal procedures were carried out in accordance with current institutional guidelines on animal experimentation in the EU (directive 2010/63/EU), in France (French Veterinary Department-A44276) and the French Ministry of Research in compliance with the commonly-accepted "3Rs".
The study consisted of a succession of 3 separate experiments (Xp1, Xp2 and Xp3) with 22, 12 and 11 individual litters respectively, as shown in Figure 1. In the three experiments, all dams received CTL diet during gestation. Delivery occurred at G21 also considered as day 0 of lactation (L0) and dams were randomly assigned either to the CTL or FEN diet throughout lactation. Pups were pooled by sex and randomly adopted by CTL or FEN dams. Litters were homogenized to 12 pups with a sex ratio of 1/1 and balanced weight of 84 ± 2 g ( Figure 1). The 3 experiments differed by the time at which animal follow up was stopped: post-natal day 21 for Xp1, L18 for Xp2 and L12 for Xp3. In Xp1 and Xp2, milk flow between dams and their pups was recorded between L11 and L18 using the deuterated water enrichment method and milk sample was collected at L18. In the Xp1, as previously described [14], pups were weaned at L20 to start the involution phase for dams' mammary glands. At post-natal day 21, also assimilated to first involution day (Inv1), dams were fasted for 4 h before sacrifice and sampling of plasma and tissues. In Xp2, after milk sampling at L18, food and pups were removed from the cage for 2.5 h and 1 h respectively, before the beginning of dams' sacrifice to avoid bias linked to hormone secretion during suckling. Finally, in Xp3, milk was sampled at L11 and dams were sacrificed at L12 after 1.5 h fasting and 1 h removal of pups ( Figure  1).  [14,15] between L11 and L18. Tissues sampled were mammary and pituitary glands.  [14,15] between L11 and L18. Tissues sampled were mammary and pituitary glands.

Dam-Litter Dyad Follow-Up and Milk Flow Measurement
Throughout lactation period, the dam weight, food and water consumption and total litter weight were recorded every two days from L0 to L11, and every day from L11 to L18. Pup weight was calculated as litter weight divided by litter size (12) and weight gain by subtracting weight al L0 from daily weight.
Milk flow was measured by the water turnover method as previously described [14,15]. Briefly, at L11, after 4%, isoflurane anesthesia, mothers received an injection of 4.92 ± 0.10 g/kg body weight D 2 O (99.9 mole% D 2 -enrichment, Sigma-Aldrich, Lyon, France) into a tail vein. Plasma samples were collected from dams by tail snip from L11 to L15 and pooled urine samples from their own whole litters were collected from L12 to L18 by stimulating lower belly with an iced cotton bud. The D 2 O enrichment of both plasma and urine samples was measured using the Fourier Transform infrared spectrophotometer α II ® (Brucker, Rheinstetten, Germany). Milk flow calculation was performed using a bi-compartmental model [15] as no significant intra-litter sexual dimorphism was observed in previous experiments [14]. Flow constants of the model (h −1 ) were then determined with the SAAM II ® software (The Epsilon Group, Charlottesville, Virginia, US) and the dams' total body water (g) using the intercept of dams D 2 O concentration curve. The water flow, from dams to their own litter (g/h), was calculated as the product of total body water by flow constant from dam to the litter and was expressed in g/day by multiplying by 24.

Milk, Plasma and Tissue Sampling
After 20 min dam/pups' separation and intraperitoneal injection of oxytocin (1 unit of Syntocinon ® ; Sigma-Tau, Ivry-sur-Seine, France) to stimulate milk ejection, dams were anaesthetized with 4% of isoflurane, and milk was manually collected and stored at −20 • C until analysis. Blood samples were collected before animal sacrifice via intra-cardiac puncture and placed into Ethylene Diamine Tetraacetic Acid (EDTA) tubes (Pfizer-Centravet, Plancoët, France), centrifuged at 1132 g for 15 min at 4 • C, and plasma were stored at −20 • C until analysis. Dams were sacrificed by intra-cardiac injection of 0.5 mL Exagon ® (Richter pharma, Wels, Austria). Left inguinal mammary gland and pituitary gland were removed, sampled and immediately frozen in liquid nitrogen before storage at −80 • C.

RNA Extraction and Purification
The entire pituitary gland (weight 19 ± 5 mg) or 28 ± 8 mg (mean ± standard deviation (SD)) of the inguinal mammary gland were placed into 2 mL Eppendorf tubes with glass beads in carbonic ice, and were then crushed using Precellys Evolution Homogenizer ® (Ozyme, Saint-Cyr l'Ecole, France) at 4000 rpm for 1 min. Total RNA was isolated from each tissue sample and purified following the NucleoSpin ® RNA/Protein (Marchery-Nagel, Düren, Germany) extraction/clean up protocol. RNA was eluted in 60 or 100 µL of RNase-free water. RNA concentration was measured at 260 nm and sample quality determined by ratios of 260/280 and 260/230 nm (Supplementary Table S1) using a NanoVue Plus ® (GE Healthcare, Little Chalfont, UK). As both ratios were in the range of 2.0-2.3, this guaranteed of sufficient pure RNA with low contamination by proteins and extraction solvents [31].

Reverse Transcription (RT)
Complementary DNA (cDNA) was prepared from 1.5 µg of purified total RNA using SuperScript ® IV reverse transcriptase (Invitrogen, Villebon-Sur-Yvette, France) according to the manufacturer's instructions. Moreover, some samples were randomly chosen to perform negative RT (RT -) by replacing SuperScript ® IV with water, to verify the absence of non-specific cDNA amplification in samples. Reverse transcription was performed with a SimplyAmp ® thermal cycler (Applied Biosystem, Villebon-Sur-Yvette, France) with thermal protocol of 23 • C for 10 min, 50 • C for 10 min and 80 • C for 10 min. Finally, the cDNA obtained were diluted to 1:40 before quantitative analysis.

Primers Design and Validation
A total of 54 genes linked to milk constituents secretion and its regulation and lactogenic hormones secretion were studied [19,32,33] (Appendix A, Table A1). Primer pairs used for qPCR were designed based on the sequence obtained from the online gene section of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/gene/) and using open-source PerlPrimer v1.1.21 software (http://perlprimer.sourceforge.net.) [34] according to the following parameters: primer melting temperature ranging from 58 to 62 • C with a maximal difference of 2 • C, amplicon size ranging from 100 to 150 bp, primer length ranging from 20 to 24 bp, span intron/exon boundary and overlap intron/exon boundary by 7 bases to avoid the possibility of genomic DNA contamination. The selected primers had, if possible, no expandable dimer and a total dimerization with the lowest possible force (dG < 3). Then, to check for their selectivity and the absence of potential genomic DNA amplification, selected primers were blasted using Primer-BLAST of NCBI (https: //www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi), against Rattus norvegicus "Refseq mRNA" and "Refseq representative genome". Finally, all the primers were validated by measuring their efficiency (E) to amplify mammary gland or pituitary gland cDNA and by ensuring they did not amplify RT -. Efficiency was measured by preparing a pool of all samples of 1:40 dilution of cDNA which was then serially diluted 5 times to 1:2. The curve of cycle quantification value (Cq) function of log 2 (dilution) was linearly fitted and PCR efficiency (E) was calculated according to the following formula: Sample cDNA amplifications were compared to the five RTamplifications by calculating ∆Cq (RT-) = mean Cq (RT-) -mean Cq (sample) , setting a Cq value of 50 assigned in case of no amplification (Supplementary Table S1). All genes with their official symbol (NCBI), main functions, primer sequences, amplicon length and efficiency are listed in Appendix A Table A1. Overall, the mean amplicon size was 124.7 ± 28.6 (mean ± SD). For mammary and pituitary glands, efficiencies were greater than 90 and ∆Cq (RT-) greater than 10, indicating good amplification capacity and negligible genomic DNA amplification.

Quantitative Polymerase Chain Reaction (qPCR) Analysis
The qPCR was executed in Hard-Shell ® 96-well PCR plates (Bio-Rad, Marne-la-Coquette, France) on a CFX connect TM Real-Time PCR detection instrument (Bio-Rad). Reactions were performed with a total volume of 15 µL, including 5 µL of 1:40 dilution of cDNA, 250 nM of forward and reverse gene-specific primers and iTAQ Universal SYBR ® Green Supermix (Bio-Rad). The qPCR conditions were as follows: initial denaturation at 95 • C for 5 min followed by 40 cycles of 95 • C for 10 s, 60 • C for 30 s and 72 • C for 10 s. To check the specificity of PCR products, a melting curve was performed by increasing temperature from 65 to 95 • C with a 0.5 • C increment and 10 s per step.
Relative gene expression was calculated using the 2 −∆∆Cq method as previously described [35]. For mammary gland, arithmetic mean Cq of three housekeeping genes (HKG): Actb, Uxt and Rps9, was used to calculate ∆Cq (Cq Gene of interest -Cq Housekeeping genes ). Actb is a classically used HKG while Uxt and Rps9 are known to be very stable genes in lactating mammary gland [36,37]. For pituitary gland, arithmetic mean Cq of Actb, Vcp and Rps9 was used as HKG to calculate ∆Cq. Vcp was shown to be very stable in the brain in an RNAseq study of the lab (unpublished data) and the Cq values of the three genes were strongly correlated (r from 0.68 to 0.81). ∆∆Cq values (∆Cq Target group -∆Cq Reference group ) were calculated using the CTL diet at L12 as a reference group to highlight the gene expression evolution over the lactation period. Values of −∆∆Cq were considered to be representative to the Log2(Fold change) assuming that mean primer efficiency was close to 100% and relatively similar between genes (CV of 4.7% and 5.9% in mammary and pituitary glands, respectively). Finally, fold change expression of the FEN diet was calculated for each lactation period with 2 −∆∆Cq formula using, for each day, the CTL diet as a reference group.

Statistical Analysis
To study the global impact of fenugreek on genes responsible for milk synthesis and regulation in the mammary gland, multivariate analysis was performed, using SIMCA P ® version 15.02 (Umetrics AB, Kinnelon, NJ, USA). The unsupervised Principal Component Analysis (PCA) and the supervised Partial Least Squares Discriminant Analysis (PLS-DA) were used [38] to cluster the mammary gland samples based on the selected gene expression profiles (−∆∆Cq values). PCA reduced the large set of variables (43 genes) into two principal components (PCA 1 and 2). PLS-DA enhanced the separation between the dietary groups (CTL and FEN) by summarizing the data into a few latent variables that maximized covariance between the response (dietary group) and the predictors (gene expression). The performance of the PLS-DA models was evaluated by the proportion of the total variance of the variables explained by the model (R2(X)), the proportion of the total variance of the response variable (dietary group) explained by the model (R2(Y)), the prediction ability of the model reflected by both a 7-fold cross-validation of the data (Q2) and the reliability of the PLS-DA model based on an analysis of variance (ANOVA) assessment of the cross-validatory (CV) predictive residuals (p-value of CV-ANOVA) of the model. The corresponding loading plot was used to determine the genes most responsible for separation in the PLS-DA score plot. Based on the PLS-DA results, genes were plotted according to their importance in separating the dietary groups in PLS-DA score plot and each gene received a value called the variable of importance in the projection (VIP). The variables with the highest VIP values (above 1) were the most powerful group of discriminators involved in the separation of CTL and FEN groups. Finally, in order to highlight diet-specific effect on mammary gland gene expression, lactation time effect was masked by building new PCA and PLS-DA models on a global matrix of normalized variables obtained using mean centering and standardization of variables by lactation point.
For univariate analysis, parametric tests were favored to maximize statistical power, such as the Student's test used for diet comparison. Their validity was admitted as the residual distribution of variables was not significantly different from normal distribution according to the Shapiro-Wilk test and variance was homogenous according to the Brown-Forsythe test. Relative genes' expressions in tissues and plasma hormone concentrations were analyzed using two way ANOVA with diet and period factors followed by Tukey's post-hoc test for period pairwise comparisons (3 levels) and Sidak's post-hoc test for diet pairwise comparisons (2 levels), as suggested by Kim [39]. Values of −∆∆Cq were chosen for statistical analysis of gene expression because distribution of −∆∆Cq values, unlike that of 2 −∆∆Cq values, was similar to normal distribution that enabled the use of parametric tests [40,41]. To determine the strength of the link between two variables, Pearson's correlation test was used. All univariate tests were performed with the GraphPad prism ® software (San Diego, CA, USA), version 6. In all tests, p < 0.05 was considered as statistically significant.

Fenugreek Supplementation Enhances Lactation Performances in a Model of Rat Nursing a Large Litter
Lactation performance was measured in experiments 1 and 2 (Xp1 and Xp2). As no significant interaction effect was found between experiments and diet effects using an ANOVA, data from Xp1 and Xp2 were pooled. Dam water intake was similar in control diet (CTL) and fenugreek-supplemented Genes 2020, 11, 1208 7 of 28 diet (FEN) groups throughout the lactation period. Conversely, maternal food intake, pup final weight gain and milk flow from mother to pups all increased significantly, by 10.0%, 7.1% and 15.3%, respectively (Table 1).

Different Patterns of Gene Expression in the Mammary and Pituitary Glands and of Plasma Hormone Concentration throughout Lactation and Involution
Multivariate and univariate statistical analysis were performed to investigate the temporal pattern in gene expression in mammary gland from mid-lactation to the beginning of involution.
The genes with the highest expression in mammary gland were Csn2, Wap, Mtco1 and Lalba with ∆Cq ≤ −6 and codes for major milk proteins and actors of mitochondrial oxidative phosphorylation (Tables 2 and 3). Plin2 and Fas, two genes involved in lipid synthesis, were also strongly expressed in the mammary gland with ∆Cq < −1. Conversely, Cpt1a, Acads, Mtor, Igf1r Esr1, Insr, Lxra, Scd1, Pgm1 and Dgat2 were the least expressed genes in mammary gland over the entire lactation period, with ∆Cq > 4 (Tables 2-5). In pituitary gland, Prl and Gh1 were both very intensively expressed, with ∆Cq < −6 highlighting the high level of secretion of prolactin and GH during lactation (Table 6). In contrast, Oxt, Trhr and Vipr2 were the less expressed genes with ∆Cq > 4.     Following the unsupervised Principal Component Analysis (PCA) on the un-normalized complete set of 43 genes, the PCA score plot, on components 1 and 2, revealed natural longitudinal classification with differences corresponding to the three lactation time points and a clear-cut separation between involution (Inv1) and lactation period (L12 and L18) (Figure 2a). The PCA loading scatter plot (Figure 2b) reflected relationships between variables in this model, deciphering three clusters that we labelled "kinetics 1, 2 and 3". In agreement with this longitudinal pattern, the univariate statistical analysis showed that expression of all genes evolves significantly between the middle and end of lactation, except for Ghr and Mtor (Tables 4 and 5).
Most of the genes studied exhibited a significant decrease in their expression between L18 and Inv1 ranging from 1.46-fold (p = 0.010) for Atp5f1a to 36.91-fold (p < 0.001) for Fasn and belonged to "kinetic 1" and "kinetic 2" patterns.
In the "kinetic 1" pattern, genes similarly expressed or only slightly modulated between L12 and L18, suggesting they had reached peak expression in this interval (Figure 2c). This pattern was represented mostly by genes coding for regulation of milk nutrient synthesis, energy metabolism, and galactose uptakes (Stat5, Akt1, Spot14, Pparg, Prlr, Pdha1, Acacb, Mtco1, Cs, Ugp2 and Atp5f1a). All these genes were highly correlated with each other, as validated by Pearson's test (Figure 2f and Supplementary Figure S1) and mainly with Prlr (r ≥ 0.8), considered as the barycenter of the "kinetic 1" cluster in the PCA-loading scatter plot (Figure 2b). If not in the "kinetic 1" ellipse in the PCA (Figure 2b) genes related to de novo fatty acid synthesis and glucose uptake such as Lpl, Fasn, Acaca, Fabp3 and Glut1 were also strongly correlated with Prlr (r ≥ 0.7) (Supplementary Figure S2). In "kinetic 2", overall gene expression rose from 1.70-fold (p = 0.042) for Csn2 to 2.62-fold for Aqp3 (p = 0.001) between L12 and L18, suggesting that peak expression of these genes occurred closer to L18 than L12 (Figure 2d). This cluster mainly included genes related to milk constituent secretion (Slc7a5, Wap, Csn2, B4galt1, Lalba, Plin 2, Aqp3, InsR, Esr and Dgat1) as illustrated in Figure 2b,d,g, with Insr, closest to the barycenter of the "kinetic 2" cluster ( Figure 2b), and strongly correlated with all these genes ( Figure 2g). Otherwise, Esr1 was probably linked to protein expression as its expression was very strongly correlated with all main milk protein genes (r ≥ 0.85).

Fenugreek Supplementation Impacts Mammary Gland Metabolic Pathways Differently during Lactation and at Involution
A supervised Partial Least Squares-Discriminate Analysis (PLS-DA), built on the same original dataset as PCA, showed that the stage of lactation had more impact than diet on the first two components (Supplementary Figure S1). Thus, a mean-centering and standardization by lactation stage of all variables was performed to assess the effect of diet per se. The resulting PLS-DA score plot revealed a clear-cut separation of diet groups (Figure 4b) overall. Most of the genes studied contributed to diet separation with a VIP over 0.8. Acacb, Fabp3, Plin2, Glut1, B4galt1, Acads, Cs, Scd1, Lalba, Csn2, Slc7a5, Mtco1 and Igf1r were the genes the most associated to the FEN diet while Sod1, Fabp4, Dgat2, Pparg, Prlr, C3 and Pgm1 were the most associated to the CTL diet (Figure 4a,c).

Fenugreek Supplementation Impacts Mammary Gland Metabolic Pathways Differently during Lactation and at Involution
A supervised Partial Least Squares-Discriminate Analysis (PLS-DA), built on the same original dataset as PCA, showed that the stage of lactation had more impact than diet on the first two components (Supplementary Figure S1). Thus, a mean-centering and standardization by lactation stage of all variables was performed to assess the effect of diet per se. The resulting PLS-DA score plot revealed a clear-cut separation of diet groups (Figure 4b Fenugreek modulated mammary gland gene expression mostly at L12 and Inv1. By focusing on each lactation time point with the PLS-DA, the predictive ability of the L12-( Figure 4e) and Inv1 (Figure 4h)-PLS-DA models was more important (R2(Y) cum and Q2 cum values above 0.9 and 0.7, respectively) than that of L18-PLS-DA model (Q2 cum = 0.39). At L12, the genes contributing the most to the separation of diets were Acacb, Fabp3, Fasn, Glut1, Pgm1, B4galt1, Slc7a5, associated with the FEN diet, but also Fabp4, Dgat2 and C3 associated to the CTL diet (VIP > 0.8, Figure 4d,f). Moreover, most of the genes overexpressed in the FEN diet at L12 were related to the synthesis of milk constituents, energy metabolism and the regulation of milk secretion (Tables 2-5). At Inv1, some genes contribute strongly to the separation of both diets such as Mtco1, Cs, Aqp1, Akt1, Cat, Plin2 and Csn2 associated with the FEN diet, but also Lpl, Dgat2, Pgm1, Ugp2, Pparg and Sod1 associated with the CTL diet (Figure 4g,i). Genes associated with the FEN diet were mostly related to milk secretion and energy metabolism.
At L18, diets were poorly separated from each other by PLS-DA and no difference was observed in mammary gene expression between CTL and FEN diets.

Fenugreek Supplementation Increases the Expression of Mammary Genes Involved in Lipid, Lactose and Protein Synthesis
The effect of fenugreek was studied individually on genes involved in the synthesis of milk constituents ( Table 2). The FEN diet significantly increased the overall expression of genes related to milk lipid synthesis (Lpl, Fas, Acaca, Fabp3), lactose synthesis (Glut1, Ugp2, B4galt1, Lalba) and protein synthesis (Lalba and Csn2) ( Table 2). Moreover, Acacb, Fabp3, Plin2, Glut1, B4galt1, Scd1, Lalba, Csn2 and Slc7a5 were strongly associated with the FEN diet in the overall PLS-DA (VIP > 1, Figure 4a,c) and highly predictive of FEN group (regression coefficient ≥ 0.5).

Fenugreek Supplementation Increases Expression of Genes Involved in Fuel Metabolism, Particularly Fatty Acid β-Oxidation
The overall expression of Cpt1a, Acads and Cs was significantly increased with the FEN diet compared to the CTL diet, with a particular overexpression of lipid β-oxidation-related genes, Cpt1a and Acads, at L12: 3.84-fold and 2.77-fold (p < 0.001), respectively (Table 3). Moreover, Acads, Cs and Mtco1 were among the most powerful discriminators associated with the FEN diet in the overall PLS-DA (VIP > 0.8, Figure 4a,c).

Fenugreek Increases the Expression of Regulatory Factors Involved in Protein and Lipid Metabolism
The FEN diet had no significant overall effect on the gene expression of factors regulating milk protein and lipid synthesis (Table 4) but induced an early transitory significant overexpression of Lxra and Mtor at L12 (2.38-fold, p = 0.01 and 1.84-fold, p = 0.032, respectively) contributing, with Pparg, to discrimination of the FEN diet group in the PLS-DA model at L12 (VIP > 0.8, Figure 4d,f) Additionally, Pparg and Akt1 were significantly upregulated at L12 compared to L18 (1.85-fold, p = 0.046 and 1.62-fold, p = 0.035, respectively) in the FEN diet, whereas similar expression was observed for the CTL diet (p ≥ 0.98).

Fenugreek Upregulates the Expression of Insulin, GH and IGF-1 Receptors
The effect of fenugreek was studied on receptors of the main lactogenic hormones (Table 5). Igf1r expression was the most enhanced by the FEN diet throughout lactation (p = 0.033; Figure 4a,c), specifically at L12 (2.60-fold; p = 0.012 and VIP > 0.8, Figure 4d), with a similar trend observed for Insr and Ghr expressions (VIP > 0.8, Figure 4d). Interestingly, normalized gene expression of Igf1r was strongly correlated with those of Ghr and Insr overall lactation (r = 0.76 and 0.71 respectively, Supplementary Figure S3).

Fenugreek Stimulates Oxytocin Expression at the Pituitary Level
The FEN diet increased the overall expression of Oxt, C3 and Trhr, mostly at L18 for Oxt and C3 (p = 0.029 and p = 0.004, respectively), and tended to increase the expression of Vipr2. In addition, Drd2 was overexpressed in the FEN group at Inv1 (p = 0.050), whereas it was on average under-expressed in this group during lactation (interaction: p = 0.040).

Fenugreek Supplementation Increases Plasma Insulin Concentration
The FEN diet had no significant overall impact on the concentration of these hormones ( Figure 4). However, plasma estrogen concentration tended to be higher in the FEN compared to CTL group (p = 0.082), particularly at Inv1 (+35.5%, p = 0.056). Moreover, variation of insulin concentration was not homogenous between lactation and involution (p = 0.035 for the Brown-Forsythe test), impairing our ability to detect any effect of supplementation during lactation. By centering and reducing insulin concentration by time period, a rise in insulin concentration was observed in the FEN group at L12 (+77.1%, p = 0.043), whereas no difference was observed at L18 and Inv1 (interaction: p = 0.040). Finally, a significant interaction was observed for prolactin concentration explained by the precipitous drop between L12 and L18 in the CTL group (94.1%, p < 0.001), whereas it remained steady in the FEN group (p = 0.995). However, the large prolactin concentration variations observed during lactation (mean CV of 102.1% and interval of (2.87; 528.2) ng/mL) makes it difficult to draw any conclusion from these values.

Discussion
The findings of the current study confirm that fenugreek increases milk production in a model of rat nursing a large litter. They further demonstrate that the impact of dietary fenugreek supplementation on milk production is associated with a stimulation of milk macronutrient synthesis and mammary fuel metabolism. These effects were mainly observed at L12 and Inv1 but not at L18, suggesting that fenugreek promotes an earlier peak of lactation and maintains milk production longer, whereas it does not seem to have further beneficial effect on milk macronutrient synthesis once peak lactation is reached. Our findings finally suggest that the effect of fenugreek on milk production at mid-lactation may be mediated by a stimulation of insulin secretion and a modulation of the insulin/GH/IGF-1 axis, while its action on maintaining lactation until the first day of mammary involution could be due to its estrogenic effect.

The Time Course of the Gene Expression in Mammary Gland and Lactogenic Hormones Suggests Peak Lactation Is Closer to L18
A noticeable result is the importance of lactation stage on mammary gland gene expression, as illustrated by the clear-cut separation of lactation time points in both PCA and PLS-DA models (Figure 2a and Supplementary Figure S1). The first day of involution (Inv1) presented clearly different gene expression pattern from those at mid-and late-lactation.
Most genes, particularly those following kinetics 1 and 2 (Figure 2c,d), were downregulated during involution, implying they were overexpressed at mid-and late-lactation, highlighting their central role for milk production and regulation [32]. Genes following kinetic 1 were related to the first steps of milk lipid and lactose synthesis [19], such as glucose and fatty acid uptake (Glut1, Lpl), fatty acid de novo synthesis (Fasn, Acaca) and galactose synthesis (Ugp2), to mammary energetic metabolism (Cs, Mtco1) and to milk synthesis regulatory factors (Stat5, Akt1, Pparg). These genes reached their peak expression between mid-and late-lactation and were probably under primary regulation of Prlr whose expression was highly correlated with all (r > 0.8). Genes following kinetic 2 were associated with the final steps of milk constituent secretion [19], such as triglyceride synthesis (Scd1, Dgat1), lipid droplet formation (Plin2), lactose synthesis (B4galt1, Lalba), main milk protein synthesis (Csn2, Wap) and water inflow (Aqp1, Aqp3). They seemed to be mostly controlled by Insr (r ≥ 0.75 with all) but also Esr1, mainly for milk protein genes (r ≥ 0.85). To the contrary, genes following kinetic 3 (Figure 2e) were upregulated at Inv1, suggesting that their expression was repressed during lactation. These genes are related to metabolic pathways repressed during lactation such as mitochondrial β-oxidation (Acads, Cpt1a) [32,33,42] or have a secondary role on regulation (Igf1r, Lxra) [33,43,44] or on lipid synthesis (Fabp4, Dgat2) [33,45] compared to other regulatory factors (Prlr, Srebf1) or other family members (Fabp3, Dgat1). Moreover, in the present study, the drop in the expression of de novo fatty acid synthesis genes at involution is correlated (Supplementary Figure S2) with the concomitant increase in the expression of β-oxidation genes likely related to energy requirements for mammary tissue remodeling. Finally, Oxtr was overexpressed in involuting mammary gland probably because myoepithelial cells were preserved longer than lactocytes and contracted to facilitate removal of degenerating epithelial cells [46].
Interestingly, genes related to milk constituent secretion followed the kinetic 2 pattern, and as such, were overexpressed at L18 compared to L12. These genes are probably the more strongly correlated with actual milk secretion, suggesting that, among the three periods studied, milk production was likely at its maximum at L18. This is consistent with (i) the time course of plasma insulin, leptin and IGF-1, which all declined to their minimal concentration at L18. Indeed, the plasma concentration of these hormones is known to decline during lactation partly due to increased hormone extraction by the mammary gland [47][48][49] and they are negatively correlated with milk production [50]. (ii) The highest milk consumption by rat pups reported to occur at L15-L16 as they only start to consume solid food significantly at L17 and drinking water at L19 [51].

Fenugreek Supplementation Stimulates Milk Macronutrient Synthesis
Dietary fenugreek supplementation increased the overall expression of genes involved in fatty acid uptake (Lpl), de novo synthesis (Fasn, Acaca, Acacb), transport (Fabp3) and, to a lesser extent, lipid secretion (Plin2). Fenugreek also stimulated the expression of genes involved in glucose uptake (Glut1), galactose (Pgm1, Ugp2) and lactose (B4galt1, Lalba) synthesis. Finally, fenugreek increased the expression of amino acid transporters (Slc7a5) and main milk proteins (Csn2, Lalba). These results suggest an increase in the synthesis of the 3 macronutrients of milk namely lipid, lactose and protein in the mammary gland [19] and are consistent with the stimulating effect of fenugreek on milk flow by 16% reported in our earlier study [15] and 15.3% in this study. More specifically, genes involved at every step of protein and lactose synthesis were significantly overexpressed, whereas only genes involved in the fatty acid synthesis were dramatically overexpressed for lipid synthesis. The expression of genes involved in triglyceride synthesis and lipid secretion was not or only slightly affected. This relative lack of conversion of fatty acids to triglycerides may account for the decrease in milk lipid concentration found in some fenugreek supplementation studies [52][53][54]. Otherwise, the stimulation by fenugreek of genes involved in mitochondrial energy metabolism (Cs, Atp5f1a, Mtco1) and particularly those involved in β-oxidation (Cpt1a, Acads) suggests that increased ATP synthesis took place in the mammary gland to meet the increased energy requirements associated with increased milk production [22]. Nevertheless, contrary to Liu et al. [55], we were unable to observe a significant increase in water inflow in response to galactagogue even though mean Aqp1 and Aqp3 expressions were increased at L12 and L21.

Fenugreek Supplementation Advances Peak Lactation and Maintains Milk Production
The lack of separation according to diet on the PLS-DA plot, and the lack of any difference in the expression of genes related to milk constituent synthesis and fuel metabolism in FEN and CTL groups at L18, suggests that the milk synthesis peak observed at late-lactation mainly in control dams was not impacted by fenugreek supplementation. However, the time course of some genes from L12 to L18 differed between diets: either the gene expression increased in the CTL diet and did not change in the FEN diet (Lpl, B4galt1, Slc7a5, Csn2 and Aqp1), or the expression did not change in the CTL diet and decreased in the FEN diet (Acaca, Fabp3, Ugp2, Cpt1a, Acads). In both cases, gene expression was higher earlier in the FEN group. This is consistent with the clear separation between both diets observed in the PLS-DA model built on gene expression at L12 (Figure 4d) and the fact that greater overexpression was observed at L12 for most of the genes in the FEN group. Finally, the PLS-DA model built on gene expression at Inv1 also properly separated both diets and the expression of the genes Lalba, Csn2, Plin2, Slc7a5, Aqp1 and Esr1 was maintained to a slightly greater level in the FEN diet. Together these results suggest that fenugreek supplementation advanced the peak of milk production, observed near L18 for control, around L12 and maintain it to a level similar to control at L18. Thus, fenugreek supplementation may enable the dam to maintain a milk production slightly higher than control during involution.

Fenugreek Stimulates Milk Production Mainly by Modulating Insulin and Oxytocin Secretion
The main effect of fenugreek on the expression of genes related to milk synthesis was observed at L12. During the same period, fenugreek supplementation increased dams' plasma insulin by 77.1% and the expression of Insr by 1.79-fold, suggesting an effect of fenugreek on insulin secretion and its action in mammary gland. Receptors to GH and IGF-1 (GHR and IGF1R) were also remarkably overexpressed in the mammary gland at L12 (Table 5). These 3 hormones are closely linked by the insulin/GH/IGF-1 axis [56]. Briefly, GH is secreted by pituitary gland and stimulates IGF-1 secretion in the liver [56,57], which shares strong homology with insulin. Both hormones can bind to IGF1R and insulin receptors (INSR) and to hybrid formed by both receptors with more or less affinity. These receptors act through similar pathway mainly by activating the AKT/mTOR pathway [56]. In addition, insulin stimulates expression of GH receptor in peripheral tissues [58] but also modulates the expression of IGF1R and INSR [59,60].
We therefore hypothesize that fenugreek acts mainly by increasing insulin secretion by maternal pancreas β-cells. This is consistent with the hypoglycemic effect of fenugreek and its capacity to stimulate insulin secretion through the potential action of trigonelline and 4-hydroxyisoleucin [23,24]. Together with a potential direct action of fenugreek compounds [24], the rise in plasma insulin under lactation-specific conditions, i.e., hypo-insulinemia and insulin-resistance in peripheral tissues except mammary gland, could explain the Igf1r, Insr and Ghr overexpression in the mammary gland [58,60,61]. The increase of IGF-1 and insulin binding to their receptor and to their hybrid receptor leads to an activation of the AKT/mTOR pathway, as suggested by overexpression of Mtor at L12 (by 1.84-fold) and of Akt1 to a lesser extent (Table 4) and results in activation of lactose and protein synthesis [19]. GHR could also act secondarily on protein and lactose synthesis through the JAK2/STAT5 pathway [62]. Concerning lipid metabolism, Lxra and Pparg are the most overexpressed regulatory factors with fenugreek supplementation (Table 4). Both are stimulated by insulin [63][64][65] and likely IGF-1, and Lxra expression is also stimulated by GH [64]. Lxra has been shown to stimulate de novo fatty acid in the mammary gland [66] independently of SREBF1 but not fatty acid desaturation and triglyceride synthesis [44], which is consistent with fenugreek effect found only for fatty acid uptake and de novo synthesis. Lxra was also related to activation of lipolysis in adipocyte [63,64] and may be responsible for the maintenance of β-oxidation-related genes at L12, as suggested by the strong correlation with Acads and Cpt1a (r > 0.7, Figure 2e). Pparg had a secondary milk lipid regulatory role in rodent compared to ruminant [19] but promote de novo fatty acid synthesis, triglyceride synthesis and lipid secretion [19,45,65]. Finally, Igf1 is known to be an important activator of fatty acid β-oxidation [56] that could counteract insulin inhibitory effect [20,42] on this pathway and notably at involution. It may act independently on the maintenance of Acads and Cpt1a expression, as suggested by the very strong correlation between Igf1r, Acads and Cpt1a expressions (r ≥ 0.84, Figure 2e). Suggested mechanisms associated to the galactologue effect of fenugreek are summarized in Figure 5.
Another mechanism by which fenugreek could increase milk flow may be its capacity to favor milk ejection by stimulating oxytocin secretion [22]. Indeed, fenugreek is a well-known oxytocic agent used to promote uterine contractions at delivery [23,25] and oxytocin expression was overexpressed in the pituitary gland during lactation in the fenugreek-supplemented diet, particularly at L18 (1.88-fold). Even if Oxtr expression was not significantly increased in mammary gland, fenugreek might maintain greater milk flow at L18 by promoting milk ejection ( Figure 5).
Nevertheless, we could not observe increase in mother prolactin, as suggested in previous studies [26,53], because of strong variability in plasmatic prolactin concentration. These variations were probably due to the too-short dam/pup separation time before sacrifice (60 min) [67], which made prolactin concentration more likely representative to the time since last suckling.

Fenugreek Stimulates Estrogenic Activity at the End of Lactation and Increased Energy Availability for Milk Production
At L18 and Inv1, fenugreek-supplemented dams tended to have higher plasma estrogens, and in the pituitary, had higher expression of C3, which is known to be a strong estrogenic biomarker in the uterus [68]. Moreover, Esr1 tended to be upregulated in the mammary gland at L21 (Figure 3g) and was strongly correlated with genes overexpressed at this period (r ≥ 0.74 with Csn2, Lalba, Aqp1, Plin2 and Akt1, Figure 2g and Supplementary Figure S1). As fenugreek is known to contain several estrogenic compounds [23,27], those might contribute to maintaining the expression of milk secretory genes at Inv1 partly through the mediation of Akt1 [69].
Finally, fenugreek supplementation increased global food consumption of dams by 10%, in line with earlier studies [52,70,71], principally at mid-lactation (Supplementary Figure S4) that could partly account for fenugreek lactogenic effect. Indeed, lactation is a period of high energy requirement characterized by an important hyperphagia, in particular in small animals where food intake can be increased by 450% [72]. Rats have little lipid storage during lactation and energy requirement for milk production almost exclusively comes from increased food consumption [22,73]. Thus, by enhancing food consumption, fenugreek might enhance energy availability for milk production. This increase in food intake observed after fenugreek supplementation might be mediated (i) by a family of its bioactive compounds, steroid saponins, that stimulate eating motivation [74,75] and could act directly on the hypothalamus or through stimulation of prolactin, or (ii) by a greater peak of prolactin during suckling, as suggested by overexpression of receptors of prolactin secretion activators in the pituitary gland (thyrotropin releasing hormone and vasoactive intestinal peptide) [76]. Indeed, lactation hyperphagia is under the primary control of prolactin which decreases leptin secretion in adipocytes and, in the hypothalamus, stimulates orexigenic neuropeptide Y secretion and decreases its sensitivity to leptin [22]. Moreover, milk production is also known as an important determinant of food intake in rat [72,77]. Indeed, greater milk yield leads to higher dam's energy expenditure and likely, enhanced suckling stimulus, as mammary gland was reported to be emptied at each feeding for a rodent litter-size model of 12 pups [78], resulting in greater food intake. This is in accordance with the peak dietary intake observed when fenugreek dams started their peak milk production. receptor and to their hybrid receptor leads to an activation of the AKT/mTOR pathway, as suggested by overexpression of Mtor at L12 (by 1.84-fold) and of Akt1 to a lesser extent (Table 4) and results in activation of lactose and protein synthesis [19]. GHR could also act secondarily on protein and lactose synthesis through the JAK2/STAT5 pathway [62]. Concerning lipid metabolism, Lxra and Pparg are the most overexpressed regulatory factors with fenugreek supplementation (Table 4). Both are stimulated by insulin [63][64][65] and likely IGF-1, and Lxra expression is also stimulated by GH [64]. Lxra has been shown to stimulate de novo fatty acid in the mammary gland [66] independently of SREBF1 but not fatty acid desaturation and triglyceride synthesis [44], which is consistent with fenugreek effect found only for fatty acid uptake and de novo synthesis. Lxra was also related to activation of lipolysis in adipocyte [63,64] and may be responsible for the maintenance of β-oxidation-related genes at L12, as suggested by the strong correlation with Acads and Cpt1a (r > 0.7, Figure 2e). Pparg had a secondary milk lipid regulatory role in rodent compared to ruminant [19] but promote de novo fatty acid synthesis, triglyceride synthesis and lipid secretion [19,45,65]. Finally, Igf1 is known to be an important activator of fatty acid β-oxidation [56] that could counteract insulin inhibitory effect [20,42] on this pathway and notably at involution. It may act independently on the maintenance of Acads and Cpt1a expression, as suggested by the very strong correlation between Igf1r, Acads and Cpt1a expressions (r ≥ 0.84, Figure 2e). Suggested mechanisms associated to the galactologue effect of fenugreek are summarized in Figure 5. Otherwise, we observed higher leptin and insulin level in plasma of fenugreek-supplemented dams likely caused by fenugreek stimulation of insulin secretion [79]. Yet, these hormones are known to exhibit important anorexigenic effect under normal physiological conditions [80]. However, during lactation, because of central resistance to their signals [22], insulin and leptin had only secondary roles in regulating lactation-associated hyperphagia compared to suckling stimulus and prolactin action [77,81]. Moreover, under normal conditions, orexigenic effect of fenugreek was observed despite an increase of plasma insulin [74]. Thus, effect of fenugreek on food intake was probably not related to the increase of both insulin and leptin plasma concentrations.
Whatever the cause, overall food increase in fenugreek-supplemented dams was of 10% and only of 5.3% (p = 0.041) related to dams' body weight, so it cannot alone account for the 15% increase in milk production observed in our study.

Conclusions
The current study confirms the ability of dietary fenugreek supplementation to promote milk production, and sheds light into its action at the molecular level. Indeed, milk flow was increased by 15.3% and most of the genes related to milk lipid, lactose and protein synthesis as well as energy metabolism were overexpressed. Moreover, our findings suggest that fenugreek may act by extending the duration of peak lactation until the beginning of involution rather than intensifying it. According to our data, the main mechanism of action of fenugreek on milk production is likely through activation of the insulin/GH/IGF-1 axis, which leads to a greater mammary sensitivity to these hormones. Consequently, milk protein and lactose synthesis could be mainly activated by the AKT/mTOR pathway, while milk lipid synthesis and energetic metabolism could be stimulated by LXRα, PPARγ and probably IGF1R itself. Moreover, fenugreek would stimulate milk ejection by an increase of pituitary oxytocin secretion. However, data from gene expression must be confirmed by proteomic and histologic analysis.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/10/1208/s1, Figure S1: PLS-DA score plot on non-normalized mammary gene expression separating groups according to the diet and the lactation period. Figure S2: Pearson correlation table between all mammary gland genes' expressions on row data. Figure S3: Pearson correlation table between all mammary gland genes' expressions on data normalized by the period of lactation. Figure S4: Food intake of dams supplemented or not with 1 g/kg/day fenugreek during lactation. Table S1: Validation of RNA extraction and primer design.