Identification and Validation of Autophagy-Related Genes in Primary Ovarian Insufficiency by Gene Expression Profile and Bioinformatic Analysis

Background To investigate the relationship between primary ovarian insufficiency and autophagy, we detected and got the expression profile of human granulosa cell line SVOG, which was with or without LPS induced. The expression profile was analyzed with the focus on the autophagy genes, among which hub genes were identified. Results Totally, 6 genes were selected as candidate hub genes which might correlate with the process of primary ovarian insufficiency. The expression of hub genes was then validated by quantitative real-time PCR and two of them had significant expression change. Bioinformatics analysis was performed to observe the features of hub genes, including hub gene-RBP/TF/miRNA/drug network construction, functional analysis, and protein-protein interaction network. Pearson's correlation analysis was also performed to identify the correlation between hub genes and autophagy genes, among which there were four autophagy genes significantly correlated with hub genes, including ATG4B, ATG3, ATG13, and ULK1. Conclusion The results indicated that autophagy might play an essential role in the process and underlying molecular mechanism of primary ovarian insufficiency, which was revealed for the first time and may help to provide a molecular foundation for the development of diagnostic and therapeutic approaches for primary ovarian insufficiency.


Introduction
Primary ovarian insufficiency (POI), also known as premature menopause or premature ovarian failure (POF), is characterized by cessation of menstruation before the expected age of menopause [1] and is a frequent cause of female infertility. POI is typically defined as at least a 4-month history of oligomenorrhoea and elevated plasma levels of follicle-stimulating hormone (FSH; >25 IU/L) [2]. Numerous complications are found in POI patients, consisting of not only infertility but also sexual dysfunction, vasomotor symptoms, Alzheimer's disease, osteoporosis, and cardiovascular diseases [3]. Respectively, it affects approximately 1 in 100 women at the ages of 40 while 1 in 1000 women at the age of 30 [4].
POI is highly heterogeneous in etiology, such as genetic, autoimmune, iatrogenic, and infectious factors. However, the exact causes of POI remain unknown despite of its strong genetic link, indicating that genes likely to be associ-ated with this condition are yet to be discovered [5]. Currently, the medical treatment of POI is mostly by estrogen supplementation, which has some side effects such as increasing risk of endometrial carcinoma and breast cancer [6]. Hence, it is essential to explore underlying mechanism of POI, which can contribute to elucidate the role of genes in pathogenesis and generate applicable insights to develop preventive strategies, novel diagnostic modalities, and targeted therapeutics. Lipopolysaccharide (LPS) has been revealed to play an essential role in causing hormonal imbalance, ovarian dysfunction, and even infertility [7]. It can be utilized to induce primary ovarian insufficiency (POI) in mouse through causing fibrosis, inflammation, and granulosa cell apoptosis [8]. The exposure to LPS can lead to cell programed death of granulosa cells such as apoptosis and pyroptosis, which contribute to steroidogenesis dysfunction and reproduction failure finally [9]. However, the relationship between LPSinduced POI and autophagy has not been identified.
Autophagy, a lysosomal degradative process of damaged cytoplasmic organelles or cytosolic components, occurs in all eukaryotic cells from yeast to mammals and takes part in numerous physiological processes such as differentiation, development, and aging and contributes to both innate and adaptive immunity [10]. Autophagy was proved to be involved in various diseases based on genetic studies, consisting of inflammatory diseases, neurodegenerative diseases, autoimmune disorders, and cancers. Autophagy could play essential roles in multiple types of tumor through affecting tumor growth, regulating tumorigenesis, giving an adaptive response to cancer cells, and so on [11]. For instance, autophagy participated in breast cancer through regulating tumor immune response [12], involved in uveal melanoma by regulating the expression of Beclin-1 [13,14], and plays important roles in ovarian cancer through causing cisplatin resistance via mediating ERK [15]. Autophagy is active in cells in a basal state but can also be induced in response to multiple forms of cellular stress, such as nutrient or growth factor deprivation, reactive oxygen species, hypoxia, protein aggregates, DNA damage, damaged organelles, or intracellular pathogens [16]. Basal autophagy is important for maintaining cell function by controlling the quality of proteins and organelles, while under stressful conditions, autophagy plays the principal role to supply nutrients for survival [17]. Macroautophagy starts when a double membrane cisterna envelops cytosolic material, consisting of organelles and proteins, expanding into a vesicle called the autophagosome which then fuses with the endosomal-lysosomal system to form an autolysosome [18]. Multiple autophagy proteins take part in controlling this process [19]. In the ovary, autophagy act as a cell survival mechanism, which is involved in maintaining the endowment of female germ cells prior to establishing primordial follicle pools [20]. Autophagy has also been reported to be involved in follicular atresia, which is cell and developmental stage specific [21].
The studies mentioned above strongly indicated that autophagy variants could also originate primary ovarian insufficiency in human. However, to our knowledge, there is little known about the relationship between POI and autophagy. In our study, we detect the human granulosa cell line (SVOG), which was with or without LPS induced. The expression profiles were analyzed, with an especial focus on the autophagy genes and their function. Moreover, the expression of hub genes was validated by quantitative realtime PCR (qRT-PCR).

Study Design.
In order to illustrate the data preprocessing, analysis, and validation, a schematic flow diagram of the study is presented in Figure 1.
2.2. RNA-Sequence Analysis. In our research, the human granulosa cell line (SVOG) was divided into LPS treatment group and normal control group with four samples in each group. The LPS treatment group was performed with LPS, while the normal control group was performed with PBS. Total RNA isolation was performed using RNeasy Plus Mini Kit (QIAGEN), which was followed with RNA-seq. Raw reads in fastq format were analyzed with perl and finally turned into read count of mRNA.
2.3. Identification of DEGs. According to the above groups based, we performed normalization and differential gene expression analysis using the "DESeq2" R package. The normalization was based on the "Relative Log Expression" method, which is specifically implemented in the "DESeq2" [22]. The scaling factors were calculated using the median ratio between gene abundances and the geometric mean. As a method for differential analysis of transcriptome count data, DESeq2 improves the interpretability and stability of estimation because of shrinkage estimators for fold change (FC) and dispersion [22]. The differential gene expression analysis was conducted. Then, we specified logFC > 0:5 and P value < 0.05 as upregulated genes while logFC < −0:5 and P value < 0.05 as downregulated genes, and both upregulated genes and downregulated genes were defined as differentially expressed genes (DEGs). "ggplot2" and "ComplexHeatmap" R package [23] were used to generate the volcano plot and heatmap.

Identification and Functional Analysis of ARGs.
Search for the word "autophagy" on the GeneCards (https://www .genecards.org/) to retrieve autophagy-related genes (ARGs), which were downloaded and intersected with DEGs. The intersection was autophagy DEGs.
The GO biological processes were shown in three aspects, including biological processes (BP), cellular components (CC), and molecular functions (MF) [24]. P value < 0.05 as the cut of criterion was considered statistically significant [24]. GO enrichment analysis of autophagy DEGs was implemented by clusterProfiler package [25], and items were considered as significantly different if they meet the conditions of P value < 0.05 (BH method). R package GOplot [26] was used to integrate the quantitative information by implementing high-quality and novel plotting, which provided us with a collection of multilayered and prespecified charts. Valuable information was added to each layer to display the intended message.

PPI Network Construction and Module Analysis.
STRING (https://string-db.org/) is a website about protein interaction, whose aim is to achieve a comprehensive and objective global network and present them with a unique set of computational predictions [27]. In order to explore the mutual relationship between proteins encoded by different genes, autophagy DEGs were imported into STRING website for further analysis. Next, we output the analysis results to a TSV format file and used Cytoscape software (version 3.8.0) for detail processing and module analysis. CytoHubba [28] is a plug-in downloaded from Cytoscape App Store, which can find hub genes in PPI. Therefore, we applied this plug-in to detect top 10 hub genes in PPI network with the process of MMC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, and ClusteringCoefficient. After filtering and visualizing with the Upset diagram, STRING website was used again for achieving PPI network of hub genes.
RNAInter was used to facilitate elucidating the role of RNAbinding protein (RBP), whose score was more than 0.15. hTFtarget (http://bioinfo.life.hust.edu.cn/hTFtarget) provides a comprehensive, reliable, and user-friendly resource for exploring human TF-target regulations, which is useful for a wide range of users in the transcription factors (TF) and gene expression regulation community [30]. KnockTF (http://www.licpathway.net/KnockTF/index.html) constructs a TF-differentially expressed gene network and performs network analyses for genes of interest, which can help elucidate TF-related functions and potential biological effects [31]. We predicted transcription factors through hTFtarget and KnockTF with the criterion of logFC ≤ −1 and network was constructed according to those TFs.
ENCORI (The Encyclopedia of RNA Interactomes, http:// starbase.sysu.edu.cn/index.php) showed extensive and complex RNA-RNA and protein-RNA interaction networks by analyzing a large set of Ago and RBP binding sites derived from all available CLIP-Seq experimental techniques (PAR-CLIP, HITS-CLIP, iCLIP, and CLASH) [32]. The Drug-Gene Interaction Database (DGIdb, http://www.dgidb.org) is a web resource that provides information on drug-gene interactions and druggable genes from publications, databases, and other web-based sources [33]. They were both used to explore the miRNA interaction networks and potential drugs. The results were visualized by Cytoscape.

Validation of Autophagy and Expression of Hub Genes.
Samples of human granulosa cell line SVOG, consisting of the normal control group and LPS-treated group, were collected and lysed in RIPA lysis buffer (PC101, EpiZyme, China) and then proteins were extracted. The BCA protein assay kit (NCM Biotech, China) was employed to detect the protein concentration. The proteins (15 μg) were loaded and electrophoresed on 12.5% SDS-polyacrylamide gels and then transferred to PVDF membranes. After blocking with protein-free rapid blocking buffer (PS108P, EpiZyme, China) for 30 min at room temperature, the PVDF membranes were probed with primary antibodies against LC3 (1 : 2000, CST) or P62 (1 : 2000, CST) at 4°C overnight. And the anti-β-actin antibody (1 : 5000, Abclonal, China) was used as the internal control. Then, the membranes were incubated with secondary antibodies for 90 min at room temperature. The ECL western blotting kit (NCM Biotech, China) was applied to visualize the targeted protein bands. We used the ImageJ software to evaluate the protein band densities.
Total RNA was isolated from the human granulosa cell line SVOG, which was with or without LPS induced, using TRIzol (RNAiso Plus; Takara, Japan) according to the manufacturer's instructions. We also obtained RNA samples from another human granulosa cell line KGN in order to avoid deviation. One microgram of RNA was reversely transcribed into cDNA using a PrimeScript™

Analytical Cellular Pathology
RT reagent kit (Takara). Amplification was performed using TB Green® Premix EX Taq™ II (Takara) and gene-specific primers (Sangon, Shanghai, China) on a qRT-PCR device (QuantStudio 5, Thermo Fisher Scientific, Waltham, MA, USA). β-Actin was used as an internal control. The relative expression of the genes was calculated using the 2-ΔΔCT method.
2.9. Statistical Analysis. The statistical analyses were performed using R 4.0.2 (https://www.r-project.org/, version Volcano plot (a)   Analytical Cellular Pathology 3.6.3). All the statistical tests were 2-sided, and P < 0:05 was considered statistically significant. In addition, results were analyzed statistically using the Mann-Whitney U test and Wilcoxon t-test (P < 0:05).

Identification of Differentially Expressed Genes.
After standardization of microarray results, the differentially expressed genes (DEGs) were identified, and the results illustrated that there were 29 upregulated genes and 46 downregulated genes (Figures 2(a)-2(c)). Then, we downloaded autophagy genes from GeneCards. The overlap among the autophagy genes and DEGs contained 22 genes as shown in the Venn diagram (Figure 2(d)).

GO Enrichment Analyses of DEGs.
To analyze the biological classification of DEGs, functional and pathway enrichment analyses were performed, and the results are shown in Table 1. GO analysis results showed that changes in biological processes (BP) of DEGs were significantly enriched in cGMPmediated signaling, regulation of tube diameter, regulation of blood vessel size, regulation of blood vessel diameter, regula-tion of tube size, vascular process in circulatory system, negative regulation of smooth muscle cell proliferation, positive regulation of secretion, and regulation of synaptic vesicle cycle (Figures 3(a) and 3(d)-3(f)). Changes in cell component (CC) of DEGs were mainly enriched in transport vesicle, presynapse, neuron projection terminus, clathrin-coated vesicle, excitatory synapse, exocytic vesicle, transport vesicle membrane, late endosome, coated vesicle, ruffle membrane, synaptic vesicle membrane, and exocytic vesicle membrane (Figure 3(b)), while changes in molecular function (MF) were mainly enriched in amine binding, serotonin binding, syntaxin-1 binding, neurotransmitter binding, syntaxin binding, and ammonium ion binding (Figure 3(c)).

PPI Network Construction and Hub Gene Selection.
The PPI network of autophagy DEGs was constructed (Figure 4(a)), which illustrated that there were 22 autophagy DEGs correlated with POI as well as 22 protein-protein interactions. CytoHubba was used to filtrate data for hub genes (Figure 4(b)). A total of 6 genes were identified as hub genes with significance in all 11 arithmetic, which indicated that they might play important role in the process of POI (Table 2).  Public databases were used to explore potential RBP, TF, and miRNA and predict potential therapeutic agents, which might possibly interact with hub genes. The results indicated that DKC1 was the most potential gene according to RNAInter database ( Figure 5(a)). The hTFtarget database and KnockTF database were used to predict the potential transcription factors ( Figure 5(b)). The results revealed that FOXP1 could activate PTPRN, PRKG1, and SLC2A4, while GATA1 was able to activate SLC2A4, SLC6A4, and SYT1. Furthermore, regulatory microRNAs (miRNAs) were predicted for hub genes and used to establish a potential hub gene-miRNA regulation network using the ENCORI platform ( Figure 5(c)). We uploaded hub genes to the DGIdb database, and all drug options have been approved, among which there were two kinds of drugs that could target two hub genes simultaneously, namely, HALO-PERIDOL and COCAINE. The hub gene-drug network is constructed and shown in Table 3.

Correlation between Hub Genes and Autophagy Genes.
To explore the impact of hub genes on POI, we analyzed the correlation between hub genes and autophagy genes by Pearson's correlation (Figure 6(a)). The results indicated that there was notable significant correlation between autophagy genes and some hub genes.

Expression Analysis of Hub Genes and Validation of
Autophagy. The differential expression of hub genes was revealed based on FPKM datasets. Compared to normal control, the gene expression changed in LPS-treated group, in which BSN, PTPRN, and SLC6A4 downregulated while PRKG1, SLC2A4, and SYT1 notably upregulated (Figure 7). The results were validated in human granulosa cell lines KGN and SVOG by qRT-PCR, revealing upregulated expression of PRKG1, SLC2A4, and SYT1 and downregulated gene expression of BSN, PTPRN, and SLC6A4, among which the expression changes of SYT1 and SLC6A4 in SVOG cell line were of significance (Figure 8 and Figure S1). The primer sequences of hub genes were shown in Table 4.
After treated with LPS, the expression of proteins related to autophagy significantly changed (Figure 9). The ratio of LC3 II/I increased while the expression of P62 decreased, which both indicated that autophagy was promoted after LPS treatment.

Discussion
POI is characterized by marked heterogeneity, which is with a significant genetic contribution. However, it is challenging to identify exact causative genes because of numerous unduplicated discoveries. It is essential and necessary to take stock of the field, frame the progress, outline the controversies, and evaluate future directions in elucidating the genetics of POI. Multidisciplinary approaches could help to explore molecular signatures [34], which accounts for the progression of disease and contribute to searching for specific treatment [35]. Identification of effective and novel therapy strategies based on multiomics methods can make it possible to have a better management and improve the prognosis of patients [36]. Hence, we applied various approaches to reveal the underlying mechanism of POI, consisting of bioinformatics, cellular biology, and molecular biology.
In our study, we performed RNA sequence analysis between the LPS-treated group and normal control group, whose results were then used to identify the differentially expressed genes. Functional and pathway enrichment analysis and PPI network were performed to analyze the biological classification of DEGs. According to arithmetic in CytoHubba, six genes were identified as hub genes, namely, BSN, PTPRN, SLC6A4, PRKG1, SLC2A4, and SYT1. The expression of hub genes was revealed based on FPKM datasets as well as validated by qPCR, while various characteristics of hub genes were also identified, such as potential RBP, TF, miRNA, drugs, and correlation with autophagy genes.
Compared to normal control, some hub gene expression downregulated in LPS-treated group, including BSN,   Analytical Cellular Pathology

12
Analytical Cellular Pathology increased after LPS treatment, indicating that LPS promoted autophagy in granulosa cells [37]. SYT1 is an important presynaptic vesicle protein which binds Ca2+ to regulate synaptic vesicle exocytosis and was revealed to act as an important regulator in mouse oocyte activation events consisting of cortical granule exocytosis and the generation of Ca2+ signals [38]. SLC6A4 encodes an integral membrane protein, which act as the role to transport the neurotransmitter serotonin from synaptic spaces into presynaptic neurons. SLC6A4 5HTTLPR polymorphism was revealed to be related to insulin blood levels and insulin secretion during OGTT in patients with polycystic ovary syndrome (PCOS) [39]. Recent research has implied that BSN acted in concert to control presynaptic autophagy [40]. PTPRN encodes a member of the protein tyrosine phosphatase (PTP) family, which is known to be signaling molecules that regulate a variety of cellular processes including cell growth, differentiation, mitotic cycle, and oncogenic transformation.  13 Analytical Cellular Pathology PTPRN, serves as an essential negative regulator for insulin signaling, whose ablation was reported to protect against ER stress-induced cardiac anomalies through regulation of autophagy [41]. Protein kinase cGMP-dependent 1(PRKG1) was involved in cGMP-mediated signaling [42], which was also revealed in the biological process enrichment. The activity of PRKG1 in follicular cells was demonstrated to be essential for oocyte maturation [43] and was proved to play important role in the glucose uptake of granulosa cells mediated by nitric oxide [44]. The expression of SLC2A4 might be suppressed in PCOS women, which may consequently alter endometrial function [45].   CCACCCAATTCCGAGTATGGT  PRKG1  GGACAGGACTCATCAAGCATAC  CTTCACGAGTGACATTTACCGTT  SLC6A4  TGACACACGGCACTCTATCC  AGCCAATCACTGAGAGAAGGA  PTPRN  TTGAGCATGACCCTCGGATG  GCCAGAAGTCTGCGATGGTAT  BSN  GCCCTCTATCCACCAAGGC  GTCTTGCTGGGTTCAGAAGC  SLC2A4 GCCATGAGCTACGTCTCCATT GGCCACGATGAACCAAGGAA 14 Analytical Cellular Pathology To demonstrate whether POI was correlated with autophagy, Pearson's correlation was utilized to explore between hub genes and autophagy genes. It was found that the expression levels of BSN had a high negative correlation with ATG4B and the expression levels of SYT1 also had a high negative correlation with ULK1. High expression levels of SLC2A4 were correlated with high expression levels of ATG3, while high SLC6A4 expression was related with low expression levels of ATG13. Therefore, we demonstrated that there might be various correlations between POI and autophagy, which lead to the notable expression changes of autophagy genes.
For patients with genetic risk of POI, novel methods aiming at fertility preservation can benefit them a lot through early diagnosis. Recent research, for instance, revealed that there was close correlation between POI and the mutations of BRCA [46]. BRCA can bind to various essential regulatory proteins and regulate gene expression [47], which make it possible for BRCA to play a novel role in the regulation of autophagy [48]. The assessment of BRCA mutations could be detected by next-generation sequencing or droplet digital PCR [49]. Novel techniques applied for disease early diagnosis could help patients to get a better management and prognosis. It is important to perform further study to identify the correlation between gene regulation and disease, while it is equally essential to search for and validate the application of novel techniques.
However, there were some limitations in our study. Firstly, the pathogenesis of POI is multidimensional, and it is not convincing enough to demonstrate POI only in the aspect of autophagy. Secondly, more clinical characteristics of POI patients should be included in subgroup analysis and future study. Thirdly, though we validated the expression of hub genes in mRNA expression level through qRT-PCR, it is necessary to validate the protein expression and perform more multicenter and prospective studies to evaluate the possible applications of molecular signatures in the future. In addition, further studies containing in vivo and in vitro experiments are required to elucidate the molecular mechanisms of hub genes for clinical applications.

Conclusion
In conclusion, we performed RNA sequence analysis to extract the DEGs, based on which we constructed PPI network and identified six hub genes. The expression of hub genes was not only revealed based on FPKM datasets but validated by qPCR, while numerous characteristics of hub genes were also identified, consisting of potential RBP, TF, miRNA, drugs, and relationship with autophagy genes. The results indicated that autophagy might play an essential role in the process and underlying molecular mechanism of POI.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.