UBE2E3 regulates cellular senescence and osteogenic differentiation of BMSCs during aging

GEO gene expression data

Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) was created by NCBI, containing numerous gene expression data from research institutions worldwide. We downloaded two gene expression datasets (GSE35956 and GSE35958) from the GEO database. The expression profiling arrays of GSE35956 were generated through the application of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, including five samples of human mesenchymal stem cells from osteoporosis patients and the other five bone marrow from non-osteoporotic donors after total hip arthroplasty. Besides, the expression profiling arrays of GSE35958 were also generated through the application of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, including five osteoporosis samples and four control samples.

Data processing

The original CEL files downloaded from GEO were normalized by R software using R package of “limma” (Ritchie et al., 2015). Based on the platform annotation file, the probe ID was then replaced by the corresponding gene symbol. After calculating the missing value by the KNN method, we used the limma R package to screen each data set for DEGs, by limiting the —log2fold change (FC)— > 2 and p-value < 0.01. The volcano plots and heatmaps of DEGs were generated afterwards. Finally, we intersected DEGs from the two data sets by using Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/) for subsequent analysis.

GO and KEGG pathway analysis of DEGs

We used the “clusterProfiler” package (Yu et al., 2012) in R to analyze and visualize the function of genomic coordinates. Gene ontology (GO) term analysis is a method which contains three parts: biological process (BP), cell component (CC), and molecular function (MF), while Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway is a bioinformatics Database which analyzes gene functions and enriched pathways. The GO term and the KEGG pathway analyses were performed for DEGs with the p-value and the q-value cutoff set to 0.05.

Protein–protein interaction (PPI) network construction and hub gene identification

We predicted the PPI network from the shared DEGs using the online Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org/),which was designed for protein–protein interaction network analyses. Moreover, Cytoscape was used to visualize the constructed PPI network. The top3 most important modules were obtained on Molecular Complex Detection (MCODE) app with the criteria degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2. In the end, we selected hub genes according to the top 10 nodes ranked by degree through the use of cytoHubba.

Expression and co-expression analysis of UBE2E3

The expression of UBE2E3 was downloaded from the Genotype Tissue Expression (GTEx) (https://www.gtexportal.org/) and analyzed in different normal tissues. Furthermore, several osteogenesis related genes including RUNX2, COL1A1, BMP2, FOXP1, ALPL, BGLAP and TAZ were also downloaded. We then calculated the Pearson correlation (r) and p-value in R, and analyzed the expression of UBE2E3 in the single-cell sequencing data between young and old rat (Ma et al., 2020).

Animal

The animal study was conducted under the approval of Xiangya Hospital (Central South University) Ethics Committee. Xiangya Hospital of Central South University Ethics Committee approved this research (201703272). C57BL/6 mice were obtained from the Department of laboratory Animals of Central South University, and kept under the specific-pathogen-free (SPF) (2-month mice, N = 12; 15-month mice, N = 8). They were housed under controlled temperature, and with adequate food and water and were only used to isolate BMSCs. All animals were euthanized to obtain specimens, and we followed the AVMA guidelines in the procedures of euthanizing animals.

Cell culture and transfection

BMSCs were flushed out from tibia and femur of 2-month-old (young BMSCs) and 15-month-old (old BMSCs) mice and cultured in α-MEM with 10% fetal bovine serum (Gibco), 100 units/ml penicillin and 100 µg/ml streptomycin at 37 °C, 5% CO₂ humidity environment. The siUBE2E3 and siNC (Ribibio, Guangzhou, China) were transfected into young BMSCs (from 2-month mice) through the use of lipofectamine RNAiMAX (Invitrogen) according to manufacturer’s instructions. Besides, mUBE2E3 pcDNA3.1-HA-C (Youbio Biological Technology Co., Ltd) were also transfected into old BMSCs (from 15-month mice) with lipofectamine 2000 (Invitrogen) to induce UBE2E3 expression.

qRT-PCR analysis

Trizol reagent (Accurate Biotechnology) was used to extract total RNA from cultured cells for mRNA expression analysis. The reverse transcription kit (Accurate Biotechnology) was used to reverse transcribe 1 µg RNA into first strand cDNA. Finally, we used SYBR Green PCR Master Mix (Accurate Biotechnology) to perform qRT-PCR. Primer sequences can be seen in Table S1.

Western blot

Western blot analysis was used to measure UBE2E3 expression. Firstly, total cell lysates were separated by sodium dodecyl sulfate (SDS)-polyacrylamide gel electrophoresis (PAGE) and transferred to polyvinylidene fluoride (PVDF) membranes. Next, membranes were blocked at room temperature with 5% defatted milk for 1 h. Then the membrane was incubated with anti- UBE2E3 (bs-8352R, 1:500; Bioss), anti- GAPDH (10494-1-AP,1:5000; Proteintech) overnight at 4 °C. After being washed three times with tris-buffer saline and Tween (TBST) for 10 min each, the membranes were incubated with corresponding horseradish peroxidase-conjugated second antibodies. Immunoreactivity was detected by enhancing chemiluminescence reaction.

SA-β-gal staining and observation of senescent cells

Cellular senescence-associated β-galactosidase (SA-β-gal) staining was conducted according to the instructions of the β-gal staining kit (Solarbio Science & Technology). Briefly, after PBS wash, BMSCs were fixed with 4% paraformaldehyde for 15 min, and incubated overnight at 37 °C with staining solution.

Osteogenic differentiation assay

BMSCs were cultured in 12-well plates at 2.5×10⁵ cells/well with osteogenic differentiation media containing 10mM β-glycerol phosphate, 0.1uM dexamethasone, and 50 µM ascorbate-2-phosphate. Alkaline phosphatase (ALP) staining, ALP activity detection (Beyotime Biotechnology) and Alizarin Red staining (ARS) (Cyagen Biosciences Inc) were conducted according to manufacturer’s instructions.

Immunofluorescence

Immunofluorescence was conducted as previous reported instructions (Plafker & Plafker, 2015). Firstly, cells were fixed with 4% paraformaldehyde at room temperature for 15 min. Then the cells were blocked with 5% Bovine Serum Albumin (BSA) for 60 min. After aspirating the blocking buffers, we added the diluted primary antibody (Cell Signaling Technology, 12721, 1:200) to incubate at 4 °C overnight. Then the cells were incubated with Alexa488 nm-conjugated secondary antibodies. Finally, the nucleus was stained with Dapi.

Statistical analysis

All data were displayed as mean ± SD. SPSS 22.0 was used for statistical analysis. Two-tailed Student’s t-tests were conducted for two group comparisons. Statistical significance was considered at P < 0.05 (∗P < 0.05, ∗ ∗P <0.01, ∗ ∗∗P < 0.001).

Identification of DEGS in osteoporosis

We selected two gene expression profiles (GSE35956 and GSE35958) associated with osteoporosis from the GEO database. After normalizing the profiles, we first conducted gene differential analysis by comparing bone marrow of osteoporosis patients and non-osteoporotic donors (Figs. 1A, 1B, 1E, 1F). By limiting —log FC— ≥ 2 and P value < 0.01, a total of 883 downregulated DEGs were identified from GSE35956. In another GSE35958 database, 411 DEGs were identified, including seven genes upregulated and 404 genes downregulated. The volcano plots and heatmaps showed the differential expression of these genes from GSE35956 (Figs. 1C, 1D) and GSE35958 (Figs. 1G, 1H) separately. A total of 176 shared DEGs intersected by the two DEGs sets were then obtained. In general, 176 downregulated DEGs were considered to be associated with osteoporosis for subsequent analysis, while no DEGs were upregulated.

Functional enrichment analyses of DEGs

In view of the small number of DEGs selected above, hardly can significant enrichment results be obtained. In that case, we have appropriately relaxed the limitation of the selection of DEGs as —log FC— ≥ 1 and P value < 0.05. After intersecting the separate DEGs from GSE35956 and GSE35958, we obtained 4,684 DEGs to conduct biological annotation (Fig. 2A). The GO term enrichment analysis results were composed of: BP, CC, and MF (Fig. 2B). The BP result showed significant enrichment in RNA splicing, proteasomal protein catabolic process and RNA splicing, via transesterification reactions. In the CC group, the DEGs were mainly enriched in nuclear speck, spliceosomal complex as well as mitochondrial inner membrane. Moreover, as for the MF, the DEGs were enriched in ubiquitin-like protein transferase activity, catalytic activity, acting on RNA, helicase activity and so on. We also performed the KEGG pathways analysis which identified that the downregulated DEGs were mainly enriched in herpes simplex virus 1 infection, pathways of neurodegeneration-multiple diseases, amyotrophic lateral sclerosis, and Alzheimer disease (Fig. 2C).

Construction of the PPI network and identification of hub genes

In order to predict the interactive relationships among DEGs, we constructed the PPI Network of the 176 DEGs by using the STRING database. After screening, we included a total of 89 nodes and 131 edges in the PPI Network, which was presented by the Cytoscape software (Fig. 2D). Moreover, we then obtained the top three most important modules by using MCODE app in Cytoscape (Fig. 2E). We selected ElonginC (TCEB1), Cell Division Cycle 23 (CDC23), UBE2E3, Cullin-2 (CUL2), Anaphase-promoting complex subunit 10 (ANAPC10), Deoxyuridine 5′-triphosphate nucleotidohydrolase, mitochondrial (DUT), E3 ubiquitin-protein ligase (HERC2), Kelch repeat and BTB domain-containing protein 6 (KBTBD6), DNA topoisomerase 1 alpha (TOP1) and RNA-binding protein 25 (RBM25) as hub genes according to the top10 nodes ranked by degree by using cytoHubba (Fig. 2F). Among these 10 hub genes, we selected UBE2E3 with the following properties: (a) UBE2E3 was highly expressed in the bone marrow; (b) UBE2E3 was positively associated with osteogenesis related genes; (c) UBE2E3 expression reduced in old BMSCs compared with that in young BMSCs; (d) UBE2E3 is a kind of highly conserved metazoan enzyme that pairs with E3 ligases to conjugate Mono-Ub on substrates (Nguyen et al., 2014); and (e) Studies have shown that UBE2E3 has an close relationship with cellular senescence (Plafker et al., 2008; Plafker et al., 2018). Since osteoporosis has a close relationship with age, we selected UBE2E3 as the promising gene that might be associated with osteoporosis.

Validation of UBE2E3 expression

To investigate whether UBE2E3 plays a part in bone formation, we explored the expression of UBE2E3 in 7,862 normal tissues of GTEx. Among the 30 tissues, we found that UBE2E3 was highly expressed in the bone marrow (Fig. 3A). Furthermore, we analyzed the co-expression of UBE2E3 with osteogenic related genes including RUNX2, COL1A1, BMP2, FOXP1, ALPL, BGLAP and TAZ in normal tissues from GTEx (Fig. 3B). The correlations of UBE2E3 with all the selected osteogenesis related genes were positive, suggesting that UBE2E3 may have a close relationship with bone formation. COL1A1 had the largest coefficient of association with UBE2E3 (Fig. 3C). Besides, analysis of single cell sequencing of rat BMSCs showed that UBE2E3 expression was significantly lower in BMSCs from old rats compared with that from younger ones (Figs. 3D, 3E) (Ma et al., 2020). In addition, we found that UBE2E3 significantly decreased in BMSCs of patients with osteoporosis from GSE35956 and GSE35958 (Fig. 3F). Similarly, qRT-PCR and western blot analyses showed that UBE2E3 level was significantly lower in BMSCs from old mice compared with that from young mice, which indicated its lower expression with the development of age (Figs. 3G, 3H). Of note, as the number of passages increased, the expression of UBE2E3 also gradually decreased in mice BMSCs (Figs. 3I, 3J).

UBE2E3 regulates senescence and osteogenic differentiation of BMSCs

We first flushed out and cultured BMSCs from bone marrow of 2-month-old mice, and further knocked down UBE2E3 with siRNA, whose effect was verified by qRT-PCR and western blot (Figs. 4A, 4B). When verifying the effect of UBE2E3 on senescence of BMSCs, we found that 6 days post-siRNA administration, the number of SA-β-gal positive cells in UBE2E3 knockdown group increased significantly compared with the control group (Fig. 4C), and the transcription level of p16 (Omori et al., 2020; Xiao et al., 2020) and p21 (Zhang et al., 2019), markers of cellular senescence, also increased (Fig. 4D). In order to explore the effect of UBE2E3 on the osteogenic differentiation of BMSCs, we induced osteogenic differentiation of BMSCs in vitro. The results showed that knockdown of UBE2E3 significantly inhibited the osteogenic differentiation of young BMSCs as analyzed by ALP staining and ARS results (Fig. 4E). Furthermore, the mRNA levels of the osteogenic differentiation related genes, RUN2, ALP, SP7 and BGLAP, decreased with UBE2E3 knockdown (Fig. 4F).

Likewise, we overexpressed UBE2E3 in old BMSCs with plasmid, which was verified by qRT-PCR and western blot (Figs. 5A, 5B). We found that the number of SA-β-gal positive cells in UBE2E3 overexpression group decreased significantly (Fig. 5C), and the transcription level of p16 and p21 also reduced (Fig. 5D). On the other hand, the upregulation of UBE2E3 significantly promoted the osteogenic differentiation of old BMSCs as analyzed by ALP and ARS staining (Fig. 5E). Additionally, the mRNA levels of osteogenic differentiation related genes increased in UBE2E3 overexpression group (Fig. 5F). Collectively, UBE2E3 regulated the senescence and osteogenic differentiation of BMSCs.

UBE2E3 regulates cellular senescence and osteogenic differentiation through controlling Nrf2 distribution and activity

Kendra et al. found that UBE2E3 could regulate the nuclear translocation and activity of Nrf2 (Plafker & Plafker, 2015). Moreover, studies showed that Nrf2 played an important role in regulating cellular senescence and osteogenic differentiation (Fei et al., 2021; Jie et al., 2020; Yoon, Choi & Lee, 2016; Yuan et al., 2021). Hence, we examined the correlation of UBE2E3 and Nrf2 downstream genes: NQO1, GCLC, GCLM in GTEx, and interestingly found they are positively correlated (Figs. 6A–6D). In order to verify whether UBE2E3 affects Nrf2 nuclear translocation in BMSCs, we used immunofluorescence to observe the cellular distribution of Nrf2 between the control group and the siUBE2E3 group. We found that in siUBE2E3 cells, the nuclear NFR2 was significantly reduced (Fig. 6E). Furthermore, we found that the transcription of Nrf2 downstream target genes: GCLC, GCLC and NQO1 reduced 3 days after knockdown of UBE2E3, which suggested the reduced Nrf2 activity (Fig. 6F). In general, these results suggested that UBE2E3 might regulate cellular senescence and osteogenic differentiation of BMSCs through controlling Nrf2 distribution and activity (Fig. 7).