After Lan et al. (2016); Ji et al. (2018), and Wan et al. (2014) sequenced normal and degenerative NP tissues, respectively, they uploaded circRNA (GSE67566), miRNA (GSE63492/GSE116726), and mRNA (GSE56081) microarray datasets to GEO database. Detailed information for each dataset is shown in Table 1. In terms of GSE67566/GSE63492/GSE116726, the raw data were read and analyzed using the limma package in R (Ritchie et al., 2015), as well as normalized and log2-transformed. By default, the false-positive results of adjusted P-value were corrected using Benjamini and Hochberg false discovery rate (FDR). We identified differentially expressed miRNAs (DEMs) with the criterion of the absolute value of log2 fold change (FC) > 2 and −log₁₀ (FDR) > 2 based on the analysis of GSE116726. The GSE56081 dataset was obtained from Lan et al. (2016) analytical result (FC > 2 or < −2, P < 0.05).

The upstream miRNAs of ESR1 were predicted by Targetscanhuman 7.2² (Agarwal et al., 2015), mirDIP³ (Tokar et al., 2018), starBase⁴ (Li et al., 2014), miRTarBase⁵ (Chou et al., 2018), miRDB⁶ (Chen and Wang, 2020), and miRWalk 3.0⁷ (Dweep and Gretz, 2015) databases, and GSE63492/GSE116726 datasets. Targetscanhuman 7.2 (Agarwal et al., 2015), starBase (Li et al., 2014), miRTarBase (Chou et al., 2018), miRDB (Chen and Wang, 2020), miRWalk 3.0 (Dweep and Gretz, 2015), and miRanda⁸ (John et al., 2004), and GSE56081 dataset were used to predict miR-874-3p targets genes. Additionally, the upstream circRNAs of miR-874-3p were predicted via circbank⁹ (Liu et al., 2019), starBase (Li et al., 2014) databases, and GSE67566 dataset to select IDD-related circRNAs.

Based on the miR-874-3p targets genes predicted by miRTarBase database (Chou et al., 2018), the Cytoscape software was utilized to display these genes (Otasek et al., 2019). Furthermore, Gene Ontology (GO) enrichment analyses was conducted using the Search Tool for the Retrieval of Interacting Genes (STRING)¹⁰ (Szklarczyk et al., 2019), and the predominant enrichment pathways were further visualized by SangerBox tool¹¹. The P < 0.05 was regarded as statistically significant. In addition, protein–protein interaction (PPI) network was constructed, and the degree centrality of the nodes in the PPI network was speculated through the cytoHubba plug-in in Cytoscape software (Chin et al., 2014; Szklarczyk et al., 2019), of which the higher nodes degrees were considered as the hub genes (Zotenko et al., 2008).

The specific method was described in our previous study (Li et al., 2021). Human NPCs were purchased from ScienCell Research Laboratories (ScienCell, Cat. #4800, United States), which were isolated from the NP of human intervertebral disc. NPCs were cultured in Nucleus Pulposus Cell Medium (Cat. #4801, ScienCell, United States) containing 10 mL fetal bovine serum, 5 mL NPC growth supplement, and 5 mL penicillin/streptomycin solution and then were incubated at 37°C in a humidified environment with 5% CO₂. The medium was changed every 2 days. The NPCs were passaged once a week, and well-grown NPCs were taken for subsequent experiments. To simulate the microenvironment of IDD, TNF-α, and IL-1β (10 ng/mL, Proteintech) were employed to stimulate NPCs for 24 h.

The miR-874-3p mimic, miR-874-3p inhibitor, and their corresponding negative controls (NCs) were obtained from JIAMAY BIOLAB (Beijing, China). The empty vector: pcDNA3.1 + Circ Mini (5,607 bp) and overexpression vector: pcDNA3.1 + Circ Mini-circ_0040039 (6,333 bp) and pcDNA3.1 + Circ Mini- circ_0004354 (5,765 bp) were designed and synthesized by HyCell Biotechnology (Wuhan, China). As for NPC transfection, culture plates were incubated at 37°C in a humidified environment with 5% CO₂. CircRNAs plasmids or miR-874-3p mimic or inhibitor or corresponding NCs were transfected into NPCs with Lipofectamine 8000 (Beyotime, China) based on the manufacturer’s protocols. After 48-h transfection, NPCs were collected to conduct the next experiments.

Total RNAs was extracted from NPCs using TRIzol Reagent (Life Technologies, Thermo Fisher Scientific, United States) according to the manufacturer’s protocols. First, 1 μg total RNA and 1 μL Geneseed^® Enzyme Mix (Geneseed, Guangzhou, China) were used to reverse into 20 μL complementary DNA (cDNA) through Geneseed^® II First Strand cDNA Synthesis Kit (Geneseed, Guangzhou, China). Next, 10 μL Geneseed^® quantitative polymerase chain reaction (qPCR) SYBR^® Green Master Mix (Geneseed, Guangzhou, China), 0.5 μL forward (F) primer (10 μM), and 0.5 μL reverse (R) primer (10 μM) were made up and used to conduct quantitative reverse transcriptase (RT)–PCR on ABI7500 system (Applied Biosystems, CA, United States). All specific primers were shown as follows: (1) GAPDH: F1: AGAAGGCTGGGGCTCATTTG, R1: GCAGGAGGCATTGCTGATGAT; (2) ESR1: F2: 5′-ACCCTCC ATGATCAGGTCCA-3′, R2: 5′-AGATCTCCACCATGCCCT CT-3′; (3) miR-874-3p: F3: ATGGTTCGTGGGCTGCCCTGGC, Com R3: GTGCAGGGT CCGAGGT, RT3: GTCGTATCCAG TGCAGGGTCCGAGGTATTCGCACTGGATA.

CGACCtcggtccc; (4) U6: F4: CTCGCTTCGGCAGCACA, R4: AACGCTTCACGA ATTTGCGT, RT4: GTCGTATC CAGTGCAGGGTCCGAGGTATTCGCACTGGATA CGACCAAATATGGAAC. Among them, GAPDH was used as circ_0040039, circ_0004354, and ESR1 control, whereas U6 was used as miR-874-3p control. Their relative expression levels were measured based on the 2^–ΔΔ Ct method described by Livak and Schmittgen (2001).

The well-grown NPCs were inoculated into six-hole cell culture plates at a density of 5 × 10⁵ cells per well. Then, 200 μL diluted RNAs–Lipofectamine 8000 (Beyotime, China) complex was added to the cell wells that had been replaced with 800 μL serum-free medium. The NPCs were then cultured for 0, 1, 2, and 3 days at 37°C incubator. For Cell Counting Kit-8 (CCK8) assay, 10 μL CCK8 solution was added to each well and for incubation for another 1.5 h. The NPC growth was evaluated by CCK8 detection kit according to manufacturer’s protocols (Yeasen, Shanghai, China). The absorbance was determined at OD 450 nm. NPC growth rates were calculated based on the formula: Dayn OD value/Day0 average OD value (same processing sample).

The well-grown NPCs were inoculated into 6 hole cell culture plates at a density of 5 × 10⁵ cells per well. After circ_0040039 or circ_0004354 overexpression vector were transfected into NPCs using Lipofectamine 8000 (Beyotime, China), the NPC apoptosis rates were evaluated by annexin V–APC apoptosis detection kit according to manufacturer’s protocols (keyGEN, KGA1024, China). Annexin V–APC is matched with 7-AAD to distinguish NPCs in different stages of apoptosis. The NPCs were stained with 5 μL annexin V–APC and 5 μL 7-AAD, and then the data were analyzed with FlowJo VX10 software. On the scatterplot of the bivariate flow cytometry (FCM), annexin V + /7-AAD + (Q2) represented the late apoptotic and necrotic NPCs; annexin V + /7-AAD- (Q3) represented the early apoptotic NPCs, whereas annexin V–/7-AAD− (Q4) represented living NPCs.

MiRanda database (John et al., 2004) was used to predict the potential binding sites of miR-874-3p with ESR1 mRNA 3′-UTR. Luciferase reporter vectors: psiCHECK2–Firefly luciferase–Renilla luciferase containing ESR1–700-bp wild-type (WT) sequences or corresponding mutant (MUT) sequences, were constructed by Geneseed Biotech Co. (Guangzhou, China). Human embryonic kidney (HEK) 293T cells were plated on 24-well plates at a density of 1 × 10⁵ cells per well. Subsequently, 1 μg vector and 100 μL miR-874-3p mimic or mimic NC were cotransfected into HEK-293T cells using 2 μL Lipofectamine 8000 (Beyotime, China). After 48-h transfection, the relative luciferase activity was measured using the Dual Luciferase Assay Kit (Promega E1910, Madison, WI, United States) according to the manufacturer’s directions. The activation degrees of the target reporter genes were calculated between different samples according to the obtained ratio of the relative light unit (RLU) value measured by Renilla luciferase divided by the RLU value measured by Firefly luciferase.

The specific method was described in our previous study (Li et al., 2021). RIPA lysis buffer containing phenylmethanesulfonyl fluoride (Beyotime, Shanghai, China) was used to extract the total protein from NPCs. The protein concentrations were quantified using the Micro Bicinchoninic Acid Protein Assay kit (Beyotime, Shanghai, China). After making sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) gels, the proteins were isolated through SDS-PAGE and then transferred to polyvinylidene difluoride (PVDF) membranes (Bio-Rad, CA, United States) at 350 mA for 70 min. Subsequently, the PVDF membranes were blocked by 5% non-fat milk and incubated overnight at 4°C with primary antibody, including anti-ESR1 antibody (diluted 1:1,000; Abcam, ab32063) and anti–β-actin antibody (diluted 1:5,000; Proteintech, 66009-1-Ig), followed by incubation with a secondary antibody. Phosphate-buffered saline with Tween-20 was utilized to wash the PVDF membranes. Finally, the signals were tested by BeyoECL Star Luminescence kit (Beyotime, Shanghai, China) and a chemiluminescence system (Bio-Rad, CA, United States).

All the experiments were performed at least three times. GraphPad Prism software 6 version was used to analyze the data. The statistical significances between the two groups were compared using unpaired Student’s t-test, where the differences among more than two groups were assessed using one-way analysis of variance followed by Tukey multiple-comparisons test. The P < 0.05 was considered as statistically significant.

A large number of studies have indicated the involvement of a multitude of DEMs in regulating IDD via repressing their target genes. However, which target genes were the most important remains unanswered. To address this question, the key differentially expressed genes (KDEGs) of these DEMs were summarized and are listed in Supplementary Table 1. A total of 97 different target genes were reported as IDD-related KDEGs. Given that the STRING website can provide experimental and predicted PPI information, and PPI is the most appropriate tool for studying the potential interrelationship among multiple genes, this study mapped the 97 KDEGs into the STRING website, followed by the analysis of their interaction using cytoHubba plug-in in Cytoscape software. The results unveiled that estrogen receptor α (ESR1) ranked the highest and was a hub gene, suggesting that ESR1 could regulate a series of IDD-related genes, encompassing protective factors, such as SIRT1, Sox9, HIF-1α, and IGF1R, as well as catabolic factors, such as IL-6, MMP2/9, and CASP3 (Figures 1A,B). The Kyoto Encyclopedia of Genes and Genomes analysis of the estrogen signaling pathway revealed that ESR1 might regulate the expression of related genes (Figure 1C) to mediate various signaling pathways, including canonical mitogen-activated protein kinase, PI3K-Akt, and estrogen pathways, thereby affecting cell cycle progression, growth, apoptosis, and other pathological processes (Figure 1D). Furthermore, Cai et al. (2020) demonstrated that the mRNA and protein expression levels of ESR1 were significantly decreased in patients with IDD diseases. Collectively, these results predicted that ESR1 might be one of the most KDEGs in IDD.

MiRNAs can degrade mRNAs and inhibit their translation via directly binding to the 3′-UTR of their target mRNAs (Pasquinelli, 2012; Ji et al., 2018). The upstream miRNAs of ESR1 were predicted and analyzed by bioinformatics analysis. The datasets used in this study were obtained from human NP specimens. Two overlapped IDD-related miRNAs were predicted by merging miRDB, TargetScan, miRTarBase, mirDIP, and miRwalk databases and GSE116726/63492 datasets (Figure 2A). Volcano plots revealed two DEMs in GSE116726; the expression of miR-874-3p was lower than that of miR-130b-3p in IDD (Figure 2B). Conversely, ESR1 was predicted to the target gene of miR-874-3p by intersecting different algorithms, including the GSE56081 dataset (Figure 2C). The potential binding sites of miR-874-3p with ESR1 mRNA 3′-UTR were predicted using the miRanda database (John et al., 2004). An ESR1 fragment with WT or MUT complementary binding sites for miR-874-3p was established and inserted into psiCHECK2 luciferase reporter vectors to confirm further the interaction between miR-874-3p and ESR1 (Figure 2D). MiR-874-3p mimic significantly repressed the luciferase activity of the ESR1-WT vector, whereas such overexpression could not change the activity of the ESR1-MUT vector, revealing that miR-874-3p could directly bind to the 3′-UTR of ESR1 (Figure 2E). Moreover, miR-874-3p mimic repressed and miR-874-3p inhibitor increased the mRNA level of ESR1 (Figure 2F). Thus, miR-874-3p was determined as a key miRNA in this study.

miRTarBase: The experimentally (luciferase reporter assay, Western blot, microarray, and next-generation sequencing experiments) validated miRNA–target interaction database (Chou et al., 2018) was used to predict the target genes of miR-874-3p. The 77 miR-874-3p target genes predicted using the miRTarBase database were then visualized using Cytoscape software (Figure 3A). As shown in Figure 3B, the analysis result of PPI revealed that ESR1 was the most key hub gene of the miR-874-3p target genes. Subsequently, GO functional enrichment analysis for these target genes was conducted using the STRING website, which predominantly included three aspects: biological process (BP), molecular function (MF), and cellular component (CC). The bubble diagram (Figure 3C) and GO chord diagram (Figure 3D) of the GO analysis of BP indicated that miR-874-3p might be mainly involved in regulating different signaling pathways through targeting ESR1/signal transducer and activator of transcription 3 (STAT3)/poly(ADP-ribose) polymerase 1 (PARP1)/cyclin-dependent kinase 9 (CDK9)/histone deacetylases 1 (HDAC1), encompassing negative regulation of macromolecule metabolic process, negative regulation of gene expression, negative regulation of cellular metabolic process, and cell population growth. MF included organic cyclic compound binding, heterocyclic compound binding, RNA binding, nucleic acid binding, and single-stranded RNA binding (Figure 3E). The most enriched in CC were nucleoplasm, intracellular organelle, membrane-bound organelle, cytoplasmic ribonucleoprotein granule, and protein-containing complex (Figure 3F). These results indicated that miR-874-3p might modulate NPC growth and apoptosis to mediate IDD through binding to ESR1 or other mRNAs in the cytoplasm.

Accumulating evidence has uncovered that miRNAs can be adsorbed or even repressed by circRNAs via a ceRNA-dependent mechanism (Piwecka et al., 2017; Cheng et al., 2018; Wang et al., 2018; Xie et al., 2019; Xiang et al., 2020). The upstream circRNAs of miR-874-3p were predicted and analyzed to explore further the novel unidentified circRNAs affecting miR-874-3p function. Circbank is a comprehensive database of human circRNAs containing beyond 140,000 annotated circRNAs from different sources, which can be used to predict the upstream circRNAs of miRNAs (Liu et al., 2019). Eight overlapped circRNAs were predicted by intersecting starBase and circbank databases and the GSE67566 dataset (Figure 4A). Furthermore, a circRNA-miR-874-3p-ESR1 interaction network was constructed and visualized using Cytoscape software (Figure 4B). The two most upregulated IDD-related circRNAs circ_0040039 and circ_0004354, both derived from the syntrophin β2 gene, were predicted to bind to miR-874-3p together (data not shown). Thus, circ_0040039 and circ_0004354 were selected for further investigation. Subsequently, the empty vector (Figure 4C) and overexpression vector of circ_0040039 (Figure 4D) and circ_0004354 (Figure 4F) were constructed. The expression levels of circ_0040039 (Figure 4E) and circ_0004354 (Figure 4G) significantly increased after transfecting their overexpression vector into NPCs. Unexpectedly, both circ_0040039 and circ_0004354 elevated (but not repressed) the miR-874-3p expression level; the role of circ_0040039 was more significant (Figure 4H). In addition, circ_0004354 slightly elevated ESR1 expression without any statistically significant difference, whereas circ_0040039 remarkably repressed the ESR1 expression level (Figure 4I). Figure 4J displayed that circ_0040039 was the most upregulated circRNA in IDD group through the analysis of GSE67566. Furthermore, Western blotting assay demonstrated that circ_0040039 inhibited the protein expression level of ESR1 (Figure 4K). Taken together, these data suggested that circ_0040039 might regulate the miR-874-3p-ESR1 pathway via a stabilization mechanism rather than a canonical ceRNA mechanism, as previously reported (Piwecka et al., 2017; Chen et al., 2019).

Considering that the elevated expression of TNF-α and IL-1β is a hallmark trait during NPC degeneration (Risbud and Shapiro, 2014; Fontana et al., 2015; Oichi et al., 2020; Wang et al., 2020), many researchers used them to simulate the microenvironment of IDD in vitro (Cheng et al., 2018; Wang et al., 2018, 2020). Consistent with the predicted result, the expression level of circ_0040039 significantly increased in proinflammatory cytokine–treated NPCs (Figure 5A). On the contrary, circ_0004354 expression significantly decreased in IL-1β–treated NPCs (Figure 5B). The expression level of miR-874-3p significantly increased but not decreased under the treatment of IL-1β, whereas its expression was not altered in response to TNF-α alone or both TNF-α and IL-1β treatments (Figure 5C). Surprisingly, only IL-1β remarkably enhanced ESR1 mRNA expression, whereas using TNF-α and IL-1β at the same time slightly inhibited its expression without reaching statistically significant differences (Figure 5D). However, this result was inconsistent with previous study. Recently, Song et al. (2021) validated that TNF-α can inhibit ESR1 expression in NPCs. Another literature has indicated that ESR1 silencing can elevate IL-1β and TNF-α expression (Sheng et al., 2018). The difference in experimental results may be related to the state of the NPCs and the different experimental conditions. We cannot rule out that IL-1β might act as a buffer to transiently enhance miR-874-3p and ESR1 expression, thereby delaying the development of IDD. The specific mechanisms of TNF-α and IL-1β do not affect or even increase ESR1 mRNA levels in NPCs, and its biological significance needs future investigation to elucidate. Based on these results, IL-1β was used to simulate the IDD microenvironment for further investigation.

CCK8 and FCM detection assays were performed in circ_0040039-overexpressing NPCs to validate the biofunction of circ_0040039 in NPCs. Compared with empty vector and control, circ_0040039 significantly promoted NPC apoptosis (Figures 6A,B) and repressed NPC growth (Figure 6C) in response to 20 ng/mL IL-1β treatments. Given that circ_0040039 promoted miR-874-3p (Figure 4H) but repressed ESR1 expression (Figures 4I,J), and miR-874-3p repressed ESR1 expression in NPCs (Figure 2D), it was speculated that circ_0040039 promoted NPC degeneration possibly via activating the miR-874-3p-ESR1 signaling pathway.