text:
Expanded CUG Repeat RNA Induces Premature Senescence in Myotonic Dystrophy Model Cells

Myotonic dystrophy type 1 (DM1) is a dominantly inherited disorder due to a toxic gain of function of RNA transcripts containing expanded CUG repeats (CUGexp). Patients with DM1 present with multisystemic symptoms, such as muscle wasting, cognitive impairment, cataract, frontal baldness, and endocrine defects, which resemble accelerated aging. Although the involvement of cellular senescence, a critical component of aging, was suggested in studies of DM1 patient-derived cells, the detailed mechanism of cellular senescence caused by CUGexp RNA remains unelucidated. Here, we developed a DM1 cell model that conditionally expressed CUGexp RNA in human primary cells so that we could perform a detailed assessment that eliminated the variability in primary cells from different origins. Our DM1 model cells demonstrated that CUGexp RNA expression induced cellular senescence by a telomere-independent mechanism. Furthermore, the toxic RNA expression caused mitochondrial dysfunction, excessive reactive oxygen species production, and DNA damage and response, resulting in the senescence-associated increase of cell cycle inhibitors p21 and p16 and secreted mediators insulin-like growth factor binding protein 3 (IGFBP3) and plasminogen activator inhibitor-1 (PAI-1). This study provides unequivocal evidence of the induction of premature senescence by CUGexp RNA in our DM1 model cells.

INTRODUCTION
Myotonic dystrophy type 1 (DM1, OMIM #160900) is the most common form of muscular dystrophy in adults (Johnson et al., 2021). The clinical manifestations of DM1 are characterized by multisystemic symptoms, including muscle wasting, myotonia, cardiac conduction defects, cataracts, frontal balding, cognitive impairment, and endocrine defects (Thornton, 2014;Meola and Cardani, 2015). DM1 is caused by the expansion of CTG repeats in the 3′-untranslated region (UTR) of DMPK (Brook et al., 1992;Fu et al., 1992;Mahadevan et al., 1992). The RNA expressed from the mutant allele exerts a toxic gain of function due to the presence of an expanded CUG repeat (CUG exp ), which forms ribonuclear foci in the nucleus and disrupts the regulation of alternative splicing by affecting several RNA-binding factors, including muscleblind-like splicing regulator one and CUGBP Elav-like family member one proteins (Lee and Cooper, 2009). Aberrant splicing events in CLCN1, SCN5A, and INSR have been associated with myotonia, arrhythmias, and insulin resistance, respectively, in DM1 (Savkur et al., 2001;Charlet et al., 2002;Freyermuth et al., 2016). However, although more than hundreds of mis-splicing events have been identified in DM1 tissues (Nakamori et al., 2013;Freyermuth et al., 2016;Otero et al., 2021), the aberrant splicing events responsible for many multisystemic DM1 symptoms remain unclear, suggesting the involvement of additional factors in DM1 pathogenesis.
Some of the multisystemic symptoms in DM1, such as muscle wasting, cataract, cognitive impairment, and frontal baldness, resemble the appearance of accelerated aging (Meinke et al., 2018). Cellular senescence is a crucial driver of the aging process and is defined as a state of irreversible cell cycle arrest induced by stress or specific physiological processes, such as telomere erosion, oncogene overexpression, oxidative stress, mitochondrial dysfunction, and inflammation (Hernandez-Segura et al., 2018;Di Micco et al., 2021). It is characterized by morphological and metabolic changes, chromatin reorganization, and a senescence-associated secretory phenotype (SASP) (McHugh and Gil, 2018). Several studies have shown that DM1 patient-derived cells reduce proliferative capacity and induce cellular senescence (Furling et al., 2001;Bigot et al., 2009;Thornell et al., 2009;Renna et al., 2014). Other studies have suggested the involvement of oxidative stress and mitochondrial dysfunction in DM1 pathogenesis (Sahashi et al., 1992;Usuki and Ishiura, 1998;Siciliano et al., 2001;Toscano et al., 2005;Garcia-Puga et al., 2020). However, the mechanism of cellular senescence in DM1 remains still unclear, because of the limitations in assessing cellular senescence in human primary cells. The characteristics of patient-derived cells vary considerably in different genetic and environmental backgrounds (Chang et al., 2002;Maier and Westendorp, 2009;Dolivo et al., 2016). Furthermore, the proliferative lifespan of human primary cells depends on the donor's age (Kaji et al., 2009). In this study, we developed a cell model of DM1 conditionally expressing abnormal RNA containing CUG exp to overcome the issue of variability in primary cells of different origins. Then, we investigated whether CUG exp RNA induced cellular senescence and which senescence process was key in DM1 pathogenesis.

Cell Culture and Transfection
Plasmid pLC16 was used for conditional transcription of expanded CTG repeats, as previously described (Nakamori et al., 2011b). Briefly, pLC16 consists of a cytomegalovirus/ chicken β-actin enhancer/promoter, followed by a floxed selection-stop cassette, a downstream complementary DNA (cDNA) sequence for hygromycin resistance, and, finally, the human DMPK 3′-UTR, modified with restriction sites for insertion of expanded CTG repeats (Supplementary Figure  S1A). The selection-stop cassette contains cDNA encoding a puromycin resistance protein (puro) followed by a triple-stop transcription terminator. The expanded CTG repeat was synthesized in the repeat donor plasmid pDWD using cell-free cloning by amplification of dimerized expanded repeats. The expanded CTG repeat was inserted into the DMPK 3′-UTR sequence in the pLC16 plasmid construct.
DM1 fibroblasts (GM05281) were purchased from Coriell Institute. DM1 myoblasts were obtained from biopsy specimens as described previously (Nakamori et al., 2017). Human fetal lung fibroblasts (TIG-3 cells, passage 20) were obtained from the Japanese Collection of Research Bioresources Cell Bank. They were cultured in Dulbecco's modified Eagle medium with low glucose (Gibco) supplemented with 10% fetal bovine serum. Then, the TIG-3 cells were cotransfected with pLC16 containing 800 CTG repeats and phiC31 integrase using Nucleofector Technology with program U-023 (Lonza). Stably transfected cells were selected with puromycin (0.5 μg/ml). The expression of expanded CUG RNA was induced by Cre recombinase-mediated excision of the puro-transcription-terminator cassette. Cells with recombination were selected using hygromycin B (50 μg/ml). The proliferative capacity of the cells was assessed by calculating CPDL, as previously described (Bigot et al., 2009).

Fluorescent in situ Hybridization (FISH)
FISH was performed as previously described (Nakamori et al., 2011b). The resultant fluorescence images were obtained using a BZ-X710 fluorescence microscope (Keyence).

SA-β-Gal and BrdU Activity and Proliferative Capacity
SA-β-gal activity was determined using a Senescence Cells Histochemical Staining Kit (Sigma-Aldrich), according to the manufacturer's protocol. Images were obtained with a BZ-X710 fluorescence microscope. Proliferating cells were detected using a BrdU Immunohistochemistry Kit (Abcam), according to the manufacturer's protocol. Briefly, the cells were incubated with 10 μM BrdU for 24 h at 37°C. BrdU-positive cells were analyzed using ImageJ software (National Institutes of Health).

Quantitative RT-PCR
Total RNA was extracted from the DM1 model cells using an RNeasy Mini Kit (Qiagen). Then, the RNA was reverse transcribed to cDNA using a Superscript III First-Strand Synthesis System (Invitrogen), according to the manufacturer's protocol. qPCR was performed using TaqMan Gene Expression assays on an ABI PRISM 7900HT Sequence Detection System (Applied Biosystems). Gene expression was determined using TaqMan primers and probes for CDKN1A, TP53, CDKN2A, IGFBP3, SERPINE1, NQO1, HMOX1, MMP1, and MMP3 (Applied Biosystems). Relative mRNA expression was normalized to 18S rRNA. The level of endogenous DMPK plus transgene-derived mRNA was determined as described previously (Nakamori et al., 2011a). Alternative splicing of MBNL1 exon 5 and MBNL2 exon 5 was analyzed as described previously .

RNA-Seq Analysis
Whole transcriptome RNA-seq analysis was performed on CUG exp− and CUG exp+ cells using a NovaSeq 6,000 System (Illumina). The raw data were evaluated by NGQC software (Novogene) and trimmed to remove adaptor contaminants and low-quality reads. Clean reads were aligned to the reference genome sequence using TopHat v2.0.12 (Trapnell et al., 2009). The estimated transcript abundance was calculated, and the count values were normalized to the upper quartile of the fragments per kilobase of transcript per million mapped reads using HTSeq v0.6.1 (Anders et al., 2015). GO enrichment analyses were conducted using goseq v2.12.0 (Young et al., 2010) and Metascape (http://www.metascape.org). Sequencing data have been deposited in Gene Expression Omnibus under accession number GSE196265.

Enzyme-Linked Immunosorbent Assay for Secreted IGFBP-3
The levels of secreted IGFBP-3 were measured using a Human IGFBP-3 Quantikine ELISA Kit (R&D Systems), according to the manufacturer's instructions. Recombinant human IGFBP-3 protein (R&D Systems) was used as the standard control.

Telomere Length Measurement
Genomic DNA was extracted from the DM1 model cells using a Gentra Puregene Cell Kit (Qiagen). Then, telomere length was measured using a Relative Human Telomere Length Quantification qPCR Assay Kit (ScienCell), according to the manufacturer's protocol. Briefly, the extracted genomic DNA was added to a reaction containing a primer pair (telomere or single copy reference) and Power SYBR Green PCR Master Mix (Applied Biosystems). PCR was performed on an ABI PRISM 7900HT Sequence Detection System (Applied Biosystems).

Comet Assay
Comet assays were performed, as previously described (Olive and Banath, 2006). Briefly, cell suspensions were mixed with a 1% low-gelling-temperature agarose solution and spotted onto agarose-precoated glass slides. The slides were gently submerged in an alkaline lysis solution for 4 h at 4°C and then transferred to an electrophoresis solution and electrophoresed at 20 V for 25 min. Next, the slides were submerged in rinse buffer for 30 min at room temperature and then incubated with 2.5 μg/ ml SYBR Green Ⅰ Nucleic Acid Gel Stain (Invitrogen) for 20 min.
The results were analyzed using ImageJ software and scored as the percentage of tail DNA and Olive tail moment.

Immunofluorescence
The DM1 model cells were plated on Lab-Tek II chamber slides (Thermo Fisher Scientific) and incubated for 24 h. After washing with PBS, the cells were fixed in 4% paraformaldehyde for 15 min and permeabilized with 0.3% Triton for 5 min. Then, cells were blocked with 5% BSA for 30 min and incubated with mouse antiphospho-Histone H2A.X (Ser139) antibody (1:500; 05-636, Merck Millipore) overnight at 4°C. After washing with PBS, cells were incubated with goat-anti-mouse Alexa 488 secondary antibody (1:500) for 1 h. Cells were washed with PBS and mounted with Vectashield Hardset mounting medium with DAPI (Vector Laboratories). The fluorescence images were obtained using a BZ-X710 fluorescence microscope (Keyence).

Abasic Site Measurement
The number of AP sites was assessed using a colorimetric DNA Damage Quantification Kit-AP Site Counting (Dojindo), according to the manufacturer's instructions. The results were measured on a Multiskan FC Microplate Photometer (Thermo Fisher Scientific).

Flow Cytometry Analysis of ROS Production and Mitochondrial Membrane Potential
ROS production was determined by flow cytometry using MitoSOX Red Mitochondrial Superoxide Indicator (Invitrogen), and the mitochondrial membrane potential was measured by flow cytometry using Rhodamine 123 (Invitrogen). Briefly, the DM1 model cells were seeded at a density of 80% confluency on 12-well plates. Then, H 2 O 2 (200 µM) was added to induce ROS production, and the cells were incubated at 37°C for 1 h. After detaching the cells from the wells using trypsin, they were incubated with 5 µM MitoSOX or 10 μM Rhodamine 123 for 30 min at 37°C. Finally, the cells were analyzed on a BD FACSCanto II flow cytometer (BD Biosciences). The total cell population was defined according to the forward versus side scatter dot plot, and data for the live cells only were extracted for analysis. The median fluorescence intensity (MFI) of at least 10,000 cells was analyzed using FlowJo software (BD Biosciences). Unstained cells were used as a control.

ATP Measurement
The levels of cellular ATP were quantified using an Intracellular ATP Assay Kit v2 (TOYO B-Net), according to the manufacturer's instructions. The emitted luciferin luminescence was quantified using a Glomax 20/20 Luminometer (Promega). Results were corrected for protein concentration using ATP-extracted samples.

Statistical Analysis
Data were presented as the mean ± standard deviation (SD) from at least three independent biological replicates. Statistical Frontiers in Genetics | www.frontiersin.org March 2022 | Volume 13 | Article 865811 significance was tested using the Student's t-test. p values <0.05 were considered statistically significant.

Expanded CUG Repeat RNA Induces Cellular Senescence in DM1 Model Cells
We developed a DM1 cell model to investigate the effect of CUG exp RNA on cellular senescence. First, we transfected normal human lung fibroblasts (TIG-3 cells) with plasmid pLC16 containing 800 CTG repeats (day 6). Then, we treated the cells with Cre recombinase, and obtained DM1 model cells conditionally expressing CUG exp RNA after 18 days of selection using hygromycin B (day 47). This cell model allowed us to control CUG exp RNA expression in primary cells with the same genetic background so that we could assess the mechanism of cellular senescence by CUG exp RNA expression more accurately than in patient-derived cells.
In the cells expressing CUG exp RNA following induction by Cre recombinase, the total expression level of transgene and endogenous DMPK was increased 3.8-fold compared with that of endogenous DMPK and DM1 fibroblasts and to a similar level in DM1 myoblasts (Supplementary Figure S1B). RNAfluorescence in situ hybridization (FISH) experiments showed RNA foci formation, as observed in DM1 patient cells ( Figure 1A; Supplementary Figure S1C). Splicing misregulation in MBNL1 and MBNL2 was also induced by CUG exp RNA expression (Supplementary Figure S1D). We calculated the cumulative population doubling level (CPDL) by culturing the DM1 model cells to compare the proliferative capacity of cells expressing CUG exp (CUG exp+ cells) with those not expressing CUG exp (CUG exp− cells). The induction of abnormal RNA slowed cell proliferation and prematurely terminated cell division ( Figure 1B; Supplementary Figure S2A). We also confirmed that the induction of CUG exp RNA reduced the proliferative capacity in two other independently established TIG-3 cell lines conditionally expressing CUG exp RNA (Supplementary Figure  S2B). However, the integration sites of the transgene and the expression level of CUG exp RNAs are different in each cell line. Therefore, to avoid the effect of variation in genetic background, we further investigated the mechanism of the premature arrest of cell proliferation in a representative cell line with a sufficient evaluation period after the induction of CUG exp RNA expression in the same genetic background. Furthermore, we confirmed that Cre induction and hygromycin B selection did not affect cell proliferation in TIG-3 cells transfected with plasmid pLC16 containing no CTG repeat (Supplementary Figure S2C).
We examined the activity of senescence-associated βgalactosidase (SA-β-gal), a biomarker of senescent cells, to clarify whether the premature arrest was associated with senescence. There were 9.4-fold and 4.3-fold increases in SAβ-gal positivity in CUG exp+ cells compared with CUG exp− cells on day 51 and day 63, respectively (p = 0.0008 and 0.00006, respectively, Figure 1C). Furthermore, CUG exp+ cells showed significantly lower levels of bromodeoxyuridine (BrdU) positivity (indicating cells with proliferative activity) compared with CUG exp− cells (p = 0.0083, Figure 1D). Thus, CUG exp expression induced early arrest of cell division and senescence in the DM1 model cells.

Expanded CUG Repeat RNA Alters Gene Expression Profiles Related to Cellular Senescence
We performed RNA sequencing (RNA-seq) analysis in CUG exp+ and CUG exp− cells at day 51 to identify the pathway that induced cellular senescence when CUG exp RNA was expressed. Upregulated and down-regulated genes in CUG exp+ cells had some overlap with those in DM1 myoblasts or myotubes reported in a previous study, despite the cell type difference (Todorow et al., 2021) (Supplementary Table S1). Gene ontology (GO) enrichment analysis was performed to identify relevant biological pathways. Genes for extracellular matrix (ECM) organization, positive regulation of the apoptotic process, and regulation of cell growth were highly enriched among the upregulated genes in CUG exp+ cells (Figure 2A). The genes involved in ECM organization included MMP1 and ADAMT genes ( Figure 2D). These genes encode ECM-degrading enzymes, which are important molecules comprising the SASP (Mavrogonatou et al., 2019). Among the ECM-degrading enzymes, the matrix metalloproteinase (MMP) family is predominantly linked with cellular senescence (Mavrogonatou et al., 2019). For example, MMP1 and MMP3 expression was reported as upregulated in senescent fibroblasts (Coppe et al., 2010). We evaluated MMP1 and MMP3 expression in the DM1 model cells at days 51 and 63 by quantitative reverse transcription-polymerase chain reaction (qRT-PCR) and found Frontiers in Genetics | www.frontiersin.org March 2022 | Volume 13 | Article 865811 6 that the expression of these genes was significantly increased in CUG exp+ cells (p = 0.0087 and 0.023 at day 51, and p = 0.0002 and 0.0007 at day 63, respectively, Figure 3A), suggesting that CUG exp RNA alters the expression of several ECM-degrading enzymes related to cellular senescence. Genes for muscle structure development and response to oxidative stress were significantly enriched among the downregulated genes ( Figure 2B). The genes involved in the response to oxidative stress and antioxidant expression, such as NQO1 and HMOX1, were significantly decreased (p = 0.016 and 0.005 at day 51, and p = 0.0026 and 0.0031 at day 63, respectively, Figures 2D, 3A).
Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that the p53 signaling pathway was highly enriched among the upregulated genes in CUG exp+ cells ( Figure 2C). This pathway involves SERPINE1, encoding plasminogen activator inhibitor-1 (PAI-1), and IGFBP3, a downstream target of PAI-1-induced senescence, which were both markedly upregulated (p = 0.0003 and 0.0014 at day 51, and p = 0.0095 and 0.00002 at day 63, respectively, Figures 2D, 3A). These genes were reported not only as markers of senescence but also as inducers of senescence in human and mouse fibroblasts (Debacq-Chainiaux et al., 2008;Vaughan et al., 2017). The protein levels of PAI-1 and insulin-like growth factor binding protein 3 (IGFBP3) were also significantly increased in CUG exp+ cells at day 54 (p = 0.0026 and 0.021, respectively, Figure 3B). We then assessed the activity of Akt, a critical downstream target downregulated by PAI-1 and IGFBP3 (Balsara et al., 2006;Elzi et al., 2012). Akt phosphorylation at Ser473 was mildly decreased in CUG exp+ cells at day 54, although this was not statistically significant (p = 0.065, Figure 3B). IGFBP3 has also been reported as a secreted mediator of cellular senescence with paracrine and autocrine activity (Vassilieva et al., 2020). The IGFBP3 levels in the culture medium of the DM1 model cells showed a 4-fold increase in CUG exp+ cell cultures at day 54 (p = 0.00013, Figure 3C). These results indicated that the PAI-1-IGFBP3 pathway was activated by CUG exp RNA expression in the DM1 model cells.

Expanded CUG Repeat RNA Activates Cell Cycle Checkpoint Inhibitors
PAI-1 and IGFBP3 were reported as induced by p53 activation (Grimberg et al., 2005;Vaughan et al., 2017). p53 is activated in response to DNA damage and various stressors that characteristically promote cellular senescence to regulate the cell cycle (Di Micco et al., 2021). Cell cycle checkpoint genes, such as TP53 encoding p53, CDKN1A encoding p21, and CDKN2A encoding p16, induce cellular senescence despite the presence of metabolic activity (Hernandez-Segura et al., 2018;McHugh and Gil, 2018). To investigate the influence of CUG exp RNA on cell cycle arrest, we measured the expression levels of cell cycle-related genes by qRT-PCR in CUG exp− and CUG exp+ cells at days 51 and 63. TP53 and CDKN1A were significantly increased in CUG exp+ cells (p = 0.021 and 0.011 at day 51, and p = 0.021 and 0.0009 at day 63, respectively, Figure 3A). CDKN2A was mildly increased in CUG exp+ cells at day 51 (p = 0.076, Figure 3A) and significantly increased at day 63 (p = 0.017, Figure 3A). The protein levels of p21 and p16 were significantly increased in CUG exp+ cells at day 54 (p = 0.018 and 0.00043, respectively, Figure 3B). p53 was also mildly increased in CUG exp+ cells, although the increase did not reach statistical significance (p = 0.33, Figure 3B). p21 and p16 are sufficient to establish cell cycle arrest in an independent and interdependent manner (McHugh and Gil, 2018;Di Micco et al., 2021). Our results indicated that CUG exp RNA expression increased the expression of cell cycle checkpoint inhibitors, resulting in cell cycle arrest in the DM1 model cells.

Expanded CUG Repeat RNA Causes DNA Damage in a Telomere-Independent Manner
Next, we investigated the regulatory mechanism responsible for the increased levels of cell cycle inhibitors p21 and p16 and secreted mediators IGFBP3 and PAI-1 in the DM1 model cells.
Telomere shortening is an important molecular mechanism of cell cycle arrest (Di Micco et al., 2021). Shortened telomeres are recognized as DNA breaks, and they promote cell cycle arrest, known as replicative senescence (Hernandez-Segura et al., 2018;Di Micco et al., 2021). We used qPCR to measure the telomere length in the DM1 model cells at days 43, 49, and 61 to confirm whether the cellular senescence induced by CUG exp occurred in a telomere-dependent manner. Although the telomere length in CUG exp− and CUG exp+ cells was shortened in later passages, there was no significant difference in the telomere length between both cell types ( Figure 4A). This suggested that telomere shortening was not associated with accelerated cellular senescence caused by the toxic CUG exp RNA.
Telomere-independent senescence (premature senescence) occurs when DNA damage response (DDR) is triggered in response to DNA damage, even in normal proliferating cells (Di Micco et al., 2021). We performed a comet assay on the DM1 model cells to evaluate the DNA damage induced by CUG exp . The parameters reflecting the extent of DNA damage, such as the percentage tail DNA and Olive tail moment (OTM), were significantly increased in CUG exp+ cells at day 53 (p = 0.00022 and 0.0056, respectively, Figure 4B). Furthermore, γ-H2AX foci, indicative of DNA damage, were significantly increased in CUG exp+ cells (p = 0.019, Figure 4C). Next, we evaluated apurinic/apyrimidinic (AP) sites (markers of base excision repair) to confirm DDR activation in response to DNA damage in the DM1 model cells. The AP site lesions were increased in later passages and were more significantly increased in CUG exp+ cells at day 49 (p = 0.0015, Figure 4D), suggesting that abnormal CUG exp RNA increased DNA damage and DDR in the DM1 model cells.

Expanded CUG Repeat RNA Impairs Mitochondrial Function and Increases Reactive Oxygen Species
We then investigated the cause of DNA damage in CUG exp+ cells. ROS are important in the induction of cellular senescence because they damage DNA via a variety of mechanisms, including oxidized DNA bases, AP sites, and double-strand breaks (Davalli et al., 2016). Considering the possibility that excessive ROS production may cause the DNA damage by toxic repeat RNA, we measured mitochondrial ROS by flow cytometry. H 2 O 2 , which induces intracellular ROS generation, was added to evaluate the response to oxidative stress (Weng et al., 2018). In our study, mitochondrial ROS were significantly increased in CUG exp+ cells at day 55 under normal conditions, and the difference was more pronounced following H 2 O 2 exposure (p = 0.018 and 0.036, respectively, Figure 5A).
Mitochondria are the major sites of ROS production (Turrens, 2003) and mitochondrial dysfunction contributes to ROS overproduction (Korolchuk et al., 2017). Therefore, we evaluated mitochondrial dysfunction as a cause of ROS production. Rhodamine 123, which accumulates in activated mitochondria that have a high membrane potential, was unaffected by CUG exp RNA expression ( Figure 5B). However, although the mitochondrial membrane potential remained unchanged, intracellular ATP levels were significantly decreased in CUG exp+ cells at day 55 (p = 0.0082, Figure 5C). Decreased cellular ATP production results from mitochondrial bioenergetic dysfunction, and insufficient mitochondrial energy production leads to excessive ROS production (Korolchuk et al., 2017). Thus, our results indicated that CUG exp RNA caused mitochondrial dysfunction, contributing to premature senescence via ROS production.

DISCUSSION
The presentation of multisystemic symptoms that resemble features of aging in patients with DM1, such as muscle FIGURE 4 | Expanded CUG repeat RNA causes DNA damage in a telomere-independent manner. (A) Measurement of relative telomere length by qPCR. Telomere length was measured at three different stages (at days 43, 49, and 61) in cells not expressing CUG exp (CUG exp− ) (blue) and cells expressing CUG exp (CUG exp+ ) (red). Data are presented as means ± SD of three independent experiments. (B) Representative images of DNA damage in comet assays at day 53 (left). DNA damage was evaluated according to the percentage tail DNA and Olive tail moment (right). Data are presented as means ± SD of three independent experiments. *p < 0.01, **p < 0.001. (C) Representative images of DNA damage in γ-H2AX assays at day 53 (left). Bar graph showing percentage of γ-H2AX positive cells (right). Data are presented as means ± SD of three independent experiments. *p < 0.05. (D) Quantification of DNA damage response by AP site measurement. The AP site lesions per 100,000 base pairs were measured in the cells before transfection (day 5) and at day 49 in cells not expressing CUG exp (CUG exp− ) and cells expressing CUG exp (CUG exp+ ). Data are presented as means ± SD of three independent experiments. *p < 0.01.
Although several studies of DM1 patient-derived cells indicate the involvement of cellular senescence (Furling et al., 2001;Bigot et al., 2009;Thornell et al., 2009;Renna et al., 2014), the senescence process remains unelucidated. In this study, we established a primary cell model of DM1 with the conditional expression of CUG exp RNA. We demonstrated that CUG exp expression induced premature senescence in the DM1 model cells. Additionally, our results suggest that premature senescence in the DM1 model cells is associated with mitochondrial dysfunction and ROS production, leading to the upregulation of secreted mediators IGFBP3 and PAI-1 and cell cycle regulators p16, p21, and possibly p53 ( Figure 6). Cellular senescence is defined as the permanent cessation of cell proliferation. It is induced by two distinct mechanisms: replicative senescence (triggered by telomere loss) and premature senescence (caused by a telomere-independent mechanism) (Kuilman et al., 2010). Our results suggest that the expression of CUG exp RNA induces cellular senescence regardless of telomere length, which correlates with the findings of another study reporting reduced proliferative capacity in muscle satellite cells from patients with DM1, even though the telomeres did not reach a critical size (Bigot et al., 2009).
Telomere-independent premature senescence is triggered by various cellular stressors, including oxidative stress, oncogene activation, and DNA damage agents such as ionizing radiation or Quantification of intracellular ATP production by firefly luciferase luminescence at day 55. Data are presented as means ± SD of three independent experiments. **p < 0.01. RLU: relative light units.
FIGURE 6 | Proposed mechanism of premature senescence by CUG exp RNA expression. RNA toxicity by CUG exp causes mitochondrial dysfunction, which contributes to ROS overproduction, triggering DNA damage and DDR. Accumulated DNA damage activates cell cycle regulators, such as p53, p21, and p16. ROS-induced DNA damage upregulates the secreted mediators PAI-1 and IGFBP3 in a p53-dependent and -independent manner, which drive premature senescence. These multiple steps in the process of senescence interact with one another, inducing premature senescence.
Frontiers in Genetics | www.frontiersin.org March 2022 | Volume 13 | Article 865811 chemotherapeutic substances (Kuilman et al., 2010;Hernandez-Segura et al., 2018). Exposure to these stressors induces persistent DNA damage and DDR, resulting in premature senescence. Our study focused on oxidative stress, a typical internal factor that causes premature senescence, because the use of our DM1 cell model allowed us to exclude the influence of external factors. Oxidative stress is characterized by ROS overproduction (Davalli et al., 2016). ROS directly damage DNA, causing the induction and maintenance of cellular senescence (Passos et al., 2010;Davalli et al., 2016). The majority of ROS are produced in mitochondria, and mitochondrial ROS production increases with mitochondrial metabolic dysfunction (Turrens, 2003;Korolchuk et al., 2017). Mitochondrial ROS were slightly but significantly increased in CUG exp+ cells solely by induction of abnormal RNA expression even in physiological condition. Previous studies demonstrated that even slight changes in mitochondrial ROS production are associated with DNA damage and mitochondrial dysfunction in fibroblasts (Kobashigawa et al., 2015;Lofaro et al., 2020). Our findings of increased mitochondrial ROS, decreased intracellular ATP, and DDR activation in CUG exp+ cells indicate that CUG exp RNA promotes ROS-mediated premature senescence by mitochondrial dysfunction. ROS homeostasis is maintained by antioxidants. Nuclear factor-erythroid 2 related factor 2 (Nrf2) transcriptionally regulates several antioxidant genes, such as HMOX1 and NQO1, to decrease ROS levels (Kasai et al., 2020). On the other hand, excess ROS activates p53, which suppresses the Nrf2-regulated antioxidant genes (Drane et al., 2001;Faraonio et al., 2006). Nrf2-dependent antioxidant factors were decreased in our CUG exp+ DM1 model cells, suggesting that Nrf2, which is normally activated to eliminate excessive ROS, is suppressed by p53, exacerbating the vulnerability of our DM1 model cells to oxidative stress. Moreover, oxidative stress has been clinically associated with cataracts (Pendergrass et al., 2005), frontal hair loss (Shin et al., 2013), and metabolic dysfunction (de Almeida et al., 2017), which are commonly observed in patients with DM1. Furthermore, antioxidant capacity has been reported as impaired in DM1 patients' serum (Toscano et al., 2005;Kumar et al., 2014). These findings suggest that ROS-induced premature senescence leads to the features of aging observed in DM1.
The mechanism of cell cycle arrest by CUG exp RNA is largely unknown. The cell cycle is chiefly regulated by the p53-p21 and p16-pRb pathways. Although these pathways have complex interactions, p16 is mainly involved in maintaining cellular senescence by mitogenic stress and sustained DDR (Mijit et al., 2020;Di Micco et al., 2021). Muscle satellite cells derived from patients with congenital DM1 have been reported to induce p16-dependent premature senescence (Bigot et al., 2009). In our study, p16 and p21 expression was significantly increased in CUG exp+ cells at both the transcript and protein levels. Moreover, p53 expression was significantly increased in CUG exp+ cells at the transcript level and mildly increased at the protein level. The p53-p21 pathway, which is necessary for senescence induction, is upregulated following DDR activation (Mijit et al., 2020;Di Micco et al., 2021). Thus, our data indicate that abnormal CUG exp RNA can induce senescence by activating the p53-p21 pathway and cause irreversible senescence by p16 upregulation.
PAI-1 and IGFBP3 are major factors affecting cellular senescence, mainly as downstream targets of p53 (Grimberg et al., 2005;Vaughan et al., 2017;Vassilieva et al., 2020). In our study, the transcript and protein levels of PAI-1 and IGFBP3 were significantly increased in CUG exp+ cells. PAI-1 and IGFBP3 are regulated in both p53-dependent and -independent manner in response to ROS-induced DNA damage (Grimberg et al., 2005;Eren et al., 2014). Hence, CUG exp RNA may increase the expression of PAI-1 and IGFBP3 not only through p53 but also through another signaling pathway activated by ROSinduced DNA damage. Further, PAI-1 and IGFBP3 induce cellular senescence with activation of cell cycle inhibitors. For example, IGFBP3 overexpression upregulates p21 in some cancer cells (Wu et al., 2013), and inhibition of PAI-1 activity reduces p16 expression (Eren et al., 2014). Therefore, the excessive ROS, DNA damage, and DDR observed in CUG exp+ cells may increase PAI-1 and IGFBP3, leading to premature senescence via activation of cell cycle inhibitors. Additionally, Akt is involved in PAI-1-and IGFBP3-induced cell growth suppression and cellular senescence (Balsara et al., 2006;Elzi et al., 2012). We observed a slight reduction in Akt activity in CUG exp+ cells, which is consistent with a previous study using fibroblasts from DM1 patients (Garcia-Puga et al., 2020). Our results suggest the possibility that abnormal CUG exp RNA can inactivate Akt, causing PAI-1-and IGFBP3-mediated senescence. Furthermore, recent studies have shown the suppression of Akt signaling in DM1 skeletal muscle, so reduced Akt activity may be associated with skeletal muscle atrophy (Crawford Parks et al., 2017;Sabater-Arcis et al., 2020;Ozimski et al., 2021). Skeletal muscle atrophy is generally associated with mitochondrial ROS (Powers, 2014), and PAI-1 inhibits the regeneration of damaged skeletal muscle (Rahman and Krause, 2020). Thus, the mechanism of cell proliferation inhibition via PAI-1 and IGFBP3 by ROS-induced DNA damage may also contribute to muscle atrophy in DM1. However, it should be noted that the mechanism of premature senescence observed in our DM1 model fibroblast cells could be different in DM1 myoblast cells. Even in myoblasts, a different mechanism of premature senescence was reported in DM1 and DM2 (Renna et al., 2014).
Previous studies have suggested that cellular senescence in DM1 is associated with mitochondrial dysfunction, oxidative stress, DNA damage, and activation of cell cycle inhibitors (Hasuike et al., 2022). Our study demonstrated the direct effects of abnormal CUG exp RNAs on these factors by using the DM1 model cells conditionally expressing CUG exp . However, we were unable to evaluate either the characteristics of CUG exp+ cells at a pre-senescence stage or the therapeutic effect by targeting CUG exp RNA or ROS, since CUG exp+ cells exhibited delayed cell proliferation soon after the cell model establishment (day 47). Hence, it is not clearly determined whether the senescence inducers such as mitochondrial dysfunction and ROS production observed in CUG exp+ cells are the cause or result of premature senescence. Generally, senescence occurs in a bidirectional manner and via various positive feedback mechanisms. For example, cell cycle regulators, such as p53, p21, and p16, produce ROS as a downstream signal transduction factor without oxidative DNA damage, while ROS itself induces p53 (Davalli et al., 2016). Furthermore, the positive feedback loops of mitochondrial damage, ROS production, and DDR activation via the p53-p21 pathway are necessary and sufficient to maintain cell cycle arrest (Passos et al., 2010). Our data suggest that the activation of these positive feedback loops maintains senescence in CUG exp+ cells. Several studies have shown impaired mitochondrial function in DM1 (Gramegna et al., 2018;Garcia-Puga et al., 2020). Thus, excessive ROS production by mitochondrial metabolic dysfunction may trigger the series of senescence pathways induced by CUG exp RNA.
In conclusion, our findings indicate that abnormal expanded CUG repeat RNA leads to mitochondrial dysfunction, ROS production, and DDR activation, and induces premature senescence. Various drug targets have been investigated to eliminate the influence of senescent cells and SASP (Di Micco et al., 2021). Interventions targeting the senescence-inducing factors resulting from CUG exp RNA may be potential therapeutics for symptoms that resemble accelerated aging in DM1.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, accession ID: GSE196265.

AUTHOR CONTRIBUTIONS
MN designed the research. MN and YH performed the research, analyzed the data, and wrote the manuscript. HM supervised the research. All authors contributed to the article and approved the submitted version.


Upregulation of Centromere Proteins as Potential Biomarkers for Esophageal Squamous Cell Carcinoma Diagnosis and Prognosis

Esophageal squamous cell carcinoma (ESCC) has a high incidence and low survival rate, necessitating the identification of novel specific biomarkers. Centromere-associated proteins (CENPs) have been reported to be biomarkers for many cancers, but their roles in ESCC have seldom been investigated. Here, the potential clinical roles of CENPs in ESCC patients were demonstrated by a systematic bioinformatics analysis. Most CENP-encoding genes were differentially expressed between tumor and normal tissues. CENPA, CENPE, CENPF, CENPI, CENPM, CENPN, CENPQ, and CENPR were upregulated universally in the three datasets. Survival analysis demonstrated that high expression of CENPE and CENPQ was positively correlated with the outcomes of ESCC patients. The CENPE-based forecast model was more accurate than the tumor-node-metastasis (TNM) staging-based model, which was classified as stage I/II vs. III/IV. More importantly, the forecast model based on the commonly upregulated CENPs exhibited a much higher area under the curve (AUC) value (0.855) than the currently known TTL, ZNF750, AC016205.1, and BOLA3 biomarkers. The nomogram model integrating the CENPs, TNM stage, and sex was highly accurate in the prognosis of ESCC patients (AUC = 0.906). Besides, gene set enrichment analysis (GSEA) demonstrated that CENPE expression is significantly correlated with cell cycle, G2/M checkpoint, mitotic spindle, p53, etc. Finally, in validation experiments, we also found that CENPE and CENPQ were significantly overexpressed in esophageal cancer cells. Taken together, these results clearly suggest that CENPs are clinically promising diagnostic and prognostic biomarkers for ESCC patients.

Introduction
Esophageal cancer ranks seventh globally in terms of incidence and sixth in overall mortality [1]. Esophageal cancer is classified into two main subtypes: esophageal squamous cell carcinoma (ESCC) and adenocarcinoma (EAC). In parts of Asia and sub-Saharan Africa, nearly 90% of esophageal cancer cases are ESCC [2], suggesting that ESCC is the dominant subtype. Although a variety of diagnostic methods and multiple therapies for ESCC patients have proven to be effective, the 5-year overall survival (OS) rate is still <20% [2]. Several studies have suggested that long noncoding RNA (lncRNA), microRNA-(miR-) 375, and miR-483-5p might serve as potential biomarkers for ESCC [3][4][5]. However, the outcomes of ESCC patients are still not ideal. The identification of novel diagnostic and prognostic biomarkers for ESCC is urgently needed to improve the outcomes of ESCC patients.
Centromere-associated proteins (CENPs) are a group of proteins involved in kinetochore formation. CENPs contain 18 inner kinetochore-located proteins (CENPA, B, C, H, I, K, L, M, N, O, P, Q, R, S, T, U, W, and X) and two fibrous corona-located proteins (CENPE and CENPF) [6]. Malfunctions of kinetochore can cause aneuploidy [7], a well-known hallmark of human cancers. This has prompted the reasonable hypothesis that CENPs play critical roles in tumorigenesis. Indeed, many researchers have confirmed that the dysregulation of CENPs is significantly associated with cancer prognosis and may serve as a biomarker for non-smallcell lung cancer and breast cancer [8,9]. Only two specific CENPs, CENPE and CENPF, have been reported to be possible prognostic biomarkers for ESCC [10,11]. The other CENPs need to be further investigated. The focus of this paper is to find a more accurate model for ESCC than previously reported, such as ZNF750, TTL, AC016205.1, and BOLA3 [12,13].
The Cancer Gene Atlas (TCGA) and Gene Expression Omnibus (GEO) databases have provided a lot of information about gene expression profiles and clinical data related to cancer patients. Bioinformatics analysis provides a number of strategies for cancer prevention and treatment [14]. Herein, we conducted a comprehensive analysis of the expression of all genes encoding CENPs to assess their clinical significance in ESCC using TCGA and GEO data. Differentially expressed genes (DEGs) were identified by comparing the expression differences of the CENPencoding genes between the tumor and normal samples. The effects of the expression of overlapping DEGs from the three datasets on the survival of ESCC patients were determined in a survival analysis. Based on the expression of overlapping DEGs, forecast models were established to forecast the survival of patients with ESCC. Furthermore, we explored the underlying mechanisms using gene set enrichment analysis (GSEA) and coexpression network analysis (WGCNA). Our results suggest that CENPs are promising diagnostic and prognostic biomarkers for ESCC.

Materials and Methods
2.1. Patient Profiles. Profiles of ESCC patients were downloaded from TCGA (https://portal.gdc.cancer.gov/) and the GSE38129 and GSE20347 datasets from GEO (https://www .ncbi.nlm.nih.gov/geo/). In TCGA, the clinical information data of ESCC and EAC were downloaded together. The samples with incomplete information were removed, leaving 95 ESCC patient cases. In addition, there were 30 cases in GSE38129 and 17 cases in the GSE20347 dataset. Relationships between the clinicopathological characteristics and the OS of ESCC patients were determined using univariate analysis using SPSS software (version 23.0).

Analysis of Expression
Difference. The gene expression profiles of CENPs were extracted from the three datasets. In TCGA, some of the 95 clinical samples lacked gene expression, and 81 tumor samples and 11 normal samples were left after selection. In GSE38129, there were 30 normal and 30 tumor samples. In GSE38129, there were 17 normal and 17 tumor samples. Differences in CENP expression between normal and tumor samples were analyzed by calcu-lating log 2 fold change (logFC), false discovery rate (FDR), and P values using edgeR and limma packages in R 3.5.1 software for TCGA and GEO data, respectively [15,16]. Genes with FDR < 0:05, and P < 0:05 were identified as DEGs. Heatmaps, Venn diagrams, and boxplots were plotted using R.
2.3. Survival Analysis. According to the gene expression level, DEGs in tumor samples were divided into low and high expression groups. Gene expression levels were ranked from high to low levels, with the top 50% as the high expression group and the bottom 50% as the low expression group. OS curves were plotted using R software based on the Kaplan-Meier method [17]. Statistical significance was set at P < 0:05.
2.4. Establishment of the Forecast Model. The risk scores of each patient were calculated as the sum of the expression levels of each gene multiplied by its corresponding coefficient using multivariate Cox regression analysis in R software [18]. Based on the risk scores and survival analysis data, time-dependent receiver operating characteristic (ROC) curves were plotted using the "survivalROC," "time-ROC," and "bootstrap" package of R software. Package "survival" was used for multivariate risk regression analysis through Cox proportional hazards model. Sensitivity was the ordinate for true positive rate, and 1-specific was the abscissa for false-positive rate. A nomogram for individual forecast was generated based on the risk score of the multigene model and clinical risk factors using R software [19].
2.5. Analysis of the Mechanism. In GSEA, expression profiles of tumor samples were divided into CENPE-low and CENPE-high groups, as defined by the median expression value of CENPE. Hallmark gene sets, Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets, and oncogenic signature gene sets were used as references. Enriched gene sets were identified using GSEA-3.0.jar (http://software .broadinstitute.org/gsea/downloads.jsp). Gene sets with FDR < 0:25 and P < 0:05 were considered statistically significant [20]. WGCNA was conducted to identify genes coexpressed with CENPs [21]. The visualized network was plotted using Cytoscape 3.6.1. Correlation analysis was performed using Pearson's correlation analysis with R software.
2.6. Validation of Cell Lines. The cell lines, including KYSE 30, 410, 450, 510, 520, and HDF (human dermal fibroblasts), were purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The cells were cultured in RPMI 1640 medium (Gibco Life Sciences, USA) supplemented with 10% fetal bovine serum (Gibco Life Sciences), 100 U/mL penicillin (Gibco; Thermo Fisher Scientific, Inc., USA), and 100 μg/mL streptomycin (Gibco; Thermo Fisher Scientific, Inc.) and were incubated at 37°C in a humidified incubator containing 5% CO 2 . In quantitative real-time PCR analysis, total RNA was extracted using the RNeasy®Mini Kit (QIAGEN, USA). RT-qPCR was performed using SYBR Premix Ex Taq II (TaKaRa BIO, Japan) in a LightCycler® Real-Time PCR System (Roche, Switzerland). The thermocycling conditions included reverse transcription at 50°C for 10 min and initial denaturation at 95°C for 3 min, followed by 40 cycles of denaturation at 95°C for 15 s, annealing, and extension at 60°C for 30 s. The data were calculated using the 2−ΔΔCt method. The forward and reverse primers for CENPE were 5 ′ -CAGC AGAGAAGAATCACTTG-3 ′ and 5 ′ -GTACCATTGTA GCCTTGTATG-3 ′ and for CENPQ were 5 ′ -CAATAC CATCTCAACTTCCTG-3 ′ and 5 ′ -TGTAGTAATGCCAG ACCTTC-3′. Histograms were drawn using the GraphPad Prism 8 software. Statistical significance was set at P < 0:05. The gene expression profile of CENPE in esophageal cell lines was extracted from the GSE23964 dataset (two normal esophageal epithelial cell lines and 14 ESCC cell lines). Gene expression difference analysis was performed using the limma package in R software.

Clinical Characteristics of ESCC Patients in TCGA.
This study explored the effect of CENPs on ESCC ( Figure S1). The relationship between the clinical characteristics and OS of ESCC patients in TCGA was clarified by performing univariate Cox regression analysis. Due to the lack of survival data, the GSE38129 and GSE20347 datasets were only used in the expression difference analysis. Male sex, advanced tumor-node-metastasis (TNM) stage, and N2 and N3 stages were significantly associated with poor survival of ESCC patients (P = 0:020, P = 0:015, and P = 0:012, respectively; Table 1). Therefore, N stage, sex, and TNM stage may be potential risk factors for OS in patients with ESCC.  Figure S2. The collective findings revealed that a variety of CENPs might be promising diagnostic biomarkers for ESCC.

Correlation between Expression of CENPs and Survival of ESCC Patients.
The Kaplan-Meier OS curves were plotted to determine the prognostic value of the overlap of DEGs in ESCC patients. As shown in Figure 2(a), high expression levels of CENPE and CENPQ were significantly correlated with better outcomes in ESCC patients (P = 0:015 and P = 0:038, respectively). Although the other commonly upregulated CENP-encoding genes in three datasets (CENPA, CENPF, CENPI, CENPM, CENPN, and CENPR) did not display statistically significant survival differences between the lower-and higher-expressed groups (P > 0:05; Figure 2 (a)), they exhibited a similar trend with CENPE and CENPQ, where high expression was associated with better survival. Therefore, CENPE and CENPQ may serve as potential prognostic biomarkers for patients with ESCC.
Since the N stage, sex, and TNM stage exhibited notable relationships with OS of ESCC patients in the univariate Cox regression analysis (Table 1), The Kaplan-Meier OS curves were also drawn for the three clinical characteristics. The results further demonstrated that N stage, sex, and TNM

Prognostic Accuracies of CENPE, CENPQ, and the Other
CENPs in ESCC Patients. ROC analysis is a widely applied method to evaluate the prognostic performance of patients using the area under the curve (AUC) as an index [22]. The forecast model is significant only when its AUC value exceeds 0.60 [23][24][25]. By selecting those universally upregulated CENP-encoding genes, we established more focused forecast models to forecast OS. A single gene CENPEbased forecast model was more accurate than the TNM stag- The AUC value (0.8550 of such a multigene forecast model was satisfactory (Figure 3(c)) and was much higher than that of TNM staging, implying that the forecast model has high specificity and sensitivity for ESCC survival forecast. According to the median risk score, patients were divided into low-risk and high-risk groups. The corresponding survival curve demonstrated that low-risk patients had a higher survival rate than the high-risk group (P = 0:014; Figure 3(c)).
A nomogram is a reliable tool for the prognosis of cancer patients by incorporating and illustrating important factors for oncologic prognoses [26]. Based on the aforementioned results of the strong association of TNM stage and sex were tightly associated with patients with ESCC OS (Table 1 and Figure 2(b)), we further constructed a nomogram integrating CENP-based risk score and the two clinicopathological risk factors (TNM stage and sex). The N stage was excluded from the nomogram because it was included in the TNM stage. As shown in the nomogram, the CENP-based risk score contributed the most to forecast patients' OS, followed by sex and TNM stage (Figure 3(d)). The ROC curve showed that the true positive rate of our nomogram integrating CENPbased risk score, TNM stage, and sex could reach 90.6% (Figure 3(e)), implying the extremely high accuracy of the nomogram in forecasting individual OS of ESCC patients. In addition, to clarify whether our forecast models were superior in forecasting the survival of ESCC patients, compared with previously published biomarkers, we also performed ROC analysis of other known biomarkers, including ZNF750, TTL, AC016205.1, and BOLA3 [12,13]. The results demonstrated that our CENPE-based forecast model, CENP-based forecast model, and integrated nomogram all had higher AUC values than the four other known biomarkers (AUC values for TTL, ZNF750, AC016205.1, and BOLA3 of 0.652, 0.643, 0.623, and 0.613, respectively; Figures 3(e) and 3(f)). Collectively, our forecast models based on CENPE, CENPs, and integrating CENP-based risk score, TNM stage, and sex are promising in the prognosis of ESCC patients.

Identification of the Potential Mechanism of CENPE in ESCC Progression.
To investigate the underlying mechanism of CENPE in ESCC progression, we performed GSEA and WGCNA. GSEA is a computational method to explore whether a specific gene set is markedly enriched in a group of gene markers ranked by their relationship to a phenotype of interest [20,27]. In the experiment, the expression profiles of tumor samples were divided into CENPE-low and CENPE-high groups and then analyzed based on hallmark gene sets, KEGG gene sets, and oncogenic signature gene sets. Several cancer-related gene sets, including G2/M checkpoint, mitotic spindle, cell cycle, E2F targets, VEGF, RB/ p107, EGFR, ERB2, and p53, were significantly enriched in the high CENPE expression group (FDR < 0:25 and P < 0:05; Figure 4). To further explore the molecular  mechanism, we performed WGCNA and correlation analyses. In WGCNA, genes coexpressed with CENPs with correlation coefficients > 0:5 were selected and demonstrated in the visualized network ( Figure 5). Some of the coexpressed genes that were highly related to CENPE and/or tumorigeneses, such as TOP2A, NDC80, BRCA1, CENPF, BARD1, TTK, BRCA2, and BUB1B, were further selected to plot correlation maps ( Figure S3(a)). The results of differential expression analysis showed that BRCA1, BUB1B, and TTK were significantly upregulated in ESCC tissues based on the TCGA, GSE38129, and GSE20347 datasets, at the same time notable overexpressed in ESCC cell lines based on the GSE23964 dataset (FDR < 0:05 and P < 0:05; Figures S3(b)-S3(e)).

Target Validation in Human Cell
Lines. To validate the mRNA expression differences of CENPE and CENPQ at the cell line level, RT-qPCR was done, and expression differ-ence analysis was performed using expression profiles extracted from the GSE23964 dataset. Analysis results based on the GSE23964 dataset showed that CENPE was overexpressed by microarray assay in ESCC cell lines (logFC = 1:86, P < 0:001; Figures 6(a) and 6(b)). There were no CENPQ expression data in the GSE23964 dataset. In RT-qPCR analysis, since normal esophageal epithelial cell lines were very difficult to obtain, human dermal fibroblasts (HDF) were chosen as the ESCC control cell line with reference to other published articles [28,29]. RT-qPCR analysis results demonstrated that both CENPE and CENPQ were significantly upregulated in the five ESCC cell lines compared to the normal cell line (P < 0:05; Figures 6(c) and 6 (d)). Overall, these results validated that CENPE and CENPQ were upregulated in ESCC cell lines compared to the normal ones, consistent with their expression difference at the tissue level in TCGA and GEO databases.

Discussion
Tumor markers of esophageal cancer may have pivotal roles in evaluating tumor response to therapy [30], which could be exploited to develop early diagnostic biomarkers. Aberrant expression of CENPs has reportedly been related to several human cancers. For instance, in non-small-cell lung cancer, CENPU expression promotes cancer cell proliferation and forecasts poor survival [8]. In breast cancer, elevated expression of CENPA is associated with cancer malignant progression and is a prognostic biomarker [9]. However, the potential role of CENPs in ESCC has seldom been investigated. Herein, we systematically clarified the potential clinical functions of CENPs in ESCC patients using bioinformatics methods based on multiple datasets.
Most CENP-encoding genes, including CENPA, CENPE, CENPF, CENPI, CENPM, CENPN, CENPQ, and CENPR, were upregulated in ESCC patients in the TCGA, GSE38129, and GSE20347 datasets. Consistently, a previous study demonstrated the overexpression of CENPF in ESCC cell lines at both the mRNA and protein levels when compared to normal tissue [11]. CENPE was also upregulated in ESCC based on the TCGA dataset in a recent study [10]. However, the previous study only analyzed a single gene using a single dataset to draw the conclusion. In contrast, we obtained our results for all CENP-encoding genes based on three datasets, making it more convincing. Except for CENPE and CENPF, to the best of our knowledge, the overexpression of the other six CENPs in ESCC is described for the first time. Additionally, CENPH was overexpressed in ESCC samples compared to normal samples based on the TCGA dataset. Interestingly, CENPH is overexpressed and is prognostic in esophageal carcinoma [31].
In survival analysis, high expression of CENPE and CENPQ was significantly associated with better outcomes in ESCC patients. Similarly, some genes with high expression in pathological tissues will act against certain aberrations of pathological cells. Previous research has shown that the downregulation of CENPE causes an increase in aneuploidy, which in turn triggers an elevated level of spontaneous lymphomas and lung tumors in aged animals [32], implying that CENPE acts as a tumor suppressor. In a recent study, high expression of CENPE was closely correlated with better survival in ESCC patients but with unfavorable outcomes in EAC patients [10]. These findings indicate that CENPE might play crucial and complicated roles in the survival of cancer patients. Importantly, in the current study, CENPQ was first reported to act as a prognostic biomarker for ESCC patients. We speculate that although CENPE is highly expressed in tumor tissues, it has a good effect on survival, which may be a protective factor for ESCC. Thus, CENPE and CENPQ could serve as potential prognostic biomarkers for ESCC patients.
Since there are various factors that affect gene expression, a single gene is usually difficult to be an ideal factor to forecast. Indeed, the single gene CENPE-based model was capable of forecasting the OS of ESCC patients, while CENPQ was not. Therefore, a multigene forecast model based on the expression of CENPs was established. Satisfactorily, the multigene model exhibited an especially higher AUC value than that of TNM staging classified as stage I/II vs. III/IV, with an accuracy of 85.5% (AUC = 0:855). The stage I/II vs. III/IV based TNM staging system is a recognized benchmark for classifying the degree of spread of cancer and is a principal prognostic factor in forecasting the consequences of patients with cancer [33]. In this study, TNM staging was combined with the CENP-based risk score to construct a new nomogram-based forecast model. Surprisingly, the AUC value of the integrated nomogram reached 0.906, implying the high accuracy of our nomogram in the estimation of the individual OS of ESCC patients.
Interestingly, the CENP-based risk score was the most important factor for OS in the nomogram to forecast, suggesting that CENPs are important in the prognosis of ESCC patients.
GSEA results demonstrated that cell cycle, G2/M checkpoint, mitotic spindle, RB/p107, p53, E2F targets, VEGF, ERB2, and EGFR were significantly related to CENPE expression. CENPE is a kinesin-like microtubule motor protein that accumulates maximally in the G2 phase [34]. It plays a crucial role in the cell cycle by forming a link between the attachment of spindle microtubules to kinetochores and the mitotic checkpoint [35]. High gene expression of CENPE is positively correlated with the tumor suppressor pathway. The G2/M

11
BioMed Research International checkpoint pathway can prevent the cell from entering mitosis (M phase). Hence, high CENPE expression inhibits cell division, which leads to a better prognosis in cancer patients. CENPE is regulated by E2F transcription factor 4. This regulation is important in maintaining G2-arrest and is regulated by p130/p107/Rb signaling [36]. Moreover, we found that CENPE was positively related to genes including TOP2A, NDC80, BRCA1, CENPF, BARD1, TTK, BRCA2, and BUB1B/BUBR1 by correlation analysis. CENPE, CENPF, TTK, and BUB1B are all mitotic spindle assembly checkpoint-related genes. The depletion of CENPE and CENPF has been related to the significant disruption of the cell cycle and paclitaxel resistance in ovarian cancer [37]. In addition, CENPE and TOP2A are upregulated in a number of solid cancers and are involved in mitotic cell cycle nodes in breast cancer [38]. Additionally, CENPE, TOP2A, CENPF, TTK, and NDC80 are highly expressed in the cell cycle of basal-like breast cancer [39]. Therefore, our findings indicate that CENPE affects ESCC progression, possibly by regulating cell cycle-related pathways.
In conclusion, CENPs, especially CENPA, CENPE, CENPF, CENPI, CENPM, CENPN, CENPQ, and CENPR, could serve as promising diagnostic biomarkers for ESCC. CENPE and CENPQ may be potential prognostic biomarkers for patients with ESCC. In addition, the CENPE-based model, CENP-based model, and nomogram integrating CENP-based risk score, TNM stage, and sex are especially promising in forecasting OS of ESCC patients. Mechanistically, CENPE may affect the progression of ESCC by regulating cell cyclerelated pathways by interacting with TOP2A, NDC80, and BRCA1. Further studies are required to confirm their detailed roles by performing cell and animal experiments.

Data Availability
The datasets analysed during the current study are available in The Cancer Gene Atlas (https://portal.gdc.cancer.gov/) and Gene Expression Omnibus (https://www.ncbi.nlm.nih .gov/geo/).  Figure S1: the flow chart of the study. From left to right, the leftmost row is all the original data, the contents in the dotted grid are the methods, then the right is the experimental results, and the rightmost row is the general name of the steps. Figure S2: Pearson's correlations among CENPs. Green color represented negative correlations and red color represented positive correlations. The deeper color indicated the higher correlations. Figure S3: (a) Pearson's correlations between CENPE and its coexpressed genes, including TOP2A, NDC80, BRCA1, CENPF, BARD1, TTK, BRCA2, and BUB1B. (b-d) Expression differences of the above 8 genes between the tumor and normal samples in TCGA, GSE38129, and GSE20347 datasets. * * * P < 0:001 and NS P > 0:05. (e) Histograms showing expression differences of BRCA1, BUB1B, and TTK between ESCC cell lines and normal cells based on the GSE23964 dataset (2 normal esophageal normal epithelium normal cell lines and 14 ESCC ones). Normal cell lines and ESCC ones were filled in black and grey colors, respectively. (Supplementary  Materials)