L-lactate induces specific genome wide alterations of gene expression in cultured bovine granulosa cells

Background Previously, we could show that L-lactate affects cultured bovine granulosa cells (GC) in a specific manner driving the cells into an early pre-ovulatory phenotype. Here we studied genome wide effects in L-lactate-treated GC to further elucidate the underlying mechanisms that are responsible for the L-lactate induced transformation. Cultured estrogen producing GC treated either with L-lactate or vehicle control were subjected to mRNA microarray analysis. Results The analysis revealed 487 differentially expressed clusters, representing 461 annotated genes. Of these, 333 (= 318 genes) were identified as up- and 154 (= 143 genes) as down-regulated. As the top up-regulated genes we detected TXNIP, H19 and AHSG as well as our previously established marker transcripts RGS2 and PTX3. The top down-regulated genes included VNN1, SLC27A2 and GFRA1, but also MYC and the GC marker transcript CYP19A1. Pathway analysis with differentially expressed genes indicated “cAMP-mediated signaling” and “Axon guidance signaling” among the most affected pathways. Furthermore, estradiol, progesterone and Vegf were identified as potential upstream regulators. An effector network analysis by IPA provided first hints that processes of “angiogenesis” and “vascularization”, but also “cell movement” appeared to be activated, whereas “organismal death” was predicted to be inhibited. Conclusions Our data clearly show that L-lactate alters gene expression in cultured bovine GC in a broad, but obviously specific manner. Pathway analysis revealed that the mode of L-lactate action in GC initiates angiogenic processes, but also migratory events like cell movement and axonal guidance signaling, thus supporting the transformation of GC into an early luteal phenotype. Electronic supplementary material The online version of this article (10.1186/s12864-019-5657-6) contains supplementary material, which is available to authorized users.


Background
Folliculogenesis is a finely tuned process of cellular differentiation. The most crucial step of differentiation is the folliculo-luteal transition that is initiated by the pre-ovulatory LH surge. Besides the release of a fertilizable oocyte in particular in the bovine this transition involves a deep transformation of somatic cells of the follicular wall into luteal cells. This is essential for regulating the ovarian cycle and supporting an ongoing pregnancy. In this transition phase the follicle is completely remodeled from the vesicle-like estradiol (E2)-producing dominant follicle into the compact progesterone (P4)-producing corpus luteum (CL). In the bovine, cells of the granulosa and theca layers migrate and largely intermingle with each other during CL formation [1]. This remodeling of the follicle is preceded and accompanied by a profound and meticulous regulation of gene expression in particular in the granulosa cell layer. Especially genes involved in steroidogenesis have been shown to be strongly regulated by the LH surge [2][3][4][5]. CYP19A1, encoding the key enzyme of estradiol synthesis (aromatase) is massively down-regulated together with the gonadotropin receptors FSHR and LHCGR. On the other hand several genes are highly up-regulated due to LH, namely RGS2 (regulator of G protein signaling 2), VNN2 (vanin 2) and PTX3 (pentraxin 3). VNN2 and PTX3 are involved in processes of inflammation. Moreover, PTX3 has been demonstrated to be essential for female fertility by organizing stable extracellular matrix architecture for an intact cumulus oophorus complex [6][7][8]. RGS2 interacts with the Gα subunit of G proteins by blocking Gα-mediated signaling [9] and was shown to modulate LH receptor signaling thus having an important role during the folliculo-luteal transition [10,11]. From this knowledge, typical markers for LH-dependent differentiation in the bovine could be established [3], whereby in particular the role of various growth-factors as those of the TGFbeta superfamiliy or EGF and their role during follicle differentiation were analyzed in detail [12][13][14][15]. In our previous study we could demonstrate that L-lactate, a molecule that is usually known to be connected to the energy metabolism, can act as a signaling molecule that specifically influences gene expression and thus remarkably influences GC differentiation in vitro [16]. Other studies could show that L-lactate is present at much higher levels within the follicular fluid than in serum [17,18]. Additionally, it was shown that the L-lactate levels in rats increase at the expected time of the LH-surge implicating a regulatory role that forces GC differentiation [19]. During the present study we analyzed the effects of L-lactate in a genome wide mRNA microarray approach with subsequent bioinformatic evaluation of the datasets to elucidate the underlying pathways and biological processes.

Tissue collection and cell culture
Ovaries were collected at a local abattoir irrespective of age, nutritional status or stage of ovarian cycle and transported in 1x PBS (with 100 IU penicillin, 0.1 mg/ml streptomycin and 0.5 μg/μl amphotericin). For each cell preparation granulosa cell pools were obtained by aspiration of small to medium sized follicles (< 6 mm) with a syringe and 18 G needle from 30 to 40 ovaries. By follicle aspiration almost exclusively cells from the granulosa without contaminating thecal cells were collected [4]. Living cells were counted using the trypan blue exclusion method and cryo-preserved in freezing media (fetal calf serum containing 10% DMSO; Roth, Karlsruhe, Germany). Additionally, a small portion of the freshly aspirated sample-pool was preserved as "pre-cultured sample" for qPCR analysis in liquid nitrogen. For seeding, the GC were thawed and immediately transferred into α-MEM and centrifuged at 500 x g for 3 min to ensure fast removing of the freezing media. Afterwards the cell pellet was diluted in α-MEM containing L-Glutamine (2 mM), sodium bicarbonate (10 mM), BSA (0.1% w/v), HEPES (20 mM), sodium selenite (4 ng/ml), transferrin (5 μg/ml), insulin (10 ng/ml), non-essential amino acids (1 mM), penicillin (100 IU/ml) and streptomycin (0.1 mg/ml), FSH (20 ng/ml; Sigma Aldrich, Steinheim, Germany), R 3 IGF-1 (50 ng/ml; Sigma Aldrich), and androstenedione (2 μM; Sigma Aldrich). Improved attachment of plated GC was achieved by coating the wells with Collagen R (0.02%; Serva, Heidelberg, Germany). Cells were cultured at a density of 1.0 × 10 5 living cells per well. Cells were additionally treated with sodium L-lactate (30 mM; Sigma Aldrich), sodium chloride as vehicle control (30 mM; Sigma Aldrich) or left untreated. If not stated otherwise all reagents were purchased from Merck Millipore (Berlin, Germany). GC were cultured for 8 days at 37°C and 5% CO 2 and two-third of the media with or without L-lactate or vehicle was exchanged every other day.

RNA preparation, cDNA synthesis and quantitative realtime PCR
Total RNA was isolated with the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Measuring the RNA concentration was done with a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Bonn, Germany). Afterwards, cDNA synthesis was performed using the SensiFAST cDNA Synthesis Kit (Bioline, Luckenwalde, Germany) from 150 ng RNA. Validation of the microarray data was accomplished with quantitative real-Time PCR (qPCR). Therefore, Sen-siFAST SYBR No-ROX (Bioline) was used with gene-specific primers as listed in Additional file 1: Table  S1. The amplification was done in duplicate from 0.2 and 0.4 μl cDNA in a total volume of 12 μl in a LightCycler 96 instrument (Roche, Mannheim, Germany). The ensuing cycle conditions were used: pre-incubation at 95°C for 5 min, 40 amplification cycles of denaturation at 95°C for 20 s, annealing at 60°C for 15 s, extension at 72°C for 15 s and single-point fluorescence acquisition for 10 s. At the end of each run, the melting point was analysed to check for amplification of the right products. Additionally, PCR products were checked by agarose gel electrophoresis (3%, stained with Roti-GelStain, Roth, Karlsruhe, Germany). As external standards for quantification cloned and sequenced products were used. Therefore, five different dilutions of the respective standards (5 × 10 − 12 -5 × 10 − 16 g DNA/reaction) were freshly prepared and co-amplified. To check for appropriate reference genes in this experimental setup the following commonly used genes were examined: B2M (Beta-2-microglobulin), GAPDH (Glyceraldehyde-3-phosphate dehydrogenase), RPLP0 (Ribosomal protein lateral stalk subunit P0) and TBP (TATA box binding protein). The two most stable reference genes were obtained by using the geNORM algorithm implemented in the NormqPCR package for R [20], revealing TBP and B2M as the most stable genes. For normalization the geometric mean of both was used.

Microarray profiling, bioinformatic evaluation and statistics
The RNA of samples from the three different culture conditions (untreated, L-lactate and vehicle control, 5 samples per group, n = 15) were subjected to mRNA microarray analysis. The quality of RNA was checked in a Bioanalyzer Instrument (Agilent Technologies, St. Clara, CA, USA) revealing a RIN factor ranging from 9.5 to 9.9 thus showing negligible degradation of individual samples. The Bovine Gene 1.0 ST Array (Affymetrix, St. Clara, CA, USA) was used for the analysis. Amplification, labelling and hybridization was performed with the "GeneChip Expression 3' Amplification One-Cycle Target Labeling and Control Reagents" (Affymetrix) according to the manufacturer's protocol. Hybridization was done overnight in the GeneChipR Hybridization Oven (Affymetrix) and visualized with the Affymetrix Gene-Chip Scanner 3000. Raw data were processed with the Expression Console (V1.4.1.46; Affymetrix) for normalization, background reduction and a gene-level summary using the RMA method (Robust Multichip Average). Additionally, a principal component analysis (PCA) was performed and plotted in R [21]. Results of the array have been submitted to the GEO database (GSE121408). Subsequent analysis was done with the Transcriptome Analysis Console 3.0 (TAC3.0, Affymetrix) to check for differentially expressed genes under different conditions. Analysis of Variance (ANOVA) was used to calculate the p-value and was additionally corrected for FDR (False Discovery Rate, Benjamini-Hochberg method) integrated in TAC3.0. Levels of significance for differentially expressed genes were set at fold change (FC) > |1.5|, ANOVA p < 0.05 and FDR < 0.05.
Statistic evaluation of qPCR values was performed with the SigmaPlot 11.0 Statistical Analysis System (Jandel Scientific, San Rafael, CA, USA). Threshold for statistical significance was set at p < 0.05.
Bioinformatic analysis was done with Ingenuity Pathway Analysis (IPA, Qiagen) using a modified list of the treatments L-lactate vs. vehicle control containing 2429 transcript clusters (FC > |1.2|, p < 0.05, FDR < 0.05). Out of the modified list 2193 genes could be mapped in IPA, settings for the pathway analysis were restricted to genes with FC < |1.5|, p < 0.05, FDR < 0.05. All genes having a FC between |1.2| and |1.5| were characterized as abundant in the dataset but did not influence the pathway analysis.

Results
Raw data from the microarray analysis were initially analyzed by principal component analysis (PCA) to reduce the multidimensionality of the dataset. Individual samples of the dataset were plotted and revealed the greatest variability between the different groups of culture conditions on the x-axis with a variation of 30.9% ( Fig. 1a and Additional file 1: Table S2). Individual samples of the same culture condition clustered tightly together. Although we detected a difference between the untreated GC and the cells treated with the NaCl vehicle control, the L-lactate treatment was most distant compared to both.
On the Bovine Gene 1.0 Array Chip 26,288 transcript clusters are represented. Comparing the different groups "untreated", "L-lactate" and "NaCl" diverse numbers of differentially expressed clusters can be observed (Fig. 1b). The higher number (788) of differentially regulated clusters (representing 735 annotated genes) of the L-lactate treatment compared to the untreated cells is in line with the PCA, demonstrating strongest effects between L-lactate treatment compared to both controls ("untreated" and "NaCl"). Least changes were observed between "untreated" and "NaCl" vehicle treated cells. Comparing L-lactate vs. NaCl vehicle control treatment revealed 487 affected clusters (representing by 461 annotated genes).

Validation with qPCR
To validate the mRNA microarray expression datasets qPCR analysis of selected transcripts was performed. Except for SLC16A1 and SLC16A7, which both are not significantly regulated by L-lactate, nearly all selected genes analysed, showed high correlations between the qPCR and microarray datasets ( Table 1). Principal Component Analysis of the qPCR dataset revealed the highest variation of 71.7% between the L-lactate treatment and both other culture conditions "NaCl" and "untreated" controls ( Fig. 2a and Additional file 1: Table S3), with these clustering closely together. Freshly isolated, non-cultured cells were clearly separate (Fig. 2a, black dot) from cultured samples, but clearly showed more closeness to the control than the lactate treated samples considering PC1. Differential transcript concentrations in L-lactate vs. NaCl vehicle treated controls could be validated by qPCR for all selected genes (Fig. 2b). Even the fold change for most analysed genes was similar, except for AHSG (8.45 vs. 2.86), HAS2 (7.31 vs. 16.12) and TXNIP (21.97 vs. 100.69).

Differentially expressed genes
Comparing the effects of lactate treatment with the NaCl vehicle control 487 transcript clusters (= 461 annotated genes) could be identified as differentially expressed. From these, 333 transcript clusters were assigned as up-and 154 clusters as down-regulated by lactate. The top 15 up-and down-regulated genes were listed in Tables 2 and 3. The highest up-regulation was found in case of TXNIP with a fold change (FC) of 21.97 followed by the non-protein coding gene H19 (FC 12.36). PTX3, one of our formerly established markers for pre-ovulatory differentiation, was also found under the top 15 of up-regulated genes. The most down-regulated gene was VNN1 (FC -2.82) in contrast to VNN2 which was up-regulated as expected according to our previous data (FC 1.58, Additional file 1: Table  S4). CYP19A1, another important GC marker was also among the top 15 down-regulated genes (FC − 2.29). Interestingly, the dataset revealed remarkably higher scores of up-regulation (FC > 21) than down-regulation (FC > − 2.8).
Pathway analysis by IPA 2193 out of 2429 (i.e. 90.3%) differentially expressed transcript clusters could be assigned to specific genes and mapped to specifically affected pathways and biological functions. Pathway analysis identified the "cAMP-mediated signaling pathway" as well as "Axonal Guidance Signaling" and "TGF-β Signaling" as significantly affected although no prediction regarding activation or inactivation could be made ( Fig. 3 and Additional file 1: Table S5). Further analysis identified TNF, beta-estradiol, progesterone and Vegf as major upstream regulators, which might be involved as activating factors and thus responsible for the observed changes of the expression profile (Table 4 and Additional file 1:  Table S6). Interestingly, the Regulator Effects analysis of IPA identified a putative activation of the functions "proliferation", "vascularization", "angiogenesis" or "cell movement" whereas the biological function "organismal death" was predicted to be inhibited (Fig. 4). This is in accordance with the observation that no significant regulation of pro-apoptotic factors like CASP4, CASP8 or TP53 (FC -1, 1.18 and − 1.3) could be observed. In this effector analysis AREG and EGR2 were identified as upstream regulators leading to the activation or inactivation of these functions. However, both were not among the top Upstream Regulator candidates (Table 4).

Discussion
Our data clearly show that any treatment with either lactate or the NaCl vehicle control significantly changed the global gene expression profiles of cultured GC when compared to untreated controls. However, it is also obvious that L-lactate induced strongest effects: 487 and 788, but only 266 genes were affected by NaCl treatment as compared to untreated cells. The effect of NaCl alone is clearly visible in the PCA of the microarray dataset that separates the vehicle control and untreated cells, revealing the sensibility of the cell culture model to the composition of media. In addition, these data also clearly underline the need to meticulously observe the respective culture conditions in order to ensure reproducibility between experiments also in other cell culture models. However, focusing on markers of GC differentiation specific effects of L-lactate are even more clearly visible without any unspecific changes in the gene expression profile. In any case these data clearly justify our experimental approach to use NaCl treated cells as valid controls to exclude transcripts from the analysis that were affected just due to increased NaCl concentrations. In contrast, cryo-preservation before culture does not alter gene expression. A previous study comparing GC cultured directly or after cryo-preservation revealed no significant differences [22]. However, it is also clear that culturing by itself remarkably changed the gene expression profile of GC (see Fig. 2), thus underlining that cell culture models like that used during the present study can only partly mimic the in vivo situation. This limits the extrapolation of the data obtained in vitro to in vivo conditions.    In earlier studies we could show that various parameters affect GC differentiation in vitro as cell density, hypoxia and the supplementation of L-lactate [16,23,24]. In this study we investigated the global change of gene expression by comparing cells treated with L-lactate or NaCl vehicle control. Former established marker genes for differentiation, e.g. CYP19A1 or RGS2 showed a specific down-or up-regulation due to L-lactate treatment as described before [16]. Interestingly, the expression of GC identity marker FOXL2 was not affected by L-lactate thus indicating that L-lactate treatment does not change their identity throughout the culture period as found during treatment with oleic acid [25]. TXNIP, coding for the thioredoxin interacting protein, was remarkably up-regulated (FC 21.97) in L-lactate treated GC indicating a role in the L-lactate induced differentiation process. In contrast, when GC were cultured at high cell density we observed the opposite effect and a tremendous down-regulation of TXNIP (FC -79.5) [26]. Possibly, the regulation of TXNIP might be a sensor for glucose usage and metabolism as it regulates glucose uptake with increased expression reducing glucose uptake [27,28]. The strongly reduced expression in high density GC culture model might therefore mirror the need of the cells for increased glucose uptake under these "glucose-deficient" conditions, whereas increased expression in the present L-lactate supplementation model could be a consequence of the ample supply with an alternative energy source thus reducing the need for glucose uptake.
Also H19 gene expression was observed to be significantly up-regulated in L-lactate treated GC (FC 12.4). H19 is an imprinted gene of which only the maternal allele is transcribed into a long non-coding RNA [29], known to counteract/regulate the transcription of the paternally imprinted gene IGF2, an early growth factor that influences the size of the offspring at birth [30]. H19 expression is mostly abundant in fetal organs, although a moderate expression of H19 was found in adult ovarian tissue [31]. In an earlier study it was shown that steroid hormones can induce the expression of H19, which is therefore highly expressed in hormone sensitive organs [32]. Additionally it was proposed that H19 expression is high when the organ or tissue is under extensive re-modulation on physiological and morphological levels. Hence, the massive up-regulation of H19 in lactate treated GC might reflect an initiation of tissue re-organization as it can be found during the folliculo-luteal transition phase.
As a top down-regulated gene in L-lactate treated cells VNN1, a GPI-anchored protein with pantetheinase activity, was identified (FC − 2.8). As a regulator of the tissue response to oxidative stress VNN1 modulates the glutathione store [33]. In VNN1 knockout mice reduced inflammation and apoptosis could be observed [33]. Within the follicle an increase in VNN1 expression levels was proposed as an indicator for follicle growth but might also reflect atretic follicles [34,35]. The down-regulation of VNN1 in our cell culture model therefore suggests that GC under conditions of increased L-lactate do not have any disposition to atresia. This is also in line with the Regulator Effector Network Analysis predicting inhibition of "organismal death" (Fig.4).
MYC was shown to be down-regulated in the L-lactate treated GC compared to the vehicle control (FC − 2.4). Myc acts as a ubiquitous transcription factor, which targets several genes thus enhancing their expression [36]. It was also stated that MYC expression declines during differentiation, when a finely tuned reprogramming takes place. Otherwise, the enhancement by MYC would lead to uncontrolled proliferation. In this context the down-regulation of MYC in GC indicate cellular differentiation processes taking place under increased L-lactate conditions.
Top-affected pathways identified by the pathway analysis were "cAMP-mediated Signaling", "Axonal Guidance Signaling" as well as "TGF-β Signaling". cAMP associated pathways contribute to multiple biological processes either under physiological or pathological conditions [37]. Mainly two different intracellular recipients were identified: proteinkinase A (PKA) and the exchange protein directly activated by cAMP (Epac) [37][38][39]. Earlier it was proposed that both PKA and Epac are involved in the process of luteinization activated by LH [40][41][42][43]. In our former study of density-driven differentiation in bovine GC we could also highlight the involvement of the "cAMP-mediated Signaling" pathway [26]. Thus, the results of the present study reflect a LH induced like differentiation of bovine GC upon L-lactate treatment.
The "Axonal Guidance Signaling" pathway likely affected by L-lactate is known to be involved in GC differentiation in connection to cumulus expansion [44,45]. In particular NTN1 (netrin-1) was identified during neuron formation. But emerging evidence also exists that netrin-1 is a critical component in vascular regulation [46,47], as well as in promoting angiogenesis [48,49], which was also postulated for Netrin-4 in placenta [50]. NTN1 was down-regulated with a fold change of − 2.1 in L-lactate treated GC compared to the control, which would suggest that vascular regulation or angiogenic processes are not induced. Interestingly, netrin-1 was found to be present in follicular fluid as well as in the theca and granulosa cell layer of swine antral follicles and was proposed to have anti-angiogenic functions [51]. However, whether netrin-1 acts as angiogenic factor or not is still a matter to debate [52,53]. Our data suggest that NTN1 is an anti-angiogenic factor regarding the down-regulation and the putative activation of angiogenesis. But the activation of angiogenesis is not only related to down-regulation of netrin-1, moreover, other more prominent factors are involved as well, e.g. AREG or CCND1. It seems that the final answer of NTN1 function regarding angiogenesis is not that straightforward in granulosa cells and needs further investigations.
On the other hand the Slit/Robo pathway is associated with the process of axonal guidance. Several subunits of SLIT and ROBO were present in the microarray dataset but no differential expression of these could be observed. Nonetheless, SLIT and ROBO expression could be detected in human luteinized GC or in the CL and are regulated by steroid hormones [54]. Previously, we discussed the involvement of NMDA receptors in mediating the L-lactate effects as it was shown for neurons [16,55]. In the microarray dataset we could identify the expression of several subunits of NMDARs in bovine GC, however without any differential regulation of these. The highest expression could be identified of GRIN2D (signal intensity of 4.6-4.8) and GRIN2C (signal intensity of 4.0-4.4) similar to another study of bovine GC in vivo and in vitro [56]. NMDARs are important receptors in axonal guidance and synapse formation [57]. However, if L-lactate signaling might be mediated via NMDA receptors in bovine granulosa cells still has to be elucidated.
The upstream regulator analysis implemented in IPA revealed beta-estradiol and progesterone having an activating influence on the L-lactate treated cell culture model. Classically, estradiol elicits a positive feedback on the hypothalamus that regulates the GnRH secretion. GnRH controls the release of the gonadotropins FSH and LH by diverging pulse frequencies. Additionally, a negative feedback mechanism on FSH secretion at the pituitary is also known [58,59]. Both actions of estradiol trigger the LH surge leading to ovulation. Interestingly, progesterone was also proposed in having an impact on the L-lactate-mediated changes observed in our cell model, although its concentrations did not change during culture. Progesterone, on the other hand, is a critical parameter for the establishment of an active corpus luteum indicating a transition of the GC's phenotype towards luteinization. However, whereas expression of the key gene of progesterone synthesis HSD3B1 is very high in fully luteinized GC (i.e. large luteal cells), shortly after the LH surge, during the folliculo-luteal transition phase, its expression is even slightly reduced as compared to that in GC isolated from large dominant follicles [3,4]. The observation that expression of HSD3B1 was almost unchanged after L-lactate supplementation (FC 1.04) indicates that the cells are not fully luteinized, but may have yet adapted only to an early post-LH, but pre-ovulatory phenotype. Vegf as upstream regulator as well as TGF-β signaling indicate the activation of angiogenic factors. It is commonly known that angiogenic processes contribute to ovulation and the later formation of the corpus luteum [60,61]. Also the IPA Effector Analysis revealed "angiogenesis" or "vascularization" as molecular functions to be activated involving the upstream regulator AREG. The function "organismal death" could be identified as inactivated, thus indicating that L-lactate treatment does not affect the viability of the cultured GC. In addition, the transcription of apoptosis markers like BAX or BCL2 was not induced in L-lactate treated cells thus supporting our assumption that the cells are not driven towards atresia. The function "cell movement" could be identified as activated, which is in line with the forthcoming breakup of follicular cell layers and necessary migratory processes during the formation of the corpus luteum.

Conclusions
Taken together our data provide novel insights into a possible regulatory role of increased concentrations of L-lactate on cells of the granulosa in large follicles during the folliculo-luteal transition. Our data suggest that the biological function of L-lactate in the granulosa cell layer of growing follicle is complex and by far exceeding its role as a product of hypoxic metabolism and energy source. It seems involved in different signaling pathways thus influencing the expression of many different genes. As a commonly known pathway of folliculogenesis our data suggest PKA signaling to be associated with the L-lactate effects. However, we also gathered first hints that NMDAR signaling usually found in neuron physiology might be involved in differentiation processes induced by L-lactate.

Additional file
Additional file 1: Table S1. List of gene-specific primers used in qPCR. Table S2. Principal component analysis of microarray dataset. Table S3. Principal component analysis of qPCR . Table S4. Differential expressed clusters of lactate vs. vehicle control -treated GC. Table S5. Affected pathways identified by IPA. Table S6