Mouse SV40-MES-13 GMCs were purchased from BNCC Biological Technology and cultured in DMEM (Beijing Solarbio Science & Technology Co., Ltd.) supplemented with 10% FBS (Biological Industries Sartorius AG) and penicillin/streptomycin at 37˚C in a humidified, 5% CO₂ atmosphere. GMCs stimulated with 3.0 µg/ml LPS (MilliporeSigma) (10-11) were defined as the model group (n=3) and unstimulated GMCs were defined as the control group (n=3).

Cell proliferation rates were assessed using the Cell Counting Kit-8 (CCK-8) assay (cat. no. BB19071X; Shanghai Besto Biological Technology Co., Ltd.). GMCs were seeded (5.0x10⁴ cells/well) into 96-well tissue culture plates with DMEM medium and treated with CCK-8 solution (10 µl/well) for 1 h. The optical density (absorbance at 450 nm) was assessed using an RT-6100 ELISA reader (Rayto Life and Analytical Sciences Co., Ltd.).

To ensure the quality of paired-end sequencing reads, FastQC software (version 0.10.1) was used to evaluate the quality of the original sequencing data. The input reads were considered good data quality when the Q20 base percentage in ‘Reads’ was ≥90% and the Q30 base percentage in ‘Reads’ was ≥80%.

High-throughput RNA-seq was performed by Genesky Bio-tech Co., Ltd. RNA fragmentation was performed using Bioruptor Pico (cat. no. B01060001; Diagenode SA) sonication in RNase-free water. RNA integrity was detected by denaturing gel electrophoresis and quantified using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Inc.). The purified RNA fragments were then used to construct libraries using the TruSeq RNA Sample Prep Kit (cat. no. RS-122-2002; Illumina, Inc.). Libraries underwent quality control and were quantified using an Agilent 2100 bioanalyzer system (Agilent Technologies, Inc.). Paired-end 150 bp sequencing was performed using an Illumina HiSeq 2500 (Illumina, Inc.).

The expression levels of NONMMUG036949.2, NONMMUG089165.1, NONMMUG030447.2, NONMMUG028702.2, NONMMUG039651.2 and NONMMUG032587.2 were assessed using RT-qPCR. Total RNA was extracted from the cells using TRIzol^® (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Primers (Table I) were designed based on the cDNA sequences, which were accessed from the National Center for Biotechnology Information (NCBI) and tested using NCBI Primer BLAST, using Primer Premier 5 (Premier Biosoft International), and synthesized by Sangon Biotech Co., Ltd. The concentration and purity of the isolated RNA were determined using an OD1000+ Ultra Micro Spectrophotometer (WuYi Technologies), and reverse transcription was performed using the PrimeScript™ RT Reagent Kit with gDNA Eraser (cat. no. RR047A; Takara Biotechnology Co., Ltd.). qPCR was performed using a StepOne Plus fluorescence quantitative PCR instrument (Applied Biosystems; Thermo Fisher Scientific, Inc.) with SYBR Green qPCR Master Mix (cat. no. G3322-05; Wuhan Servicebio Technology Co., Ltd.). β-actin was used as the internal control; for quantitative results, the expression of lncRNA was expressed as fold change using the 2^-ΔΔCq method (12).

LncTar (http://www.cuilab.cn/lnctar) was used to predict lncRNA-mRNA interactions utilizing free energy minimization (13-14). LncTar utilizes a variation on the standard ‘sliding’ algorithm approach to calculate the normalized binding free energy (ndG) and predicts the minimum free energy joint structure. ndG ≤0.1 was regarded as the cutoff to determine the paired RNAs as interacting, ndG >0.1 was considered to indicate that paired RNAs did not interact. Subsequently, Cytoscape 3.8.1 was used to visualize the lncRNA-mRNA network after screening the target mRNAs of differentially expressed lncRNAs.

Interactions between miRNAs and lncRNAs were predicted using the miRNA target prediction software, miRanda (http://www.miranda.org) (15), and the three miRNAs of the highest confidence that can bind to lncRNAs were chosen. The target mRNAs of miRNAs were predicted using TargetScan 8.0 (https://www.targetscan.org/vert_80/) (16), and mRNAs with high confidence (cumulative weighted context score cutoff level <-0.8) bound to the miRNAs were screened out (17). Cytoscape 3.8.1 was then used to delineate the lncRNA-miRNA-mRNA ceRNA network.

To evaluate the role of differentially expressed lncRNAs in LPS-induced GMCs, GO biological process (BP) term enrichment and KEGG pathway enrichment were assessed (18). GO terms with P<0.05 were considered statistically significant. KEGG pathway enrichment analysis was used to identify significantly enriched signal transduction pathways or metabolic pathways (P<0.05). GO and KEGG enrichment analysis were performed using the free online data analysis platform OmicShare tools (https://www.omicshare.com/tools).

The PPI network of the mRNA involved in the ceRNA network was generated using the STRING online database tool (Version 11.5, https://string-db.org) (19); interaction pairs with a confidence value >0.4 were deemed significant and were retained. The top 10 genes in the PPI network were evaluated using the MCC, MNC and Degree algorithms through the CytoHubba plug-in for Cytoscape (20). The final hub genes were identified as those that were identified by all three of the algorithms.

The public datasets (GSE104066) utilized in this study were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo) and were initially screened based on disease (CGN), organism (Homo sapiens) and experiment type (expression profiling by array). Ultimately, the GSE104066 dataset was selected for use in CGN analysis.

Statistical analysis was performed using SPSS 22.0 software (IBM Corp.), and data are reported as the mean + SD. The data were subjected to one-way ANOVA with Tukey's post hoc multiple comparison tests. P<0.05 was considered to indicate a statistically significant difference.

The present study evaluated the proliferative ability of GMCs induced by different concentrations of LPS (0.5, 1.0, 3.0, 5.0 and 10.0 µg/ml) using the CCK-8 assay at 24 and 48 h following treatment. The relationship between cell proliferation, LPS concentration and intervention time were assessed (Fig. 1). The results of the CCK-8 proliferation assay demonstrated that LPS could effectively induce cell proliferation in GMCs, with the highest absorbance (optical density 450 nm) observed at a concentration of 3.0 µg/ml (Table SI), which indicated that 3.0 µg/ml was the optimal concentration of LPS for the induction of cell proliferation in GMCs. Furthermore, there was no significant difference demonstrated between the treatment times of 24 and 48 h (P=0.62) at a concentration of 3.0 µg/ml LPS. Consequently, the concentration of 3.0 µg/ml and the treatment time of 24 h were selected for use in subsequent experiments.

In the present study, six groups of cells were used in subsequent experiments, including three groups of LPS-induced GMCs as the model group (LPS1-3) and three groups of normal GMCs as the control group (CON1-3).

After deduplication, quality trimming and quality filtering, the sequencing data at both ends of the R1 and R2 paired reads were assessed as being of good quality (Q20 base and Q30 base were both >95%). The proportion of ‘clean reads’ retained after cleaning was >95% for each of the six samples which met the pre-determined quality requirements for sequencing. The quality control results of the sequencing data were presented (Table II; Table SII).

The violin plot demonstrated the relative abundance of lncRNAs in each sample (Fig. 2A). The violin plot demonstrated that there were significant differences in the expression of lncRNAs in LPS-induced GMCs compared with control GMCs. The threshold for differentially expressed lncRNAs was set at absolute fold change |(FC)|≥1.5 and P<0.05 (21-22). A total of 1,532 differentially expressed lncRNAs, including 594 upregulated lncRNAs and 938 downregulated lncRNAs, were identified from a total of 46,879 lncRNAs (Fig. 2B; Table SIII). A heat map of the top 50 differentially expressed lncRNAs was presented (Fig. 2C). The top 10 upregulated and downregulated lncRNAs are presented in Table III, and the radar map in Fig. 2D presents the top 10 upregulated and downregulated lncRNA expression levels in LPS-induced GMCs compared with control GMCs.

To evaluate the functions of differentially expressed lncRNAs, 556 target mRNAs of the 236 lncRNAs were predicted using LncTar (Table SIV). The highly coordinated expression between lncRNAs and target mRNAs may be due to complementary base pairing between lncRNA and mRNA (23). The lncRNA-mRNA regulatory network presented in Fig. 3 shows the interaction relationship between lncRNA and mRNA. For example, lncRNA NONMMUG029023.2 may regulate Stoml2 mRNA expression. Ptdss2 may be affected by both lncRNA NONMMUG039651.2 and lncRNA NONMMUG095401.1.

To elucidate the biological functions of differentially expressed lncRNAs in LPS-induced GMCs, GO BP term enrichment and KEGG pathway analyses were performed on all mRNAs in the lncRNA-mRNA regulatory network (Fig. 4). Within the GO BP classification, ‘cellular macromolecule metabolic process’, ‘negative regulation of biological process’ and ‘nitrogen compound metabolic process’ were the top three over-represented terms (Fig. 4A). Within the GO cellular components (CC) classification, ‘intracellular’, ‘intracellular part’ and ‘membrane-bounded organelle’ were the top three over-represented terms (Fig. 4B). Within the GO molecular function (MF) classification, ‘galactoside 2-α-L-fucosyltransferase activity’, ‘α-(1,2)-fucosyltransferase activity’ and ‘chromatin insulator sequence’ binding were the top three over-represented terms (Fig. 4C).

The KEGG metabolic pathway enrichment analysis demonstrated that ‘GnRH secretion’, ‘Melanoma’, ‘Ras signaling pathway’, ‘Glycosphingolipid biosynthesis-globo and isoglobo series’ and ‘β-Alanine metabolism’ were the top five most significantly enriched KEGG pathways (Fig. 4D). In addition, inflammation-related signaling pathways, such as the ‘Rap1 signaling pathway’ and ‘Ras signaling pathway’, and substance metabolism pathways, such as ‘Ribofiavin metabolism’, ‘β-Alanine metabolism’ and ‘Histidine metabolism’ were also significantly enriched. These results suggested that inflammation and substance metabolism disorder may be the underlying CGN pathogenesis.

To verify the reliability of the sequencing data and to support more clinically meaningful follow-up research, RT-qPCR validation of certain highly conserved lncRNAs between humans and mice was performed. The conservation information of lncRNAs was obtained from the NONCODE database (http://www.noncode.org). The conservation analysis of lncRNAs was assessed using the E-value, where a threshold of E-value <1x10^-5 was used (24). After conservation analysis among species (Table IV), NONMMUG036949.2, NONMMUG089165.1, NONMMUG030447.2, NONMMUG028702.2, NONMMUG039651.2 and NONMMUG032587.2 were selected for RT-qPCR validation (25). The RT-qPCR results demonstrated the same expression trend as the RNA-Seq results, with the six selected lncRNAs all significantly upregulated in LPS-induced GMCs compared with the control group (Fig. 5; Table SV).

In addition to the lncRNA-mRNA regulatory network, lncRNAs can also serve a role in other gene networks that regulate diverse biological processes, such as the ceRNA network (26). To construct an lncRNA-associated ceRNA network, the aforementioned six highly conserved and RT-qPCR-validated lncRNAs were used to construct ceRNA networks. The ceRNA network consisted of miRNAs of the highest confidence with lncRNAs, and mRNAs with a high confidence (cumulative weighted context score cutoff level <-0.8) (27) bound to the miRNAs, which included 6 lncRNAs, 18 miRNAs and 419 mRNAs (Fig. 6; Tables SVI and SVII). lncRNAs can be competing targets of shared miRNAs with other mRNAs and form a complex regulatory ceRNA network. For example, the lncRNA NONMMUG089165.1 may act as a sponge for mmu-miR-3960 to affect EVX1 expression. PARD3B could be affected by both the lncRNA NONMMUG039651.2/mmu-miR-7081-5P axis and the lncRNA NONMMUG028702.2/mmu-miR-7044-5P axis at the same time.

To evaluate the biological functions of the lncRNAs in the ceRNA network, GO BP and KEGG pathway analysis was performed on all mRNAs involved in the ceRNA network (Fig. 7). Within the GO BP classification, ‘system development’, ‘animal organ development’ and ‘multicellular organism development’ were the top three over-represented terms (Fig. 7A). Within the GO CC classification, ‘neuronal cell body membrane’, ‘cell body membrane’ and ‘α-β T cell receptor complex’ were the top three over-represented terms (Fig. 7B). Within the GO MF classification, ‘DNA-binding transcription activator activity, RNA polymerase II-specific’, ‘α-1,6-mannosylglycoprotein 6-β-N-acetylglucosaminyltransferase activity’ and ‘DNA-binding transcription factor activity’ were the top three over-represented terms (Fig. 7C). The KEGG metabolic pathway enrichment analysis demonstrated that ‘Collecting duct acid secretion’, ‘Circadian rhythm’, ‘Proteoglycans in cancer’, ‘Pathways in cancer’ and ‘Oxytocin signaling pathway’ were the top five most significantly enriched KEGG pathways (Fig. 7D). Furthermore, inflammation-related signaling pathways, such as ‘MAPK signaling pathway’, ‘mTOR signaling pathway’ and ‘Rap1 signaling pathway’ were significantly enriched KEGG pathways. These results suggested that inflammation may be the underlying CGN pathogenesis.

The PPI network of all the mRNAs involved in ceRNA was successfully constructed, including 396 nodes and 260 edges (Fig. 8A). The top 10 genes were evaluated in the PPI network using the MCC, MNC and Degree algorithms; the results of the three algorithms were cross-compared to identify five core genes, including IKAROS family zinc finger 1 (IKZF1), CD3 ε subunit of T-cell receptor complex (CD3E), synapsin 1 (SYN1), glutamate ionotropic receptor NMDA subtype 2B (GRIN2B) and SPI1) (Fig. 8B). According to the ceRNA hypothesis, lncRNAs compete with miRNAs, to bind to the target genes, which increases their mRNA expression. Public dataset GSE104066 from the GEO database was used to validate the expression of IKZF1, CD3E, SYN1, GRIN2B and SPI1 in CGN. The data demonstrated that IKZF1, CD3E, SYN1 and GRIN2B were significantly upregulated; however, no significant difference was demonstrated for SPI1 (Fig. 8C; Table SVIII).