s2fieldsofstudy: Biology Medicine

text:
Challenges Toward the Identification of Predictive Markers for Human Mesenchymal Stromal Cells Chondrogenic Potential

Abstract Human bone marrow derived mesenchymal stromal cells (BMSCs) represent a putative cell source candidate for tissue engineering‐based strategies to repair cartilage and bone. However, traditional isolation of BMSCs by adhesion to plastic leads to very heterogeneous cell populations, accounting for high variability of chondrogenic differentiation outcome, both across donors and across clonally derived strains. Identification of putative surface markers able to select BMSC subpopulations with higher chondrogenic capacity (CC) and reduced variance in chondrogenic differentiation could aid the development of BMSC‐based cartilage and bone regeneration approaches. With the goal to identify predictive markers for chondrogenic BMSC populations, we assessed the gene expression profile of single cell‐derived clones with high and low CC. While a clustering between high and low CC clones was observed for one donor, donor‐to‐donor variability hampered the possibility to achieve conclusive results when different donors were considered. Nevertheless, increased NCAM1/CD56 expression correlated in clones derived from one donor with higher CC, the same trend was observed for three additional donors (even if no significance was achieved). Enriching multiclonal BMSCs for CD56+ expression led to an increase in CC, though still highly affected by donor‐to‐donor variability. Our study finally suggests that definition of predictive marker(s) for BMSCs chondrogenesis is challenged by the large donor heterogeneity of these cells, and by the high complexity and plasticity of the BMSCs system. Multiple pathways and external parameters may be indeed involved in determining the chondrogenic potential of BMSCs, making the identification of putative markers still an open issue. stem cells translational medicine 2019;8:194&11

INTRODUCTION
Human bone marrow derived mesenchymal stromal cells (BMSCs) have been proposed as a promising cell source candidate for cartilage regeneration and bone repair through endochondral ossification [1-3], mainly due to their ease of isolation and abundance in human bone marrow [4,5]. However, their usage in potential). In order to define a more homogenous starting population, several studies have recently investigated different surface markers potentially able to select within the unfractionated bone marrow BMSCs subpopulations with enhanced CFU-F capacity, trilineage differentiation potential or selfrenewal properties. With this goal, several single surface markers for example, CD271 [7], CD49a [8], and CD146 [9] and surface marker combinations including CD271 + CD146 + [10], CD271 + MSCA-1 + CD56 + [11], CD271 + CD90 + CD106 + [12], and CD271 + CD140a − [13,14] were proposed to demarcate these cells. Yet, no consensus has been achieved within the scientific community on such putative marker(s). CD271 + cells isolated from fresh bone marrow were shown to have a better chondrogenic differentiation potential in vitro and at an orthotopic site compared with the cells isolated by plastic adherence, however this was demonstrated only with one donor [15]. In contrast, against some previous findings suggesting a higher chondrogenic capacity (CC) of CD105 + cells compared with bulk or CD105 − cells, Cleary et al. [16] demonstrated that there is no correlation between the CC and the CD105 expression level at different passages of multiclonal BMSCs. It is thus not clear yet whether any marker can delineate cell subsets with higher CC within the bulk population. Moreover, single BMSC derived clones showed a broad range of chondrogenic and osteogenic in vitro differentiation capacity, confirming the heterogeneity of the bulk population [2,17]. In a recent report, single cell RNA transcript levels of chondrogenic genes in bovine BMSCs did not correlate with their functional properties thus challenging a preselection of naïve BMSCs based on specific markers in order to improve chondrogenic outcome [18].
The goal of this study was thus to identify markers for subpopulations of BMSCs exhibiting a high and reproducible chondrogenic differentiation potential across different donors. To this end, we hypothesized the existence of a correlation between the transcriptomic profile of single cell-derived clones and their chondrogenic differentiation capacity. Based on differentially expressed genes identified by transcriptomic profiling, we then found CD56 to be upregulated in clones with high CC. Subsequently, we sorted multiclonal BMSCs based on their CD56 protein expression in order to verify, whether the CD56 + cell subset has improved differentiation potential compared with the negative fraction and the bulk population. Even if enriching expanded BMSCs for CD56 + expression led to an increase in CC, we showed that donor-to-donor variability drastically challenged the possibility to achieve standardized solutions to enhance BMSC CC.

MATERIALS AND METHODS
If not otherwise stated, all reagents were purchased from Sigma-Aldrich (St. Louis, MO). Cell culture media and supplements were purchased from Gibco (Grand Island, NY).

Cell Isolation and Expansion
Bone marrow aspirates and cartilage tissue biopsies were obtained from patients after informed consent during surgical procedures in accordance with the local ethical committee (University Hospital Basel; Prof. Dr. Kummer; approval date 26/03/2007 Ref Number 78/07). Human mesenchymal stromal cells were isolated from bone marrow aspirates (n = 11; mean age: 39 AE 13 years, 9:2 male : female) by plating 0.10-0.13 × 10 6 of nucleated cells/cm 2 in alpha-mimimum essential medium (MEM) with 10% fetal bovine serum (FBS), 100 mM HEPES buffer, 1 mM sodium pyruvate, 100 IU/ml penicillin, 100 μg/ml streptomycin and 0.29 mg/ml glutamate (complete medium, CM), supplemented with 5 ng/ml FGF2 (BioTechne, Minneapolis, MN, USA) and expanded in the same medium for one additional passages before using them for chondrogenic assays or sorting experiments.

Generation of Single-Cell Derived Clones
The experimental set-up for single cell-derived clones generation is depicted in Figure 1A. In details, 10,000 nucleated cells from 4 fresh bone marrow aspirates (donor 1: male, 23 years; donor 2: male, 38 years; donor 3: female, 49 years; donor 4: female, 30 years) were seeded in 96 well plates, being the expected frequency of clonogenic cells of approximately 1/10 4 . After 1 week, wells were checked for colony growth: wells without cells and wells with two or more colonies were discarded. Upon confluence, colonies were passaged first into three 12-well plate wells and subsequently, according to the cell number, into culture flasks with a cell density of 500-2,000 cells/cm 2 . When a minimum of 0.8 × 10 6 cells was reached, cells were subjected to RNA isolation for transcriptomic analysis and chondrogenic assays. Senescent and nongrowing clones were discarded.
For the clonal study of CD56 + and CD56 − expressing clones, P1-expanded cells from one donor (age: 33 years, male) were single-cell sorted into 6 × 96 well plates based on CD56 expression. Wells with growing colonies derived from one cell were expanded subsequently as described above. From the (in total 82) generated CD56 + clones, clones that did not reach confluence in 12 well plate after 20 days of expansion were discarded. After expansion, clones were tested for their capacity to differentiate toward the chondrogenic, osteogenic and adipogenic lineages, and for their expression of CD56 at mRNA level.

Adipogenic Differentiation
Cells were seeded into 12 well plates (0.04 × 10 6 cells/well, 1-3 replicates depending on the gained cell number) and cultured for 2 weeks in adipogenic differentiation medium consisting in CM based on DMEM, 4.5 g/l D-glucose supplemented with 10 μg/ml methyl-iso buthyl xanthine, 10 μg/ml insulin Actrapid HM, 0.1 mM indomethacin and 10 −6 M dexamethasone and for one additional week in CM containing insulin. After the culture, cells were fixed for 10 minutes in formalin and stained with oil red. For semiquantitative analysis, three

Osteogenic Differentiation
Cells were seeded into 12 well plates (0.04 × 10 6 cells/well), 1-3 replicates depending on the gained cell number) and cultured for 3 weeks in osteogenic differentiation medium consisting in CM based on alpha-MEM supplemented with 0.1 mM ascorbic acids 2-phosphate, 10 mM beta-glycero-phosphate and 10 −8 M dexamethasone. Afterward, the cells were fixed for 10 minutes in formalin and stained with alizarin red. For quantitation purposes, after image acquisition the alizarin red was solubilized by incubation in 10% acetic acid, subsequently at 85 C and with 10% ammonium hydroxide neutralized supernatants were optically measured at a wavelength of 405 nm.

Cytofluorometry and Fluorescently Activated Cell Sorting
BMSCs at different passages were analyzed using the fluorophore-conjugated antibodies against the following sur-  3,300 cells/cm 2 and expanded until fully confluent or not more than 12 days before testing in chondrogenic assays.

Histology
Pellets after in vitro culture were fixed in 4% paraformaldehyde, dehydrated, and embedded in paraffin. Safranin-O/fast green (SafO/FG) staining with hematoxylin (J.T. Baker) nuclear counterstaining was performed to analyze cartilage tissue formation. For semiquantitative assessment of the chondrogenesis the Bern score according to Grogan et al. [19] was determined.

Quantification of Glycosaminoglycans and DNA
Pellets were digested with 1 mg/ml protease K in 50 mM Tris with 1 mM EDTA, 1 mM iodoacetamide and 10 mg/ml pepstatin-A for 16 hours at 56 C. For glycosaminoglycan quantification, the method of Barbosa et al. [20] was used. Briefly, diluted or undiluted digested constructs (depending on SafO intensity) were incubated with 1 ml of dimethylmethylene blue assay (DMMB) solution (16 mg/l dimethylmethylene blue, 6 mM sodium formate, 200 mM GuHCL, pH 3.0) on a shaker at room temperature for 30 minutes. Precipitated DMMB-glycosaminoglycans (GAGs) complexes were centrifuged and supernatants were discarded. Complexes were dissolved in decomplexion solution (4 M GuHCL, 50 mM Na-Acetate, 10% Propan-1-ol, pH 6.8) at 60 C, absorption was measured at 656 nm and GAG concentrations were calculated using a standard curve prepared with purified bovine chondroitin sulfate. DNA content was measured by using the CyQuant Cell Proliferation Assay Kit (Molecular Probes Inc., Eugene, OR) according to the instructions of the manufacturer.

RNA Sequencing
Based on SafO/FG staining, the clones were evaluated and divided into two groups: clones with a Bern score ≥3 were classified to have a high CC and clones with a Bern score <3 as clones of low CC. Batch I analysis included 14 clones of donor 1 (6 of high CC), while batch II comprised 6 clones of donor 2 (3 of high CC), 8 of donor 3 (4 of high CC) and 6 clones of donor 4 (3 of high CC). For mRNA isolation, the protocol of the column-based quick-RNA min-prep (ZymoResearche, Orange, CA, USA) was followed. mRNA was quantified by using nanodrop and its quality was determined with the help of the 2,100 bioanalyzer (Agilent Technologies, Palo Alto, CA). Three hundred nanograms of DNase-digested mRNA was used for library preparation with TruSeq Stranded Total RNA Kit with Ribo-Zero Gold (Illumina, San Diego, CA, USA). Libraries were run on the fragment analyzer (Advanced Analytical, Heidelberg, Germania) for quality control and were adjusted to equal concentrations. Sequencing was performed with a HighSeq 2,500, with single end reads of 50 base lengths. Reads were aligned to the human genome (UCSC hg19) using the spliced read aligner STAR [21]. Sequencing and mapping quality was assessed by the qQCReport function of the R package QuasR [22]. Gene expression was quantified by the qCount function of QuasR using an exon-union model of hg19 RefSeq genes (downloaded from UCSC 2015-02-09). The R package edgeR [23] was used for differential gene expression analysis. For batch I, two different approaches were tested: first only a subset of clones was compared based on similar proliferation rates, then all clones were analyzed in an unbiased way without accounting for the proliferation rate. In both cases, an additive generalized linear model including the factor "CC" (high and low) was fitted to the raw counts (function glmFit) and differential expression between samples with high and low CC was evaluated by likelihood ratio tests (function glmLRT).
p-Values were adjusted by controlling the false-discovery rate.
All genes with an adjusted p-value <.05 were considered to be significantly differentially expressed. In case of batch II, the additive model included the additional factor "donor" (BM2, BM3, and BM4). All calculations were performed at sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel. Significantly differentially expressed genes of the subset analysis of batch I were additionally compared with the panel of surface markers (QIAGEN, Hilden, Germany) to identify differentially expressed surface markers and subjected to ingenuity pathway analysis (QUIAGEN, Hilden, Germany) to identify differently activated signaling networks.

Quantitative Real-Time Polymerase chain reaction (qPCR)
The mRNA was isolated as described in the RNAsequencing paragraph except that the DNA digestion was not performed. cDNA was generated with random primer (Promega, Madison, WI) and the kit of SuperScript reverse transcriptase III (Thermo Scientific, Waltham, MA, USA). For Quantitative Real-Time PCR (qPCR) the TaqMan assay on demand system based on FAM and Tamra (Thermo Scientific, Waltham, MA, USA) was used. Samples were analyzed for expression of CD56 (NCAM1; Hs00941830) and as reference gene GAPDH (Hs02758991) was used.

Statistical Analysis
Statistical analysis was performed with Graph Pad Prism 6 software. For comparing one variable from different donors or single cell derived clones, unpaired t tests were applied. Pearson d'Agostino test was done to check for normality of the data sets. For multiple comparison tests of non-normally distributed data sets, one-way analysis of variance (ANOVA) was performed using Kruskal-Wallis test with Dunn's multiple comparison test as post hoc. Statistical significance was considered with a p-value <.05. All the data are shown as average AE standard deviation. n indicates number of donors tested and analyzed.

Isolation of Single Cell Derived Clones and Assessment of their CC
In order to find a molecular signature of BMSCs with higher CC, we isolated and expanded single cell derived clones (n = 4 donors) and assessed their gene expression and their chondrogenic in vitro differentiation potential (Fig. 1A). Overall, in 23.2% AE 11.5% of the seeded wells a single cell derived colony could be identified. Clones were expanded for 28 AE 2 days, with a mean proliferation rate of 0.7 AE 0.1 doublings/day for the different donors (Fig. 1B, 1E). As expected, chondrogenic differentiation occurred to a different extent among the clones (Fig. 1C, 1D, 1F) Fig. 2A). In order to reduce the confounding effect of the different proliferation rates, in a second analysis only clones with a similar proliferation rate (0.67-0.70) were taken into account (9 out of 14), while outliers in terms of proliferation rate were discarded. Interestingly, the clustering became more prominent and the number of significantly differentially expressed genes rose to 1,771 (Fig. 2B). This demonstrated that clones with a closer proliferation rate shared more common gene expression patterns and that the different doubling times introduced more heterogeneity, eventually masking the putative similarities among clones with equal CC. Among the significantly upregulated genes (log2-fold-change >0.5 in clones with high CC), there were extracellular matrix proteins which are important components of the cartilage ECM, such as hyaluronan and proteoglycan link protein (HALPN1), aggrecan (ACAN) and cartilage oligomeric matrix protein (COMP; Fig. 2C). Moreover, the gene expression of several surface markers (NCAM1/CD56, VACAM1/CD106, ALCAM/CD166, ITGA1/CD49a) was significantly higher in the clones with high CC, while the gene expression of two surface markers (ENG/CD105 and CSF1R/CD115) was significantly lower expressed in the clones with high CC (Fig. 2C). Notably, NCAM1 was the mostly upregulated surface marker by the clones with high CC (3.25 log2-fold-change). The significant upregulation of NCAM1 in clones with high CC was also confirmed by qPCR (Fig. 2D). Ingenuity pathway analysis (IPA) showed that downstream targets of PDGFBB signaling were downregulated and that genes implicated in chondrogenesis were rather increased in clones with high CC (Fig. 2C). According to IPA canonical pathway prediction, clones with high CC exhibited a potential inhibitory status of EIF2 pathway, while the estrogen, mTOR, eIF4, and IL8 signaling pathways were potentially activated (Fig. 2E). RNA sequencing was then performed for clones combined from the three additional different donors (donors 2, 3, and 4; batch II), in order to determine if the genes identified in the first batch (donor 1) are robustly related to the BMSCs CC. Principle component analysis showed that clones preferentially grouped according to their donor origin, rather than to their CC. Moreover, clones from batch II clustered furthest apart from the clones of batch I (Fig. 3A). Since the number of samples per donor was relatively low (a lower number of clones were obtained for donor 2, 3, and 4 as compared with donor 1), the differential gene expression analysis was done for all the clones irrespective of their proliferation rate. Figure 3B shows the genes that were significantly differentially expressed between samples of high and low CC in batch II analysis. None of these six genes overlapped with the genes found in batch I (independent of whether all clones or the subset with similar proliferation rate were taken into account). Notably, a trend of increased expression of CD56 was apparent in the clones with high CC both by mRNA sequencing (Fig. 3C) and qPCR (Fig. 3D). However, the range of CD56 expression levels was too broad to lead to statistical significance.

Expression of Surface Markers Identified by the Transcriptomic Analysis in Expanded Multiclonal BMSCs
Aiming at finding surface marker(s) that select for BMSCs with higher CC within the bulk population, expanded multiclonal BMSCs (P1) were analyzed for the surface markers found to be upregulated in clones with high CC in the subset analysis of batch I. On average 87% AE 3% of the cells were positive for CD49a (n = 3 donors), 93.9% AE 7.1% for CD106 (n = 7 donors), while CD105 and CD166 were expressed on all of the cells (Fig. 4A). The expression of ROR2, previously correlated with chondrogenic potential of BMSC [24], was also checked. Notably, in our study BMSCs did not express ROR2 (percentage of ROR2 + cells lower than 0.1%, data not shown). The frequency of CD56 + cells instead ranged between 1% and 35%, with an average of 13% AE 11% CD56 + cells (n = 12 donors; Fig. 4B). The frequency of CD56 + cells in bulk multiclonal BMSCs slightly decreased from 19% AE 13% at P1 to 14 AE 18 at P3, while the frequency of CD56 + in the positively sorted cell subset decreased from 94% AE 2% directly after sorting to 33% AE 9% after two passages. The negatively sorted cell subset showed an increase of CD56 + cells from 0.8 AE 0.4 at P1 to 4.5 AE 2.6 at P3 (Fig. 4C). Previously, the frequency of BMSCs positive for CD56 in P2-P3 cells was determined between 24% and 89% [25]. Markers exhibiting high and stable expression profiles are unlikely to be unique markers of BMSCs with high chondrogenic potential. Instead, markers with low expression, reflecting the low presence of chondrogenic clones in the total BMSCs population, are more likely to be selective markers of chondrogenic BMSCs [26]. Therefore, we pursued CD56 for comparative studies of CC on bulk BMSC populations.

Chondrogenic Potential of Sorted Multiclonal CD56 +/− BMSCs
To verify the potential of CD56 to select for BMSCs with high CC, P1 expanded BMSCs were sorted and CD56 positive and negative subpopulations were tested in chondrogenic cultures, using unsorted BMSCs as a control (Fig. 5A). The gene expression levels of CD56 directly after sorting and after expansion were more than one order of magnitude higher in the CD56 + cells as compared with the negative ones and corresponded to the level found in the clones of batch I, while the unsorted cells showed an intermediate expression level (Fig. 5B). This demonstrates that the two sorted populations clearly differed in terms of CD56 expression. Importantly, trends in the expression levels among cell populations were maintained after expansion (Fig. 5B). Three out of the 11 BMSC donors were assessed in terms of chondrogenic potential by means of pellet culture (the other BMSC donors had to be discarded due to insufficient cell yield). For all BMSC donors, chondrogenically cultured CD56 + cells formed pellets with higher GAG contents as shown by SafO/FG staining (Fig. 5C). Coherently, the Bern score assessed by histological analysis was higher for the CD56 + cells as compared with either the negative or the bulk populations (Fig. 5D). On average, GAG/DNA content was significantly higher in CD56 + BMSC pellets as compared with CD56 − and bulk ones (Fig. 5E).

STEM CELLS TRANSLATIONAL MEDICINE
However, the GAG/DNA content of single donors (Fig. 5F) showed that the averaged results originated mainly from cells of donor 1. Thus, high donor-to-donor heterogeneity affected the robustness of the results.

Trilineage Differentiation Potential of Single Cell Derived Clones from Sorted CD56 +/− BMSCs
In order to better characterize the heterogeneity within the CD56 sorted cell populations, P1-expanded BMSCs were sorted based on CD56 protein expression and single cells were collected for clonal expansions (Fig. 6A). 10.4% of the CD56 + and 5.9% of the CD56 − derived clones survived and proliferated. Clones were expanded for 24 AE 3 and 23 AE 4 days and had a mean proliferation rate of 0.81 AE 0.11 and 0.86 AE 0.14 doublings/day, respectively (Fig. 6B). The gene expression of CD56 at the end of the clonal expansion phase did not correlate with the clone origin-namely derivation from a CD56 + or a CD56 − single cell-or with their chondrogenic potential (Fig. 6C). After expansion, clones were subjected to in vitro trilineage differentiation assays. Regarding chondrogenic potential, 1 clone out of 11 (i.e., 9.1%) and 4 clones out of 22 (i.e., 18.2%) showed high CC within the CD56 − and CD56 + derived clones, respectively, not resulting in a significant difference between the two subpopulation derived clones (Fig. 6D). The CD56 + clones were marginally osteogenic and adipogenic, with the exception of one clone. All CD56 − clones were instead osteogenic (Fig. 6E) and the majority underwent adipogenesis (Fig. 6F), as measured by alizarin red staining for calcification and oil red staining for the accumulation of fat droplets, respectively. The different potential of CD56 + and CD56 − single cellderived clones in differentiating toward osteogenic and adipogenic linages suggests a functional distinction between these two BMSCs subpopulations.

DISCUSSION
In this study, we explored the existence of predictive markers identifying BMSCs subsets with high CC, starting  from analyzing the molecular signature of clones derived from single BMSCs. Transcriptomic data from one BMSC donor gave a clear-cut difference between clones with high and low CC and hinted NCAM/CD56 as a marker to prospectively isolate a subpopulation of BMSCs with superior chondrogenic potential. However, results also underlined a high BMSCs donor-to-donor variability that challenges the possibility to identify genes regulated in a similar fashion by different donors and to simply apply a marker identified for one donor (i.e., CD56) to others. Aiming to overcome the heterogeneity of bulk BMSCs isolated by plastic adherence, we assessed whether gene expression correlated with the CC of individual single BMSC-derived clones. Differential gene expression analysis resulted in the segregation of clones with high and low CC, indicating a correlation between the gene expression profile and the CC of the underlying clone. Notably, the segregation of the two groups was even more prominent when clones with a different proliferation rate than the average (proliferation outliers) were excluded from the analysis. RNA sequencing was performed on three additional BMSC donors to assess whether the results found by the transcriptomic analysis of one donor could be extended to different donors. Overall, our analysis suggested that intrinsic variability among samples (either clone-to-clone for example, proliferation rate or donor-to-donor) exceeded a potential gene expression difference based on CC, masking the putative similarities among clones with similar chondrogenic potential. It could not be excluded that variable culture conditions such as the duration of expansion or different cellular confluence states at the time of harvest introduced additional variance. The difficulty in extrapolating predictive markers from differential gene expression analysis among single BMSCderived clones was reported previously with a bovine bone marrow donor [18]. Our results highlight the requirement for a more comprehensive analysis, considering a significantly higher number of donors and clones per donor. This stands in contrast to the study of Dickinson et al. [24], in which RNA sequencing of clones derived from a single donor was performed and suggested ROR2 as a putative marker identifying an highly chondrogenic BMSCs subpopulation. Notably, chondrogenic clones derived from the donors considered in our study did not over-express ROR2.
Interestingly, when excluding the proliferation rate as a confounding effect in the analysis of batch I, several surface markers were found to be significantly upregulated in clones with high CC and were thus considered for further analyses. Among them, CD56 was the surface marker with the highest log2-fold-change difference between clones with high and low CC. Moreover, CD56 was expressed at different moderate levels in different donors of BMSCs heterogeneous populations. Together with N-cadherin, CD56 has previously been demonstrated to be an important mediator for cell-cell adhesion during prechondrogenic condensations of chondroprogenitors in limb bud development [27][28][29]. Since the initial condensation phase is also critical for triggering in vitro chondrogenesis of adult BMSCs [27], even if no evidence support, yet the role of CD56 in this process, we hypothesized that an enrichment in CD56 could be beneficial. CD56 + sorted cells from expanded multiclonal BMSCs indeed exhibited a higher chondrogenic differentiation potential compared with the bulk or the CD56 − cells. However, differences among donors affected the robustness of these results and this could be mainly attributed to the low number of tested donors. Eight out of 11 donors had to be excluded from this analysis due to low proliferation rates of corresponding cells after sorting. For this reason, future efforts should also be given in defining optimal expansion conditions for sorted BMSCs.
To further investigate the heterogeneity within the CD56 sorted cell populations, single CD56 +/− cell-derived clones were compared in trilineage differentiation assays. A clear distinction between the two subpopulations was evident in terms of osteogenic and adipogenic differentiation capacity (high in CD56 − clones and limited/absent in CD56 + ones). Differences in chondrogenic potential were less conclusive, with a slightly higher fraction of Safranin-O positive clones in the CD56 + compared with the CD56 − one. Previously, Battula et al. sorted CD56 + BMSCs from fresh CD45-depleted human bone marrow in combination with MSCA-1 and CD271. Interestingly, the authors observed higher CFU-Fs and better in vitro chondrogenic differentiation of CD271 + MSCA-1 + CD56 + compared with the CD56 − cells, along with superior adipogenic differentiation of CD56 − cells at multiclonal and clonal levels [11]. Our results are in line with this previous study, although different sorting times were applied. In the present work, BMSCs were expanded for one passage prior to sort CD56 + and CD56 − populations, thus leading to a higher cell yield. This indicates that CD56 expression identifies a certain cell subpopulation not only in naïve BMSCs [11], but also in expanded one, extending the clinical relevance of CD56 as a prospective surface marker. Indeed, from a translational point of view, a key criterion for the applicability of a selection strategy is the possibility to couple it with the achievement of a clinically relevant cell yield.
Even if attempts of simplification through the selection of a single marker for distinguishing BMSC subpopulations may drastically increase the clinical relevance of this approach, the intrinsic complexity and plasticity of the BMSC system challenges this conclusion. First, variables related to donor-todonor heterogeneity should be taken into consideration. Results obtained from one donor were not immediately transferable to other donors, highlighting the importance of considering larger donor cohorts to achieve robust and generally applicable hints.

CONCLUSION
Moreover, lessons from different stem cell fields (e.g., cancer stem cells, hematopoietic stem cells) suggest that a high level of complexity characterizes stem cells systems. For example, it starts to be evident that tumor progenitor cells [30], as well as hematopoietic stem cells [31], are hardly classifiable through a simple combination of markers, but are most likely characterized by spatio-temporal plasticity. In a similar fashion, we here propose that the BMSC system may also be challenged by the same level of complexity: multiple intrinsically activated pathways and various external parameters (e.g., donor origin, cell derivation, isolation method, expansion protocol, and so forth) may be involved in determining the fate of BMSC subpopulations. This awareness should alert the field approximately the difficulty of identifying univocal preselection criteria following the current approaches, while possibly opening the path to innovative computational methods eventually able to process multiple variables of the BMSC system.


extfieldsofstudy: Biology Medicine

text:
CancerClock: A DNA Methylation Age Predictor to Identify and Characterize Aging Clock in Pan-Cancer

Many biological indicators related to chronological age have been proposed. Recent studies found that epigenetic clock or DNA methylation age is highly correlated with chronological age. In particular, a significant difference between DNA methylation age (m-age) and chronological age was observed in cancers. However, the prediction and characterization of m-age in pan-cancer remains an explored area. In this study, 1,631 age-related methylation sites in normal tissues were discovered and analyzed. A comprehensive computational model named CancerClock was constructed to predict the m-age for normal samples based on methylation levels of the extracted methylation sites. LASSO linear regression model was used to screen and train the CancerClock model in normal tissues. The accuracy of CancerClock has proved to be 81%, and the correlation value between chronological age and m-age was 0.939 (P < 0.01). Next, CancerClock was used to evaluate the difference between m-age and chronological age for 33 cancer types from TCGA. There were significant differences between predicted m-age and chronological age in large number of cancer samples. These cancer samples were defined as “age-related cancer samples” and they have some differential methylation sites. The differences between predicted m-age and chronological age may contribute to cancer development. Some of these differential methylation sites were associated with cancer survival. CancerClock provided assistance in estimating the m-age in normal and cancer samples. The changes between m-age and chronological age may improve the diagnosis and prognosis of cancers.

More and more diagnostic, recurrence and prognostic biomarkers of cancer have been identified (Zhou et al., 2015a(Zhou et al., , 2016. Many age-related biomarkers such as DNA methylation, telomere length changes, cell metabolic factor regulation, transcription level, and protein expression level have been discovered and further studied (Aubert and Lansdorp, 2008;Vilchez et al., 2014;Zubakov et al., 2016). DNA methylation could modulate gene expression during development and cancer progression (Wang et al., 2013;Zhang et al., 2019). In comparison to normal somatic tissues, the cancer methylome is typically characterized by a pattern of global hypomethylation coupled with site-specific promoter hypermethylation (Urbano et al., 2019). Through comparative analysis of a variety of existing age-related biomarkers, DNA methylation has been found as one of the most promising biomarkers for age prediction. DNA methylation refers to the addition of a methyl group to C at the 5th position of cytosine to promote or inhibit gene expression (Dirk, 2015;Nwanaji-Enwerem et al., 2018). In recent years, DNA methylation age (m-age) biomarkers has been able to accurately estimate the age of any tissue throughout life (Horvath, 2013). In addition, m-age biomarkers are valuable tools for evaluating tumor process, which can be a predictor of human health (Perna et al., 2016;Dhingra et al., 2018). However, the differences between m-age and chronological age of different tissues in tumor and the biological processes involved remained to be studied.
There are many studies on the age of individuals using DNA methylation level (Horvath, 2013;Galamb et al., 2016;Dhingra et al., 2018). In 2013, it was proposed for the first time to construct human multi-tissue age predictor based on DNA methylation to measure the degree of aging in human (Horvath, 2013). Subsequently, Weidnei et al. constructed age predictors using three CpG sites as characteristics in 575 healthy samples (Weidner et al., 2014;Nwanaji-Enwerem et al., 2018). Zbiec-piekarska et al. used 5 CpG loci as characteristics in 420 healthy samples to predict age (Spólnicka et al., 2018). Due to tissue specificity of methylation level, Stubbs et al. proposed to construct a multi-tissue age predictor based on methylation level in biological mice model (Stubbs et al., 2017). In addition, DNA methylation is maintained throughout life course in multiple tissues, linking many known early life factor to cancer risk (Kristina et al., 2015). The changes in DNA methylation with age occur at regulatory regions and contribute to tumor development (Johnson et al., 2017). Methylation of SLFN11 is a biomarker for poor prognosis in colorectal cancer and methylations of SLIT1, SLIT2, and SLIT3 are abnormal in gastric cancer. F2RL3 methylation is recently identified as a biomarker closely reflecting both current and past smoking exposure, causing lung cancer (Yan et al., 2015;Kim et al., 2016;He et al., 2017). However, the difference between the m-age for cancer and healthy samples is still unknown. In addition, the effect of age-related methylation characteristics on the survival risk of diverse cancer needs to be further confirmed.
In the present study, we found that DNA methylation levels were correlated with chronological age, and ageassociated methylation sites were identified in human cancers. A comprehensive age prediction model named CancerClock was constructed using these age-associated methylation sites. Further, the characteristics of these age-associated methylation sites were described from three main aspects, including GO functions, biological phenotype, and the genomic location of the feature. Moreover, the m-age of common types of human cancers was described. The differences between predicted m-age and chronological age were also identified using CancerClock. For each cancer type, we detected a series of methylation sites, which could influence m-age compared to healthy samples. Finally, these differences were analyzed using cox regression model, and their effect on tumor survival risk was analyzed. In summary, the findings of this study would provide assistance to depict the age clock of cancers, characterize the m-age of cancer samples, and clinically evaluate the cancer progression.

Clinical and Methylation Profile of Cancers
The methylation profile of IlluminaMethylation450 chip for 33 cancer types and their matched normal tissues were obtained from the cancer genome atlas (TCGA) data portal (TCGA Release 14.0, https://portal.gdc.cancer.gov/) (Hutter and Zenklusen, 2018), and each sample contained 485,512 methylation sites. The values in the methylation profile represent methylation degree of the sites ranging from 0 to 1, which is calculated as the ratio of methylated signal. The sample age ranged from 14 to 89 years. Totally, 8,692 samples including 7,988 tumor samples and 704 normal tissue samples were extracted for the present study (Supplementary Table S1). The corresponding clinical information of the samples (including age, gender, survival prognosis, cancer stage, and type of tumor samples) was also described.

Methylation Data Processing
First, we calculated the mean, standard deviation, and maximum and minimum values of age for each cancer type sample without clinical information were excluded. Methylation sites with NA values >10% in samples were removed. Then, we imputed the remaining NA values using 10-nearest neighbor method with the function knnImputation in DMwR package by R.
DNA methylation profiles of normal and tumor samples were annotated using CpG probe annotation file from GENCODE (Release 29, https://www.gencodegenes.org) (Frankish et al., 2019). All processes were performed by R software (R 3.3.3).

The Segmentation and Statistics of Age for All Cancer and Normal Samples
Human age in this study was divided into four sections according to the novel age subsection proposed by the United Nations and World Health Organization (WHO) after 1994: (Eeghen et al., 2015) Young people: under 44 years old (Taber et al., 2016), middle-aged people: between 45 and 59 years old (Mcguire et al., 2015), older-young people: between 60 and 74 years old, and (Palmer et al., 2018) old people: over 75 years old (Supplementary Table S2). The chronological ages of samples were obtained from clinical information and were divided using above age segmentations rules (Supplementary Table S3).

Identification of Age-Associated Methylation Sites
To identify age-associated methylation sites, Spearman correlation coefficients (SCCs) were calculated based on age and methylation levels based on all the normal samples (Stubbs et al., 2017). Then, a series of methylation sites were obtained with P < 0.05 and |SCC|>0.3. These methylation sites were proved to be significantly correlated with age, and can be the candidate features for constructing the model.

Construction of Age Predictor Known as CancerClock Based on Methylation Level in Normal Samples
To predict human age, we performed LASSO regression model in the generalized linear mode, which through some regression coefficients was strictly set to 0 and we obtained a model with good performance and strong explanatory power (Hepp et al., 2016). Firstly, we obtained a set of M * N methylation matrix X, in which the value in the matrix is the methylation level of the methylation site M in N samples. Each row of X represented the 1,631 age-associated sites, and the column represented the normal samples. Then, we normalized the methylation values as follows: From the formulae, i = 1,2,3. . . m, j = 1,2,3. . . n, m represents the 1,631 age-associated sites, n represents the sample size, and x ij represents the methylation level. Next, we assumed a set of linear regression models: mAge = β 0 + β 1 x 1 + β 2 x 2 + β 3 x 3 . . . + β m x m mAge is the model prediction age, vector β is the model coefficient, and vector x is the methylation expression level.
RSSs (residual sum of squares) with the least square estimation: represents the sample size, y i describes the actual age of sample i, and x ij represents the methylation level. The smaller the RSSs, the better the fitting effect. LASSO-minimality: N means the total number of samples. γ p j = 1 |β j | represents compressed penalty term. γ means adjustment, which controls the influence of RSS and compressed penalty term on regression coefficient. β 0 is excluded from the compressed penalty term, which represents the mean of response variables when x j = 0. Next, coordinate axis descent method was used to solve the model. The initial β values were all 0. When the βvalue was lower than the threshold, we set the two adjacent iterations; the vector β k would be regarded as the final feature of our model.
The CancerClock model was constructed based on GLMNET package in R 3.3.3 (Friedman et al., 2010). We evaluated the stability and accuracy of the model by means of 10-fold crossvalidation. All the healthy samples were randomly divided into 10 subsets, 9 of which were selected as the training set. The test set was used in the verification of the model.

Functional Analysis and Phenotypic Traits for Characteristics in CancerClock
With the Enrichr tool online web server (http://amp.pharm. mssm.edu/Enrichr) using default parameters, functional enrichment was performed for genes across core clusters (Kuleshov et al., 2016), and we obtained enriched GO terms (P < 0.01, FDR < 0.05). The degree of methylation at chromosomal sites affected post-translational modification of proteins, which in turn contributed to the expression of biological traits. In this study, EWAS atlas, which is online analysis software, was used in describing the selected 282 loci with biological traits. EWAS Atlas (http://bigd.big.ac.cn/ewas) ) is a database of epigenetic traits, which include the correlation of information between nearly 330,000 methylation sites and 305 different biological traits through the integration and analysis of existing literatures.

Characteristics of Genome Location Regions for Age-Related Methylation Sties in CancerClock
To explore the distribution of CancerClock features in genome, we got the annotation information of the genome region from GENCODE database (Release 29, https://www.gencodegenes. org) and we mapped the methylation features on the genome using BEDTools (V2.29.0, https://bedtools.readthedocs.io/en/ latest/content/overview.html). In addition, we defined the promoter region as the upstream and downstream 1.5 kb of TSS (transcription start site); thus, the regions included promoter, exon, 5 ′ -UTR, 3 ′ -UTR, and TSS.

Prediction of m-Age for Tumor Types
To predict human methylation age, we used 282 age-related methylation characteristics and 704 healthy samples to construct age prediction model. To investigate further whether tumor causes changes in methylation age, we used the CancerClock prediction model to predict the 33 cancer types in TCGA based on methylation level. For cancer samples, their methylation age and chronological age showed some difference. Among the age prediction of cancer samples, we defined the error value as the absolute age difference of the quarter-point after the age prediction rank, and we showed the difference between the biology age and methylation age for different cancer types. Thus, some age-influenced cancer samples, in which methylation age was significantly different from chronological age, were obtained.

Differential Analysis of Methylation Sites Between Normal and "Age-Related Cancer Samples"
To determine which age-related methylation sites were changed in tumor samples, we conducted t-test for the methylation sites in healthy samples and age-influenced cancer samples, and we obtained significantly differential methylation sites (FDR < 0.01, Fold Change value>2 or <0.5). These cancer samples were defined as "age-related cancer samples" and they have some differential methylation sites.

The Survival Analysis for Differential Methylation Sites
An integrated pipeline was constructed to explore the associations between these age-related differential methylation sites and cancer survival (Zhou et al., 2019). First, we divided the samples into two groups including hyper-methylated and hypomethylated sample groups according to the median methylation level. Next, the age-related differential methylation sites of the 16 cancer types were selected to establish a cox risk regression model for each tumor type, which was selected according to P < 0.05 (Gellar et al., 2015). The covariance matrix was composed of the cancer stage, gender, age, and methylation level. Next, the risk score for each cancer patient was calculated according to the linear combination of the methylation values weighted by the coefficient from multivariate Cox regression analysis. The median risk score was used as the cut-off point to divide the cancer patients into high and low risk groups. Finally, Kaplan-Meier survival analysis was conducted for the two groups, and statistical significance was assessed using the log-rank test. The survival results were considered significant when P < 0.05. All analyses were performed within the R 3.3.3 framework.

Correlation Between the Levels of Some DNA Methylation Sites in Human Normal Tissues and Chronological Age
Chronological age is referred to the real age of an individual, which is recorded by researchers. To study age-associated methylation levels in human over a wide range of ages and tissues, we collected 704 normal tissue samples from 14 to 89 years ( Figure 1A). We found that 145,523 DNA methylation sites across all tissues were correlated with age (a multiple testing corrected p < 0.05) ( Figure 1B). Further, we screened the correlations and we obtained 1,631 methylation sites that were significantly correlated with age, which were considered to be significantly correlated with chronological age (|cor| >0.3). Overall, we identified up to 1,609 positive correlations, accounting for 98.7% of the total number of methylation sites, which suggested that in most samples, the higher the level of age-related methylation, the older the organism in the sample ( Figure 1C). We classified tissue samples follow tissues to verify tissue specificity of the 1,631 age-related sites stated above. We performed correlation analysis with chronological age of each tissue, and we compared them with all-tissue correlation scores. The correlations of the most specific tissues were higher than the overall correlation level. This indicated that the selected sites were not only significantly correlated with age in all samples, but also more highly correlated with chronological age of samples in a single tissue ( Figure 1D).Therefore, we demonstrated that these 1,631 methylation sites were related to the chronological age of individuals not only in specific tissues but also in all samples.

Feature Selection and Model Construction of CancerClock
CancerClock model was constructed to predict the m-age in human normal samples based on adjacent normal tissue for each cancer (Figure 2A, see Methods). The LASSO linear regression model was used to screen the 1,631 age-related methylation sites in the training set using the GLMNET package of R. Finally, 282 sites were selected from 1,631 methylation site as model features ( Figure 2B, Supplementary Table S4). According to the least Mean-Squared Error (MSE, reflect the degree of difference between the estimator and true value, the smaller the MSE, the better the model fit) of LASSO linear regression model, we found that the value β 0 of the model was 34.63 when the adjustment parameter was 0.1419941 ( Figure 2C). Among the 282 methylation characteristics in the construction of model, the levels of cg08461576, cg05923914, cg27641628, cg13221458, cg05632420, and cg07103722 were significantly and negatively correlated with sample age, and the linear model coefficients of the six sites above were negative as well ( Figure 2D). We described the methylation level for these 282 methylation characteristics and we promoted them to effectively describe the chronological age of all normal samples through the heatmap (Figures 2A,D). The results showed that it was difficult to characterize the methylation age of samples by single methylation sites, and the methylation clusters might produce a better result. The scatter diagram shows the correlation between methylation sites (cg16867657, cg23606718) and chronological age in KIRC and BRCA. The overall background is the methylation sites of all samples (gray color), and the blue color represents the correlation between chronological age and methylation age in a certain type of samples such as BRCA and KIRC.
To evaluate the accuracy of the model, we used the training set and the test set to verify the model, separately. The results showed that the model responded well to the training set samples, and the correlation between chronological age and m-age was as high as 0.939 (P < 2.38e−295) ( Figure 2E). Meanwhile, for test set samples, the correlation between chronological age and m-age predicted by the model was up to 0.843 (P < 5.96e−20) (Figure 2E). The average accuracy of the model was 81%. To determine the difference between chronological age and predicted m-age in normal samples, we ordered the absolute value of the difference between the chronological age and predicted m-age in descending order, and we selected the quartile as the error value. The results showed that the age error value was 3 years, which indicated that the prediction error of over 75% of the sample age was controlled within 3 years, which is far The correlation between chronological age and predicted m-age in the training set and the correlation between chronological age and predicted m-age in the test set. lower than the age error predicted by previous researchers through telomeres.

Biological Processes and Phenotypic Traits of Age-Related Methylation Sites for CancerClock
Gene Ontology (GO) analysis of the genes for the 282 model characteristics revealed that these genes were associated with some GO terms such as "lysine catabolic process" (GO:0006554) and "lysine metabolic process" (GO:0006553) ( Figure 3A). A network consisting of GO terms and genes was constructed ( Figure 3B). These results indicated that age-related methylation sites could alter many important biology processes (Supplementary Table S5). The online analysis software of EWAS Atlas was used to analyze the biological characteristics enrichment of 282 selected features, and it was found that 128 of them were significantly enriched in age traits [-log10 (p) > 318], which was manifested as significant methylation quantitative trait loci (meQTLs) related to age. Meanwhile, 128 sites were enriched in age ontology entry (GO:0007568). In addition, we found that most methylation sites were enriched in human acute leukemia [-log10 (p) > 50] (Figure 3C, Supplementary Table S6). When we mapped these characteristics on the corresponding gene positions, we found that most of the loci corresponded to a single gene, while few loci corresponded to multiple genes ( Figure 3D). To describe the distribution of these features of CancerClock in genomic regions, we determined the position of these features in the genome. We found 37 of the 282 characteristics located in the promoter region of genes ( Figure 3E).

Predicted m-Age and Chronological Age of Cancer Patients
To depict the difference between m-age and chronological age in cancer samples, we utilized the CancerClock based on the methylation levels of all the 33 tumor types. We found that the m-age was generally different from the chronological age, but the degree of the difference was closely related to tumor types. Here, we ranked the absolute age differences between chronological age and we predicted the m-age and selected the quartile values as the age difference scores of tumors, separately ( Figure 4A, Supplementary Table S7). A number of tumors showed higher age difference score. Uterine carcinosarcoma (USC) had the highest age difference score, followed by ovarian serous cystadenoma carcinoma (OV) and uterine corpus endometrial carcinoma (UCEC). We found that the top three diseases with the largest age difference scores were all women tumors, suggesting that their chronological age is related to the pathological changes of the reproductive organ, and that female diseases affect chronological age to some extent. In addition, the distribution of m-age compared to the chronological age was described, separately ( Figure 4B). We found that the mage of different disease types was specific to their chronological age. For example, USC showed that all the 57 cancer samples were younger than the chronological age, and 96% of 54 testicular tumors showed age decline in men. According to our results, in the 121 and 184 samples with thymoma and ganglioma, 95% of the m-age was below the chronological age, respectively. In addition, in the 247 and 80 samples with renal cell carcinoma (RCC) and melanoma, 94% had younger m-age than chronological age, respectively.

Differential Levels of Some Methylation Sites Between m-Age and Chronological Age in Age-Related Cancer Samples
To explore further the difference between m-age and chronological age, we selected 25% age-related samples with the largest difference, and then performed differential analysis with normal samples in each cancer type (Supplementary Table S8).
In the age-related samples of cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), there were 77 methylation sites levels that were, respectively, different (FDR < 0.01, |log (FC)|>1) (Supplementary Table S9). Moreover, there were 38 methylation sites with significant changes in breast invasive carcinoma (BRCA) and 22 methylation sites with significant changes in kidney renal clear cell carcinoma (KIRC). The differential methylation sites of some cancers showed hypermethylation in tumor samples as follows: cervical squamous cell carcinoma (CESC) and bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), kidney renal clear cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), pancreatic adenocarcinoma (PAAD), sarcoma (SARC), and uterine corpus endometrial carcinoma (UCEC). It is suggested that hyper-methylation of these age-related differences may lead to changes in m-age and cancer development.

Association Between the Differential Methylation Sites in Age-Related Samples and Survival
Predicting survival state of cancer patients was critical and challenging (Zhou et al., 2017;Bao et al., 2019). To evaluate the influence of the methylation level on patient survival, cox risk regression model was performed for the differential methylation sites in age-related samples of each tumor type (Supplementary Table S10). The patients were divided into two groups according to median risk score. In most cancer types, these differential methylation sites in age-related samples were associated with survival ( Figure 5). The higher the risk scores of BRCA, COAD, esophageal carcinoma (ESCA), liver hepatocellular carcinoma (LIHC), and stomach adenocarcinoma (STAD), the lower the survival of the patients. The results indicated that differential methylation sites of age-related samples maybe could be an effective prognostic biomarker for cancers.

DISCUSSION
DNA methylation levels were closely associated with age and related to the occurrence of tumors . In this study, we explored the correlations between DNA methylation level and age in normal tissues. Importantly, we established CancerClock, a predictive model of age based on FIGURE 4 | The levels of some methylation sites were differential in age-related cancer samples between m-age and chronological age. (A) Age different scores between chronological age and predicted m-age among 33 cancer types. (B) The relationship between chronological age and predicted m-age for all the 33 disease types. Compared to the chronological age group, the pink color indicates that the predicted m-age group was down-regulated, the blue color indicates that the predicted m-age group was up-regulated, and the gray color indicates that the predicted m-age group remained unchanged.
DNA methylation level to depict the age clock in tumors. This predicted the age based on methylations of 282 sites from different tumor samples and it allowed us to access the m-age among methylation datasets. In addition, we found that the model was affected by a variety of biological processes, which may indicate the molecular influence of methylation on age. Age difference score in tumor samples showed the extent to which age affected different tumors. Meanwhile, methylation sites affecting the m-age of these tumors were identified. Through weighted survival analysis of cancer samples, we finally determined the impact of these age-causing sites on tumor survival.
CancerClock model is a multi-tissue age prediction model based on the methylation level of diverse cancer types in TCGA. It adopted a similar way as Horvath clock in building the age prediction model of normal samples (Horvath, 2013). More importantly, we applied this model to cancer samples and we depicted the difference between m-age and chronological age. These differences between m-age and chronological age may provide assistance to understand cancer development. Previous study also showed that epigenetic age acceleration is associated with colorectal cancer molecular characteristics and can be a significant predictor of overall survival, as well as age and tumor stage (Zheng et al., 2019). DNA methylation-based measures of biological age may be an important predictors of breast cancer risk (Kresovich et al., 2019). In the present study, we depicted this FIGURE 5 | The differential methylation sites in age-related samples were associated with survival. Kaplan-Meier survival analysis of two groups of patients with high-(blue line) and low-score (red line) groups. Survival days are shown along the X-axis. Overall survival rates are shown along the Y-axis. phenomenon in multiple cancer types and we tried to explain the mechanism by which methylation levels contribute to m-age and cause cancers.
Many studies used expression of gene and lncRNA to predict cancer development and prognosis (Zhou et al., 2015b(Zhou et al., , 2018Sun et al., 2019). In present work, m-age was considered as cancer biomarker for development and prognosis. CancerClock could be applied in predicting biological age for normal samples. The differences between biological and m-age were also could be identified. These differences could be used to explore the roles of methylation in cancer development and prognosis. In future work, more samples and experiments should be used to validate the present work.
In summary, the present study suggested that some methylation sites were associated with chronological age. Comprehensive age predicator CancerClock could predict m-age for normal samples and could find the differences between m-age and chronological age in age-related cancer samples. We further discovered the differential methylation sites between age-related cancer samples and normal samples. These differential methylation sites were associated with survival in cancers. In addition, the present study suggested that DNA methylation-based measures of chronological age might be important predictors of cancer risk.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
SN, LW, and HC conceived the project. TZ, YG, JW, and XL acquired the data. SS, YW, and SG analyzed the data. HZ, HL, and DS constructed the data resource. TZ, SN, and YG wrote the manuscript. All authors read and accepted the final version of the manuscript.


extfieldsofstudy: Physics Medicine Chemistry

text:
Reversible 2D networks of oligo(ε-caprolactone) at the air–water interface

Hydroxyl terminated oligo(ε-caprolactone) (OCL) monolayers were reversibly cross-linked forming two dimensional networks (2D) at the air–water interface. The equilibrium reaction with glyoxal as the cross-linker is pH-sensitive. Pronounced contraction in the area of the prepared 2D OCL films in dependence of surface pressure and time revealed the process of the reaction. Cross-linking inhibited crystallization and retarded enzymatic degradation of the OCL film. Altering the subphase pH led to a cleavage of the covalent acetal cross-links. The reversibility of the covalent acetal cross-links was proved by observing an identical isotherm as non-cross-linked sample. Besides as model systems, these customizable reversible OCL 2D networks are intended for use as pH responsive drug delivery systems or functionalized cell culture substrates.

Introduction
Poly(ε-caprolactone) (PCL) plays a paramount role as degradable implant material or as matrix for controlled drug delivery systems. It can be processed to nanoparticles, micelles, microspheres, electro spun mats, scaffolds or films [1][2][3][4]. End-group functionalized poly-or oligo(ε-caprolactone)s as telechelics are starting materials for building complex architectures. In such polymer systems, thermal properties and degradation profiles are adjustable by controlling the number of arms and chain lengths of the telechelics. With capability of further polymerization, one dimensional telechelics give new opportunities for creating three dimensional dendrimers, bulk polymer networks and 2D nanostructures [5,6]. Reversibly crosslinked polymer networks of PCL can be stimuli responsive. Temperature sensitive net-points can induce a shape memory response [7,8] while reversible chemical bonds promote self-healing characteristics [9,10]. Permanently cross-linked PCL networks are applied for tissue regeneration [11,12] and in printing technologies [13,14]. An important parameter in PCL networks modulated by cross-linking is ester bond cleavage by degradation [15][16][17][18][19]. The rate-determining step in PCL degradation by enzymatic hydrolysis generally slows down in the presence of netpoints formed by crystalline domains [20][21][22]. PCL based elastomers enzymatically hydrolyse completely over a wide range of temperatures at pH 7 when lacking crystalline domains. PCL elastomers combining controlled biodegradation and reasonable mechanical strength can be prepared by cross-linking a star-shaped hydroxyl-terminated poly(γ-methyl-εcaprolactone) with bis(β-lactone) [15]. Implants prepared by cross-linking PCL during ring-opening polymerization of caprolactone with bis(trimethylene carbonate) exhibit reduced crystallinity and formstability, but show suppressed enzymatic degradation [16]. Degradation kinetics are also influenced for benzoperoxide cross-linked PCL membranes in phosphate buffered saline [17]. Interestingly, PCL films prepared by photocross-linking with benzophenone or benzil exhibit faster hydrolysis due to enhancement in overall polarity of the film [18,19]. These examples show that the physico-chemical properties, mechanical stability and degradation profiles of 3D networks can be tailored in a controlled fashion by cross-linking. On the other hand, 2D polymer systems with interconnected telechelics in lateral plane can also lead to Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. generation of new material features as seen for 2D materials such as graphene and its derivatives, which modulate responses with cells, tissues or other biointerfaces [23,24]. Although graphene materials have immense potential in healthcare such as for drug delivery, bio imaging, capturing pathogens and tissue engineering, several contradictory results are reported in literature due to graphene shape, size, concentration, carbon source, and cytotoxicity [25,26]. Therefore, it is still a challenge to develop alternatives to conventional 2D materials for biological applications. We hypothesize that functionalized 2D polymer networks based on OCL telechelics with different chemical functionalities can be used as building blocks for creating bioresponsive polymer 2D networks. Besides network stability, the degradability should be altered by an external stimulus such as by pH and has to be systematically investigated. Chemical functions sensitive to pH are of specific interest, since gradients of pH exist at different physiological levels and in nature. Sites of tumour and inflammatory tissues have a slightly acidic pH when compared to their surrounding healthy tissues. Moreover, the environment of endosomes and lysosomes that carry substances inside cells is acidic. Hence, systems with pH sensitive chemical cross-links have attracted a key focus in developing drug delivery systems and in tissue regeneration [16,[27][28][29]. There are several possibilities to have OCL cross-links sensitive to pH, such as by using ortho ester [30], cis-aconityl bond [31], Schiff base formed between hydrazine and ketone groups [32] and acetal cross-links [33].
The Langmuir technique is ideally suited for preparing 2D materials of non-water or partially water soluble precursors. It allows for positioning of reactive end groups of monomers at the air-water interface for linking by different reaction mechanisms such as by irradiation [34], enzyme [35], catalyst [36], metal ions [37,38] and 'click' reactions [39]. Additionally, faster analysis of biodegradation processes are possible by this technique at molecular level [40]. While various ideally symmetric, small molecules are linked using Langmuir technique to form 2D materials [34,37,39,[41][42][43][44], only few studies described 2D cross-linking of oligomers or polymers using the Langmuir technique in recent years [5,38,[45][46][47][48][49]. For example, a semiinterpenetrating oligomeric network is achievable by embedding cross-linked methacryloxypropyl-terminated polydimethylsiloxane in a cellulose acetate butyrate matrix [5]. In our concept for preparing biodegradable 2D polymer networks, we use OCL telechelics end-capped with hydroxyl groups, because of their stability and well-studied behaviour at the airwater interface [50]. Acetal cross-links in hydroxyl terminated OCL monolayers are favoured since they are well described, can be used with different types of hydroxyl groups [51] and form under the conditions of the Langmuir technique experimental setup (aqueous subphase at room temperature). Acetal bonds form reversibly and provide the required pH sensitivity for a stimulus responsive material. This is a key advantage over traditional 2D materials forming irreversibly. The lateral size of 2D sheets is limited due to the nearly inevitable formation of grain boundaries when several sheets grow simultaneously next to each other. With reversible reactions, the grain boundaries can be cured and large lateral sizes can be achieved. The chemical equilibrium determines the number of unreacted chain ends that are available for further functionalization. As a proof of principle, the dialdehyde glyoxal is used, which imparts hemiacetal and acetal links in the hydroxyl terminated OCL polymer system and can react with four hydroxyl terminated chain ends of OCL. An idealised acetalization reaction scheme is shown in scheme 1.
The carbonyl group of glyoxal protonates under water, which is attacked by the nucleophilic hydroxyl group of an OCL molecule [52,53]. After subsequent deprotonation of the intermediate species, hemiacetal bonds are formed. Hemiacetal is an unstable intermediate and can react with hydroxyl groups of an OCL chain present in the vicinity to form an acetal [51]. The acetal bonds are stable in alkaline pH while under acidic pH, cleavage of some acetals occurs to release free end groups. The density of OCL hydroxyl groups at the air-water interface can be influenced by modulating OCL molecule packing density. Here, we examine two approaches to form a covalent 2D network. With linear, hydroxyl terminated OCL diol, a network is formed because each glyoxal molecule can connect up to four linear molecules and act as a branching centre (scheme 1(a)). Here, a branching centre is a point in the network where more than two OCL chains are covalently connected. When eventually more chain-ends are connected in the branching centres, the branched macromolecular entities link and ultimately form a 2D network. When using star shaped OCL tetrols, a certain number of branching centres are already provided by the tetrols (scheme 1(b)). Hence, we expect a faster network formation for tetrols, compared to diols. The maximum cross-link density for both linear molecules and tetrols is determined by molecular weight or arm length, respectively. Here, we use tetrols with an arm length (3200 g mol −1 ) that is similar to the molecular weight of the diol (2800 g mol −1 ), and therefore, in the diol monolayer, the hydroxyl end-groups concentration per unit area is double compared to the tetrol monolayer. It is known that PCL chain segments near covalent network points have limited mobility due to steric hindrance and hence their crystallization is hindered. We will use this effect to qualitatively evaluate the progress of the cross-linking reaction by studying the surface pressure induced crystallization of the layer. Furthermore, cross-links delay the biodegradation of OCL by reducing the number of free chain ends where small fragments can be generated relatively fast, so the enzymatic degradation of OCL is studied to obtain further information on the degree of cross-linking.
To free OCL molecules from a cross-linked OCL 2D network, the formation of branching centres is reversed. This is achieved by breaking the acetal bonds. We expect that acetal bonds in the network can be cleaved to release free chain ends by reducing the pH. The reversion of the network formation after lowering of the pH is studied by π-A isotherms and pressure induced crystallization. The prepared pH triggered OCL 2D network aims in broadening applications of 2D crystallisable polymers, such as for forming pH triggered enteric drug delivery systems.

Materials
Oligo(ε-caprolactone) diol (OCL diol, trade name CAPA 2304, Solvay Caprolactones, Warrington, UK) was used without any further purification. Star-shaped polymer based on ε-caprolactone with four-arms (4 arm OCL) was synthesized using the catalyst dibutyltin oxide and the initiator pentaerythritol [54]. Glyoxal (40% w/v, Sigma Aldrich, Germany) used as the cross-linker was dissolved in up to 5 ml water in required concentration before injecting into the subphase. The aqueous subphase and water for dissolution of glyoxal was obtained by a Milli-Q Gradient Scheme 1. Reaction schemes of 2D network formation of (a) OCL diol and (b) OCL tetrol with glyoxal. Encircled area represents a branching centre. By progressive connection of the branched macromolecular entities, a 2D network is formed.
A-10 water purification system (Millipore, 18.2 MU cm, toc<4 ppb). For cleaning of trough, ethanol (HPLC grade; Bernd Kraft, Duisburg, Germany) and chloroform (HPLC grade, Roth) were used. For the used OCL diol and 4 arm OCL, the characterisation data were reported previously [55]. The average number molecular weight (M n ), determined by gel permeation chromatography (GPC) was 2800 g mol −1 and 13 000 g mol −1 , respectively. The determined M n (GPC) per end group for OCL diol and 4 arm OCL was 1400 g mol −1 and 3200 g mol −1 , respectively.

Langmuir monolayer technique
The surface pressure-area (π-A) isotherms were recorded with a Langmuir trough (KSV NIMA, Finland) placed on an active vibration isolation system (Halcyonics variobasic 40, Accurion, Germany) within a laser safety cabinet. The surface pressure (π) was measured by the Wilhelmy technique with a calibrated sensor located in the centre between the barriers. Water level compensation system was used to compensate the loss of water due to evaporation during the course of experiment to avoid false changes in the recorded surface pressure. The trough was cleaned thoroughly with ethanol followed by chloroform. Afterwards, the deionized water was filled and removed at least three times. By monitoring the surface pressure, the purity of the trough and the subphase was controlled. In the final aqueous subphase, while closing the barriers, the total change of the surface pressure was below 0.2 mN m −1 . For all OCL samples, the chloroform stock solutions had a concentration in the range of 0.2-0.4 mg ml −1 . The solution was applied drop-wise onto the air-water interface using a microsyringe (Hamilton Co., Reno, NV, USA). The chloroform was allowed to evaporate for ten minutes before the start of compression. The rate of compression and expansion of monolayer was 10 mm min −1 . The surface pressure was recorded as a function of the mean molecular area per repeating unit (MMA).
The elasticity modulus was determined using the formula, κ=A −1 (δπ·δA −1 )) where A is the area per molecule and π is the surface pressure. At the maximum of elasticity modulus, the corresponding surface pressure and mean molecular area is π i and MMA i , respectively.

Cross-linking of OCL monolayer
For cross-linking by dialdehyde, OCL molecules were compressed and withheld at a constant surface pressure. Thereafter, glyoxal (Sigma Aldrich, Germany; vapour pressure of glyoxal 40% at 20°C: 24 mbar) was injected into subphase (glyoxal/caprolactone molar ratio in the subphase=20 000) and the monolayer was withheld at a specific constant surface pressure. Glyoxal is highly water soluble and its hydrate forms, which are considered non-volatile, contain reactive carbonyl groups [56]. After the fixed period of reaction time, the barriers were opened completely to initial position of the Langmuir trough and π-A isotherms of the cross-linked layers were recorded.

Degradation by lipase
For enzymatic degradation of OCL 2D networks, lipase (Sigma, Germany) from Pseudomonas cepacia (0.007 mg ml −1 ) was injected into the subphase. For degradation of the control sample, glyoxal was also added just before the injection of the enzyme into the subphase. A t ·A 0 −1 , as function of time was recorded where A 0 is the initial surface area occupied by Langmuir film at t=0 and A t is the acquired surface area after a certain time interval, t. The degradation experiment of Langmuir film was performed at 7 mN m −1 ; 22±0.5°C.

Brewster angle microscopy (BAM) studies
BAM images were recorded with ellipsometer nano-film_ep3 (Accurion, Gottingen, Germany) equipped with a high performance CCD camera, a 10× magnification lens with a maximum lateral resolution of 2 μm, and a 658 nm class IIIB laser source.

Interfacial infrared spectroscopy
Polarization Modulation Infrared Reflection Adsorption Spectroscopy (PM-IRRAS; KSV NIMA) was used together with a Langmuir trough (KSV NIMA) with water level compensation to study the change in chemical composition of film under spectral range of 800-4000 cm −1 .

Langmuir film isotherms
The π-A isotherm of a non-cross-linked OCL diol monolayer was recorded first to evaluate the suitable conditions for the acetalization reaction. It shows the typical Langmuir monolayer behaviour reported for PCL [57], and is marked by a sudden reduction of surface pressure (collapse pressure~11.2 mN m −1 ; figure 1). Based on the values of static elastic modulus of control (non-cross-linked) sample of OCL diol monolayer, acetalization reaction was performed at three different surface pressure values 2, 5 and 7 mN m −1 , corresponding to the increasing concentration of OCL end groups per area. To avoid hydrolysis of ester bonds, mild reaction conditions (pH 5.7; room temperature) were realized. After the reaction, the barriers were opened to initial position of the trough and π-A isotherms were recorded to compare reacted OCL 2D film with the control sample.
Apparently, compared to the control sample, lowering of mean molecular area, an increased π i (table 1) and an elevated collapse pressure was observable for all cases of OCL diol monolayer reaction. Rather, a clear effect on relative reduction of MMA i was observed for all OCL 2D films prepared depending on end-groups concentration (figure 1). Such a contraction of a monolayer at the air-water interface by cross-linking was also reported for cross-linked PDMS 2D networks [5]. For OCL 2D films, a twofold decrease in MMA i compared to the control sample was observed only after OCL diol monolayer was reacted at a high packing density of 7 mN m −1 for 18 h. The acetalization reaction at air-water interface is slow and depends on the OCL molecule end-group concentration. Thus, it is expected that increasing the end-group concentration promotes cross-linking via acetal formation [58]. After reaching the collapse surface pressure, OCL 2D networks formed during 18 h reaction at 5 mN m −1 and 5 h reaction at 7 mN m −1 showed no sudden decrease, rather a plateau region. Additionally, a collapse pressure in π-A isotherm was not observable until the end of compression for the OCL 2D film formed at the high packing density at 7 mN m −1 for 18 h reaction time. The shifting of collapse pressure of all OCL 2D films to a higher value compared to the control sample proves an increase in stability of the film due to acetalization reaction. To prove the acetalization reaction, Polarization Modulation-Infrared Reflection Absorption Spectroscopy (PM-IRRAS) spectra under different reaction times were recorded. The spectrum (figure 2(a)) before injection of glyoxal had no discernible peaks and was used as the reference spectrum, i.e. all spectra were normalized with respect to the spectrum at t=0. The PM-IRRAS signal of the OCL 2D network for visualizing acetal links was weak owing to molecular level ultrathin film. Nevertheless, the acetal doublet peaks (1060 and 1100 cm −1 ) were clearly distinguishable. The peaks at 1100 and 1060 cm −1 were absent in the non-cross-linked OCL diol monolayer ( figure 2(a)) and became visible in all spectra recorded after 21 h glyoxal reaction ( figure 2(b)) with the monolayer. The bands are characteristic for acetals and attributed to symmetric COCOC (1100 cm −1 ) and antisymmetric COCOC (1060 cm −1 ) stretching.
Further evidence of occurrence of cross-linking was demonstrated by 4 arm OCL. Compared to 2D film of OCL diol, a~2 fold decrease in MMA i was reached in 2D film of 4 arm OCL in a relative less

Crystallization behaviour of OCL 2D network
The crystallization behaviour of OCL in restricted geometries is well studied and is affected by molar mass, film thickness, temperature, the field for chain diffusion and many other factors [50,57,59]. At the constant surface pressure of 10.5 mN m −1 , which is between π i (6.9 mN m −1 ) and the collapse pressure (11.2 mN m −1 ), the non-cross-linked OCL monolayer (control) shows an area reduction (figure 4(a)) due to formation of typical small sized 3D domains under BAM ( figure 5(a)). These 3D domains of OCL seen as bright regions in the BAM image are the spherulites, which are surrounded by the amorphous layer (dark region). To study the crystallization behaviour of OCL 2D networks, we cross-linked OCL monolayers at 7 mN m −1 (in absence of spherulites) for different reaction times and then compressed the layers to the elevated surface pressure of 10.5 mN m −1 . Introduction of branching structure in PCL bulk systems lowers crystallization ability and T m (melting temperature) while it increases the mechanical stability [60][61][62]. Since random cross-linking of OCL molecules in Langmuir monolayer presumably leads to branching of OCL chains to form a 2D network, we expect that a retardation of the crystallization process can provide an indirect proof of cross-linking of OCL monolayer.
For OCL 2D networks formed at 5 h reaction time at 7 mN m −1 , the normalized area reduction at 10.5 mN m −1 was decreased compared to the control sample ( figure 4(b)), with accompaniment of large sized spherulites observed in BAM ( figure 5(b)). Such an increase in lateral dimensions can be asserted to an increased molar mass of the molecules after crosslinking and therefore such large spherulites, typically formed for intermediate molar mass (~8 kDa) of OCL molecules were observed [57].
On the contrary, the relative area reduction at 10.5 mN m −1 was minimal for the OCL 2D network formed at 18 h reaction time at 7 mN m −1 ( figure 4(b)), with no visible crystals in BAM images ( figure 5(c)). This strongly projects that cross-linking impedes the surface pressure induced lateral crystallization of OCL monolayer. By connection of OCL chains by acetal cross-links, molar mass in the monolayer increases while due to formation of cross-linking netpoints, the OCL chain mobility is decreased leading to a retarded crystallization. However, lateral surface pressure induced crystallization in ultrathin films is a dynamic process and cannot be directly interpreted  and therefore can be seen only as an indirect proof of cross-linking. A further insight into the crystallization behaviour is available by recording the expansion π-A isotherms (figure 5). Expansion isotherms formed by opening the Langmuir trough barriers after the end of compression at a constant rate indicate the melting behaviour of surface pressure induced spherulites. The spherulites of the control sample dissolve during the expansion of the monolayer, where the polymer chains of the spherulites continuously unfold to render initial '2D' monolayer state. This marks a faster initial decrease in surface pressure in the expansion isotherm followed by a plateau around 8.5 mN m -1 ( figure 5(a)). Such a similar progression of the expansion isotherm was also seen by OCL 2D networks formed after 5 h cross-linking at 7 mN m −1 ( figure 5(b)). However, for OCL 2D network formed after 18 h cross-linking at 7 mN m −1 , a plateau in the expansion isotherm was not observable, while the curve follows a similar trend as its compression cycle ( figure 5(c)). Figure 6 shows enzymatic degradation isotherms of OCL 2D networks formed under different crosslinking reaction durations at the surface pressure of 7 mN m −1 , in the absence of any 3D domains in the   layer. The catalytic triad of the enzyme lipase cleaves OCL chains into water-soluble fragments, causing monolayer area reduction. For control sample, after injection of lipase into the subphase (time=0), a fast degradation was detected. It was previously shown that OCLs below molar masses of 10 000 g mol -1 present identical degradation isotherms [55]. On the contrary, the degradation isotherm for OCL 2D networks shifts substantially compared to the control sample. The slowest degradation rate was observed for OCL 2D networks formed after 18 h reaction time at 7 mN m -1 . Such a dense 2D network with modified and interconnected end groups of OCL molecules possibly hinders interaction of enzyme with OCL segments or has slow formation of water-soluble OCL fragments. Figure 7 shows the behaviour of OCL diol monolayer at acidic pH of 2.5 with and without the presence of glyoxal in the subphase. A slower area reduction of OCL diol monolayer at pH 2.5 in presence of glyoxal suggests an interplay of acidic hydrolysis of ester bonds and the cross-linking reaction phenomenon. To prove that the acetal cross-links in the OCL 2D network are reversible, the reverse cross-linking reaction was promoted for retrieval of free OCL molecules by altering the pH and π-A isotherm ( figure 8(a)) was recorded.

Reversibility of the cross-links
OCL 2D network was formed at 7 mN m -1 after 18 h of glyoxal reaction of OCL diol monolayer in aqueous subphase of pH 5.7. To carry out the reverse acetalization reaction, the pH was shifted to 2.5 for the next 9 h. The barriers were opened to initial Langmuir trough position and compression-expansion isotherms were recorded ( figure 8(a)). The π-A curve of OCL 2D network after alteration to acidic pH ( figure 8(a)) deviates from the π-A curve of OCL 2D network formed at pH 5.7 ( figure 5(c)). Indeed, it follows an identical slope with the control sample (figure 5(a)), with observance of collapse pressure at 13.5 mN m -1 . Change in crystallization behaviour further validates the reversibility of cross-links in the network. In contrast to the BAM image of OCL 2D network which shows no 3D domains at 10.5 mN m -1 ( figure 5(c)), the BAM image for OCL 2D network that underwent reverse acetalization reaction ( figure 8(b)) shows occurrence of relatively larger spherulites comparable to the ones observed in non-cross-linked OCL layer at 10.5 mN m -1 ( figure 5(a)). This implies that alteration to acidic pH caused a decrease in the covalent acetal link density of the cross-linked OCL 2D network, allowing OCL chains to fold and participate again in the lateral crystallization process.
The similar compression-expansion π-A isotherm behaviour as the control sample and occurrence of spherulites implies an increase in the conformational freedom in OCL chain segments and enhanced crystallization degree in the film by the decrease of acetal cross-link density. Therefore, acetal links orchestrated the forming of a reversible OCL 2D network that is reversible by means of pH.

Conclusion
The 2D cross-linking of hydroxyl end-capped OCL monolayer by glyoxal at the air-water interface was demonstrated by changes in π-A isotherms, film crystallization behaviour and film degradation rate. Additionally, interfacial infrared spectroscopy demonstrated the acetal bond formation in the OCL 2D film. After cross-linking, π-A isotherm displayed a higher reduction in the mean molecular area for OCL 2D networks formed under high OCL molecular packing density and longer reaction time. Moreover, π-A isotherms demonstrated that reduction in mean molecular area is achieved faster for the 2D network based on tetrafunctional OCL tetrol compared with Figure 8. Compression-expansion π-A isotherm of OCL 2D network that underwent reverse acetalization reaction by the alteration to acidic pH at 22±0.5°C. Dotted lines represent the expansion isotherm and the BAM image (scale=100 μm) is taken at the end of compression cycle. bifunctional OCL diol. Reversibility of acetal crosslinks in the OCL 2D network is demonstrated by altering the subphase pH to promote the reverse reaction. Such a film exhibited comparative crystallization and π-A isotherm of non-cross-linked OCL monolayer.
Such 2D networks of polymers can give fundamental insights on physical chemistry and rational design of bulk polymer networks. Further work on covalent functionalization of unreacted end-groups will explore how such a stimuli responsive 2D polymer network can be used for delivery of bioactive substances, as coating for medical devices and as functionalized substrates for studying cell behaviour.