Lipidomics profiling of goose granulosa cell model of stearoyl-CoA desaturase function identifies a pattern of lipid droplets associated with follicle development

Despite their important functions and nearly ubiquitous presence in cells, an understanding of the biology of intracellular lipid droplets (LDs) in goose follicle development remains limited. An integrated study of lipidomic and transcriptomic analyses was performed in a cellular model of stearoyl-CoA desaturase (SCD) function, to determine the effects of intracellular LDs on follicle development in geese. Numerous internalized LDs, which were generally spherical in shape, were dispersed throughout the cytoplasm of granulosa cells (GCs), as determined using confocal microscopy analysis, with altered SCD expression affecting LD content. GC lipidomic profiling showed that the majority of the differentially abundant lipid classes were glycerophospholipids, including PA, PC, PE, PG, PI, and PS, and glycerolipids, including DG and TG, which enriched glycerophospholipid, sphingolipid, and glycerolipid metabolisms. Furthermore, transcriptomics identified differentially expressed genes (DEGs), some of which were assigned to lipid-related Gene Ontology slim terms. More DEGs were assigned in the SCD-knockdown group than in the SCD-overexpression group. Integration of the significant differentially expressed genes and lipids based on pathway enrichment analysis identified potentially targetable pathways related to glycerolipid/glycerophospholipid metabolism. This study demonstrated the importance of lipids in understanding follicle development, thus providing a potential foundation to decipher the underlying mechanisms of lipid-mediated follicle development.

organelle dysfunction, ultimately leading to cell death and tissue damage [5,6]. Evidently, deregulation of this network contributes to the onset of pathology as lipids are tightly regulated, both spatially and temporally, in various parts of the organism [7,8]. Many lipid molecular species have been identified as potential biomarkers involved in several diseases, including, epidemic diseases, cancers, inflammation, and genetic diseases [9].
Intracellular lipid droplets (LDs) occur through the deposition of lipid complex molecular pathways. They are found in a wide variety of cell types, and are viewed as highly dynamic organelles with unique spherical structures, of which primary neutral lipids are mainly sterol esters and triacylglycerols (TGs) [10,11]. LDs affect cell and organ health and contribute to the dynamics of biological membranes, energy storage and expenditure, and conveying messages at all cellular levels, including in the nucleus [12,13]. Despite their important functions and nearly ubiquitous presence in cells, many aspects of LD biology are unknown. We have previously confirmed that de novo lipogenesis occurs in goose GCs [14] and that LD accumulation capacity depends on the stage of granulosa cells (GCs) during goose follicle development [15]. It has been reported that lipid metabolism in bovine, sheep, and human GCs is pivotal for oocyte maturation, fertilization, and early embryonic development [16][17][18][19]. Lipid content in follicles differs among species, as reported in studies performed in porcine [20], bovine [21], and canine oocytes [22], wherein LD number, size, and distribution were analyzed. However, an understanding of the biology of LDs in GCs remains limited. Recently, the role of stearoyl-CoA desaturase (SCD) as a regulator of fatty acid desaturation in goose GCs was identified, whereby metabolite changes in lipid-related metabolism were found to be closely related to follicular development [23]. Although this observation indicated that such links might exist, it is unclear whether LDs play a role in goose follicle development and which lipid species or classes might be affected. It is therefore hypothesized that lipidomic analysis of the goose GCs at the level of individual lipid molecular species might help to elucidate such a mechanism.
Lipidomics comprises the characterization of an organism's differences in individual lipid molecular species and their biological roles with respect to the expression of proteins involved in lipid metabolism and function [24,25]. Thus, to obtain a detailed view of lipid metabolism, it is crucial to analyze the full lipid profiles of the individual lipids using a lipidomics approach. This is the first time that a combination of lipidomics and transcriptomics analysis has been applied to the GC model of SCD function to identify alterations in the regulation of specific gene, lipid classes, and lipid-mediated signaling processes that are involved in LDs associated with follicle development in geese. The detailed characterization of lipid classes and species provided in this study might be useful to understand the mechanism of intramyocellular LDs, and an analysis based on integrating two-omics datasets will broaden our understanding of lipid alterations with respect to lipid-gene networks in GCs associated with follicle development.

LD distribution in goose GCs
Previous studies have demonstrated that de novo lipogenesis occurs in goose ovarian GCs, which have lipid accumulation capacity [15]. To investigate the localization of LDs in goose GCs, confocal microscopy was used to analyze LD location in goose GCs stained with BODIPY (500/510, green). GCs contained numerous internalized LDs that were generally spherical in shape and dispersed throughout the cytoplasm (Fig. 1). The dispersed LDs clumped to form clusters. The mobility of the dispersed LDs was from the perinuclear regions to the extracellular domain.

Altered SCD expression affects LD content
To further investigate the effect of SCD expression on LDs in GCs, the LDs were stained with Oil Red O after SCD-overexpression and knockdown (Fig. 2a, b). The accumulation of LD content was quantified using appropriate equations (see "Materials and methods" section for details). Compared with that in the control group, the overexpression of SCD significantly stimulated the accumulation of LD content. However, SCD knockdown significantly inhibited the accumulation of LDs (Fig. 2c).

Analysis of lipid classes
In total, 662 species of lipids were identified in the GC samples using LC-MS/MS analysis (Additional file 1: Fig. S1 Table S1). Venn diagram analysis revealed 26 overlapping lipids in the SCD-overexpression group, 35 overlapping lipids in the SCD-knockdown group, and nine overlapping lipids in the SCD-overexpression and SCD-knockdown groups (Additional file 2: Fig. S2). As shown in Table 1 and Fig. 4, 19 lipids followed the same trend, whereas seven lipids showed the opposite pattern, as compared to that of the two overexpressing groups (LN vs. LS comparison and LG vs. LS comparison). Thirty-three lipids followed the same trend, whereas two lipids showed the opposite pattern, as compared to that of the two knockdown groups (LC vs. LT comparison and LC vs. LF comparison). Seven lipids followed the same trend, whereas two lipids showed the opposite pattern between the overexpressed and knockdown groups.

Cross comparisons between and within each group
With the aid of the KEGG database, metabolic pathways are summarized in Additional file 7:  inositol phosphate metabolisms altered significantly; for LC vs. LT, glycerophospholipid, sphingolipid, glycerolipid, and inositol phosphate metabolisms altered significantly; and for LC vs. LF, glycerophospholipid, sphingolipid, glycerolipid, and inositol phosphate metabolisms altered significantly.

RNA-seq investigation of overexpression and knockdown of SCD expression in GCs
To obtain a comprehensive view of the role of the SCD gene in GCs, the six groups of cells (  The cells were after SCD knockdown. c Relative LD content. The results are shown as mean ± SD from six experiments, the data were analyzed by ANOVA and Tukey's test. The lowercase letter indicates significant differences between each group (P < 0.05)

Identification and functional enrichment of DEGs
In total, 679 DEGs were identified in the GCs after SCD overexpression and knockdown (Additional file 3:  (Fig. 6a, b).
To characterize differences in lipid-related genes between the transcriptomes of the SCD-overexpression and SCD-knockdown groups, the GO slim terms of the complex Gene Ontology annotations were used. The results showed that DEGs were assigned to 20 lipidrelated GO terms in the NC vs. OS comparison, 21 lipid-related GO terms in the OG vs. OS comparison, 63 lipid-related GO terms in the SC vs. ST comparison, and 63 lipid-related GO terms in the SC vs. SF comparison ( Fig. 7; Additional file 9: Table S4). The top 20 pathways of KEGG functional analysis showed that SCD-overexpression had a significant effect on one carbon pool by folate, glycosphingolipid biosynthesis, ferroptosis, and influenza A. Knockdown of SCD had a significant effect on sphingolipid metabolism, proteasome, mTOR signaling pathway, hedgehog signaling pathway, glycosphingolipid biosynthesis (metabolism), ether lipid metabolism, amino sugar and nucleotide sugar metabolism, and protein processing in the endoplasmic reticulum (Additional file 4: Fig. S4a, b).

Level of coordination between transcript and lipidomic data
Monitored responses in parallel on the transcriptome and the lipid level were carried out, thus allowing one to compare the level of coordination between both   Fig. 4 Heat map from the hierarchical clustering of differential overlapping lipid in the (a) SCD-overexpression group. b SCD-knockdown group. c SCD-overexpression and SCD-knockdown groups. The scaled expression by row (lipids) is shown as a heat map and is reordered by a hierarchical clustering analysis (Pearson's distance and Ward's method) on both rows and columns. The color scale indicates the relative amounts of metabolites: red, higher levels; green, lower levels phosphatidylinositol signaling system, alpha-linolenic acid metabolism, arachidonic acid metabolism, and linoleic acid metabolism. Additionally, a targeted approach using prior biological knowledge in conjunction with canonical correlation analysis (CCA) revealed that lipids from glycerolipid/glycerophospholipid metabolism were subjected to a CCA together with transcript data of all those pathways as derived from SCD knockdown in GCs (Fig. 9, Additional file 10: Table S5). We were also interested in other lipid-related pathways subjected to a CCA together with transcript data (Additional file 4: Fig. S4). However, no similar CCA results were observed with SCD overexpression in GCs.

Discussion
In the present study, we used a GC model of SCD function and performed a detailed study on cellular lipid metabolism using lipidomic profiling to evaluate the lipid molecular changes in a detailed manner. The method is based on ultra-high performance liquid chromatography combined with quadrupole-time-of-flight mass spectrometry and gas chromatography mass spectrometry, which allow for the analysis of hundreds of molecular lipids from a single sample run [26]. More than 660 lipids were identified, including phospholipids, ceramides, diacylglycerol and triacylglycerols, cholesterol esters, and fatty acid compositions (Fig. 3). Among these lipids, glycerophospholipids and glycerolipids were the predominant species. Additionally, altered SCD expression had the greatest impact on glycerophospholipid metabolism, this lipid-related pathway contains the most key lipids, with more than six representative metabolites as follows: There is great diversity in the structures and functions of the vast array of lipid species, with intracellular levels of lipids tightly regulated by a network of metabolic pathways to sustain normal cellular functions [27]. Glycerophospholipids make up phospholipids, which play both structural and metabolic roles in cells. Therefore, the level and species of phospholipids might reflect the biochemical and physiological conditions of cells [28]. Phospholipids are important for determining LD formation and are attributed to the formation of large LD fusions [29]. Glycerolipid metabolism is altered by fatty acids and, induces LD formation [30]. Blocking PC synthesis causes elevated SREBP-1-dependent gene transcription and LD accumulation in mouse liver and human cells, and the ratio of PE to PC in the LD surface tends to affect  the LD number and size [31,32]. The present study provided the lipid species composition and lipid metabolism of the GC model of SCD function. Given the important biological activities of many lipids, such information could provide potential insights into the mechanisms of internalized lipid-mediated follicle development.
We further investigated whether altered SCD expression of GCs affects gene expression involved in lipid metabolism. We observed that altered SCD expression resulted in many DEGs assigned to lipid-related GO terms, including lipid binding, lipid transport, lipid localization, and LDs. Compared with those in the  Table S4 SCD-overexpression group, more DEGs were identified and assigned to lipid-related GO terms in the SCDknockdown group (Fig. 7; Additional file 9: Table S4). These results indicated that SCD regulates the lipid metabolic processes during goose follicle development. Particularly, SCD knockdown might mobilize more genes associated with lipid metabolism. Combined analysis of lipid and gene expression profiles showed that more complicated lipid-related metabolism pathways were significantly clustered, especially in the SCD-knockdown group.
Numerous internalized LDs were dispersed throughout the cytoplasm of GCs (Figs. 1, 2). However, in steroidogenic tissues such as the ovary, the ovarian follicles of almost all mammals contain LDs, which are prominently and generally non-pathogenic, except in cases of genetic disorders of steroidogenesis [33]. We speculate that there is a relationship between LDs and energy supply or the  (Fig. 4c). Previous results demonstrated intracellular metabolite changes in which cholesterol was mostly upregulated in the SCD-overexpression group, whereas androsterone was upregulated in the SCD-knockdown group [23].
LDs make dynamic contact with nearly every other organelle. Some contact activity might exist among the ER, endosomes, mitochondria, and peroxisomes in the cell. This contact is believed to mediate the transfer of lipids (including phospholipids, sterols, and fatty acids) between the different cellular compartments [34,35]. Pharmacological inhibition of the enzymatic activity of fatty acid synthase significantly reduces the production of progesterone in bovine GCs, suggesting that ovarian steroidogenesis relies on lipid metabolism in GCs [16,36]. Additionally, LD accumulation capacity of goose GCs depends on the different stages of follicle development, with the stage of the highest accumulation capacity being F1 GCs, followed by F4-F2 GCs and pre-hierarchical GCs [15]. As ovarian follicles grow with changes in both the size and number of LDs in bovines [37], the total intracellular LD content increases gradually with follicular growth. This exponential increase was demonstrated in large antral follicles in rats and goats [38,39], with an increase in LD numbers and size in GCs within one day after ovulation in rabbits [40]. These results confirm that LDs are considered a source of energy for oocyte maturation [41,42].
SCD catalyzes the rate-limiting step in the production of monounsaturated fatty acids, which are major substrates for the synthesis of triglycerides, wax esters, cholesterol esters, and phospholipids [56]. Recent research has revealed that alterations in SCD expression change the fatty acid profile of these lipids and have a significant role in lipid metabolism [43]. The integrated optical density of the LDs is related to the amount of triglycerides in the LDs. Our results showed that altered SCD expression affected LDs in GCs. The overexpression of SCD stimulated the accumulation of LD content, whereas SCD knockdown inhibited the accumulation of LD content (Fig. 2c). Previous research has demonstrated that SCD is an important control point in lipid homeostasis in goose GCs [23]. Similar results suggest that increased hepatic SCD activity accelerates the last stage of TG synthesis and stimulates lipid accumulation [44]. These results provide direct evidence that SCD expression plays a major role in the triacylglycerol and phospholipid secretion processes. Much evidence indicates that the direct antisteatotic effect of SCD deficiency stems from increased fatty acid oxidation and reduced lipid synthesis [45] and that the inhibition of SCD decreases the proliferation and apoptosis of cells [46]. Additionally, a link between SCD activity and LD size was previously observed in the  Table S5 nematode Caenorhabditis elegans [47]; and experiments on cell lines cultured from patients with Berardinelli-Seip congenital lipodystrophy [48] demonstrated that SCD activity is necessary for large-sized LDs. Recent studies have indicated that the size of LDs might influence the balance of lipogenesis and lipolysis in adipocytes and hepatic cells [49], with decreased SCD expression resulting in reduced accumulation of oleic and palmitoleic acids and consequently fat cells with small LDs [50]. However, the functional implications of the diverse sizes of LDs in goose GCs are poorly understood.
In summary, these results provide new evidence that altered SCD expression changes GC lipid profiles, genes, and regulated relevant metabolic pathways, which might implicate SCD as one of the major genes in the regulation of lipid secretion by goose GCs.

Animals and granulosa layer isolation
Geese (from a maternal line of Tianfu goose) were raised under natural temperature and light conditions at the experimental station of waterfowl breeding at Sichuan Agricultural University. Six geese with regular laying sequences were randomly selected as experimental samples and sacrificed by cervical bleeding under anesthesia, 2 h after oviposition. A pool of ovarian follicles was immediately collected from the goose abdominal cavities and divided into prehierarchal follicles and hierarchal preovulatory follicles according to previously reported nomenclature [51]. For granulosa layer collection, the outer connective tissue was removed from the preovulatory follicles, and the preovulatory follicles were dissected to allow the yolk and adhered granulosa layer to flow out [52].

Cell culture and transfection
The granulosa layer was dispersed by incubation in 0.1% type II collagenase (Sigma, USA) in Dulbecco's modified Eagle's medium (DMEM, HyClone, USA) for 10 min in a 37 °C water bath. After incubation, cells were dispersed with a pipette and pelleted by centrifugation at 1000×g for 10 min (20 °C). The supernatant was discarded and the cells were resuspended in 3 mL of fresh basic medium without collagenase and centrifuged. The washing procedure was repeated twice. The GCs were dispersed in DMEM supplemented with 1% antibiotic/antimycotic solution (Solarbio, China) and 3% fetal bovine serum (Gibco, USA). Transient transfections based on the GC cellular model of SCD function (SCD-specific overexpression and knockdown) were performed using Lipofectamine ® 3000 and Lipofectamine ™ RNAiMAX Transfection Reagent (Invitrogen), according to manufacturer's protocol. SCD siRNA-210 and siRNA-405 were used to achieve SCD mRNA knockdown, with scrambled siRNA as a negative control. The extent of SCD mRNA knockdown was measured as a percentage compared to that the scrambled siRNA. SCD-specific overexpression was used to achieve SCD mRNA overexpression, termed pEGFP-N1/SCD, and the pEGFP-N1/empty vector (GFP vector) served as the negative control, and another control with no transfection components was also included. The extent of SCD mRNA overexpression was measured as a percentage compared to that with the GFP vector and no transfection components conditions. Transfection efficiency was determined as previously described [23].

Oil red O staining
To determine LD accumulation in GCs after transfection, the cells were stained with Oil Red O solution (Sigma Chemical Co., St. Louis, MO, USA). After 48 h of transfection, the cultured GCs were washed with phosphatebuffered saline (PBS) and fixed with 10% formaldehyde at room temperature. Then, the cells were washed twice with PBS, and were subsequently stained with 0.5 μg/mL Oil Red O solution in the dark for 1 h at room temperature. Then, cells were washed with PBS and photographed using an optical microscope system (Nikon ECLIPSE 90i) at 200 × magnification. Finally, the LDs were dissolved in isopropanol and the optical density (OD) value was determined at a wavelength of 540 nm using a microplate reader (Thermo Varioskan, USA). The relative lipid content was calculated using the following equation: relative LD content (%) = sample OD/mg × mL −1 protein.

Immunofluorescence staining and confocal microscopy
For the photographic documentation of LD localization in goose GCs, live cells were stained with 4,4-difluoro-3a,4a-diaza-s-indacene (Bodipy 500/510, Thermo Fisher), wheat germ agglutinin (WGA, Alexa Fluor 594 conjugate, Invitrogen), and Hoechst dyes (33258, Invitrogen) for labeling the LDs, membrane, and nucleus, respectively. Briefly, cells on culture plates were rinsed with PBS and stained for 5 min with 2 μg/mL Bodipy 500/510. After staining, cells were washed three times with PBS for 5 min. The cells were then stained for 5 min with 5 μg/mL WGA and 5 μg/mL Hoechst 33,342. After staining and incubation, the cells were washed thrice with PBS for 5 min each. Finally, cells were fixed in formaldehyde for 10 min. The fixed cells were imaged using a confocal laser scanning microscope (FV1200, Olympus) mounted on an IX83 inverted microscope (Olympus). For treblecolor images, a 550-nm laser line was used to image cells stained with Bodipy 500/510, a 635-nm laser line was used to image cells stained with WGA, and a 543-nm laser line was used to image cells stained with Hoechst 33258.

Sample preparation and extraction
For lipidomic analysis, each group of cells (LS, LG, LN, and two independent siRNA groups [siRNA-210 and siRNA-470, referred to as LT and LF, respectively], as well as the LC) were placed in liquid nitrogen for 2 min, then thawed on ice for 5 min and vortexed. The cells were centrifuged at 12,000 rpm at 4 °C for 10 min. The supernatants (300 µL) were homogenized with a 1 mL mixture containing methanol, MTBE, and an internal standard. This mixture was stirred for 2 min, followed by the addition of 500 µL water, stirring for 1 min, and centrifugation at 12,000 rpm at 4 °C for 10 min. A total of 500 µL of the supernatant obtained was extracted and concentrated. The powder was dissolved with 100 µL mobile phase B (0.1% acetic acid), and the supernatant (200 μL) was transferred to an LC-MS sampling vial with an inner liner, for LC-MS analysis. Quality control (QC) samples were produced by pooling equal aliquots taken from each individual sample in the analytical run together. These QC samples were used to monitor the stability and reproducibility of the analytes in the samples during the analysis.
For RNA-sequencing analysis, each group of cells was extracted from three independent biological replicates using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. RNA quality was confirmed on 1% agarose gels. RNA purity and concentration were determined using a NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) and Qubit ® RNA Assay Kit with a Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA) according to the manufacturer's instructions. RNA integrity was assessed using an RNA 6000 Nano Assay Kit (Bioanalyzer 2100 system; Agilent Technologies, CA, USA).

Lipidomics profiling
Lipidomics profiling was performed by MetWare (Wuhan, China) using an LC-ESI-MS/MS system (UPLC, Shim-pack UFLC SHIMADZU CBM A system, https:// www. shima dzu. com/; MS, QTRAP ® System, https:// sciex. com/). LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (QTRAP ® ) LC-MS/MS System, equipped with an ESI Turbo Ion-Spray interface, operating in positive and negative ion mode, and controlled by Analyst 1.6.3 software (Sciex). Qualitative analysis of primary and secondary MS data was carried out by a comparison of the accurate precursor ions (Q1), product ions (Q3), retention time (RT), and fragmentation patterns with those obtained by injecting standards at the same conditions if the standards were available (Sigma-Aldrich, USA), or was conducted using the self-compiled database MWDB (MetWare, China). The quantitative analysis of metabolites was based on the MRM mode. The characteristic ions of each metabolite were screened using a QQQ mass spectrometer to obtain the signal strengths. Integration and correction of chromatographic peaks was performed using Progenesis QI software (Waters Co., Milford, MA, USA). The corresponding relative metabolite contents were represented as chromatographic peak area integrals. In addition, accurate masses of features representing significant differences were searched against the METLIN, Kyoto Encyclopedia of Genes and Genomes (KEGG), HMDB, and LIPIDMAPS databases.

RNA-sequencing
Sequencing libraries were generated by Novogene (Beijing, China) using a NEBNext ® Ultra ™ RNA Library Prep Kit (Illumina ® , NEB, USA) following the manufacturer's instructions. Clustering of the index-coded samples was performed on a cBot Cluster Generation System by using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq (Illumina, USA) platform and 125 bp/150 bp paired-end reads were generated.

Bioinformatic and statistical analyses of data
Raw LC-MS data were filtered, identified, integrated, corrected, aligned, and normalized using Progenesis QI software (Waters Co., Milford, MA, USA). A data matrix of RT, mass-to-charge ratio, and peak intensity was obtained. The processed data were analyzed using principal component analysis and orthogonal correction partial least squares discriminant analysis (PC) using SIMCA-P14.0 software (Umetrics, Umeå, Sweden). Differentially abundant metabolites between dietary treatments were identified from variable importance in projection (VIP) from OPLS-DA and Student's t tests (VIP > 1 and P < 0.05). Metabolites were identified from public databases, including MassBank (http:// www. massb ank. jp/), KNApSAcK (http:// kanaya. naist. jp/ KNApS AcK/), Human Metabolome Database (http:// www. hmdb. ca/), and METLIN (https:// metlin. scrip ps. edu/). The KEGG database (http:// www. genome. jp/ kegg/) was used to view the enriched pathways of the different metabolites. Hierarchical clustering analysis and heat map analysis were conducted using the R package, version 3.3.1. The raw read counts were normalized considering both the different depths of sequencing among the samples and the gene GC content. Normalization was performed using the EDASeq package. We considered all DEGs with a P value < 0.05 after false discovery rate correction. The putative functions of DEGs were investigated by Gene Ontology (GO) enrichment analysis using the R package GOseq [53], in which gene length bias was corrected. GO terms with a corrected P value < 0.05 were considered significantly enriched. KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways [54]. Statistical plots were prepared using Origin version 6.1.