Functional Profiling of Chondrogenically Induced Multipotent Stromal Cell Aggregates Reveals Transcriptomic and Emergent Morphological Phenotypes Predictive of Differentiation Capacity.

Multipotent stromal cells (MSCs) are an attractive cell source for bone and cartilage tissue repair strategies. However, the functional heterogeneity of MSCs derived from different donors and manufacturing conditions has limited clinical translation, emphasizing the need for improved methods to assess MSC chondrogenic capacity. We used functionally relevant morphological profiling to dynamically monitor emergent morphological phenotypes of chondrogenically induced MSC aggregates to identify morphological features indicative of MSC chondrogenesis. Toward this goal, we characterized the morphology of chondrogenically stimulated MSC aggregates from eight different human cell-lines at multiple passages and demonstrated that MSC aggregates exhibited unique morphological dynamics that were both cell line- and passage-dependent. This variation in 3D morphology was shown to be informative of long-term MSC chondrogenesis based on multiple quantitative functional assays. We found that the specific morphological features of spheroid area, radius, minimum feret diameter, and minor axis length to be strongly correlated with MSC chondrogenic synthetic activity but not gene expression as early as day 4 in 3D culture. Our high-throughput, nondestructive approach could potentially serve as a tool to identify MSC lines with desired chondrogenic capacity toward improving manufacturing strategies for MSC-based cellular products for cartilage tissue repair. Stem Cells Translational Medicine 2018;1-12.

Despite their regenerative potential, the functional heterogeneity of multipotent stromal cells (MSCs) derived from different donors and manufacturing conditions has limited clinical translation, emphasizing the need for improved methods to assess MSC differentiation ability. Utilizing functionally relevant morphological profiling to dynamically monitor emergent morphological phenotypes of chondrogenically induced MSC aggregates, this work demonstrated that MSC aggregates derived from different donors showed unique morphological dynamics that were informative of long term chondrogenic differentiation capacity. Furthermore, transcriptomic analysis revealed that the gene expression of undifferentiated MSCs also correlated with chondrogenic potential. Overall, a well‐defined tool is presented for the early estimation of MSC chondrogenesis.

Introduction

Articular cartilage is a dynamic, stress‐bearing connective tissue that functions to facilitate the frictionless articulation of adjoining bones within synovial joints. The avascular, but highly hydrated, cartilage matrix comprises a hierarchical structure of proteoglycans and collagens organized to provide its characteristic viscoelastic and compressive mechanical properties 1, 2. The cellular constituents of cartilage are chondrocytes, which are specialized extracellular matrix‐producing cells with low proliferative capacity 3. However, damaged articular cartilage (due to trauma or disease) remains one of the most difficult tissues to repair due to its low endogenous capacity for regeneration and continues to be a leading cause of disability worldwide 4. Furthermore, the successful outcome of cell‐based therapies for cartilage repair requiring the use of chondrocytes, such as autologous chondrocyte implantation (ACI), is often marred by the propensity of chondrocytes to dedifferentiate during ex vivo expansion 5 as well as the significant donor morbidity associated with their harvest.

To mitigate such problems associated with ACI, multipotent stromal cells (MSCs) have become an attractive alternative cell source due to their capacity for self‐renewal, ease of isolation, and chondrogenic differentiation potential 6, 7. Unlike adipogenesis and osteogenesis, the chondrogenic induction of MSCs typically requires the use of a 3D culture format allowing for cell aggregation 8. Many research groups have successfully induced chondrogenic differentiation of MSCs in vitro via high‐density 3D culture platforms, which allow the formation of 3D microtissues by promoting intimate cell–cell contacts 9, 10, 11. This 3D microtissue model is thought to recapitulate aspects of developmental processes during mesenchymal condensation that precede cartilage formation 12, 13. Although widely recognized, this biologically inspired design has only recently been incorporated into tissue engineering strategies for enhancing cartilage repair in vivo 14, 15, 16. While abundant research demonstrates their therapeutic potential in the context of cartilage repair, a major challenge to successful clinical translation that remains is MSC functional heterogeneity 17. Significant cell‐to‐cell heterogeneity exists within MSC populations derived from different donors and manufacturing processes, which results in substantial differences in functional capacity and has thus limited their overall therapeutic effectiveness and clinical potential 18.

In order to better understand and address MSC functional heterogeneity, functionally relevant morphological profiling (FRMP) has been previously used to identify morphological features predictive of MSC functional capacity 19. Indeed, effective analytical tools such as FRMP have successfully linked phenotypic readouts such as cell morphology to MSC osteogenic potential 20, 21. Under osteogenic stimulation in 2D, MSCs displayed more uniform actin fiber rearrangement, greater cell spreading area, decreased aspect ratios, as well as increased focal adhesions 22. These morphological and cytoskeletal “profiles” of MSCs could then be used to predict the capacity of an MSC line to undergo osteogenic differentiation. An application of this methodology demonstrated that certain morphological features of MSCs exposed to osteogenic stimuli were highly predictive of their ability for eventual mineral deposition 20. It was further revealed in another investigation that distinct morphological changes of MSCs exposed to the inflammatory cytokine interferon‐γ were highly correlated to MSC ability to suppress T cell activation, highlighting the utility of using functionally relevant stimuli for developing predictive bioassays 23. While the morphological profiling of single cells in 2D can reveal heterogeneous responses to external stimuli, there still remains significant need for quantitative methods to assess MSC functional heterogeneity within physiologically relevant 3D model systems. The adaptation of a functionally relevant approach for the estimation of MSC chondrogenic capacity via the interrogation of their 3D aggregate morphology would provide a useful tool for bridging the gap between understanding MSC heterogeneity and clinical translation for cartilage repair.

Here, we use a nondestructive approach to quantify high‐dimensional morphological features of chondrogenically stimulated MSC aggregates and show differences in morphological dynamics related to cell lines from different donors and cell line passage. We hypothesize that such differences in MSC aggregate morphology, which represents a functionally relevant 3D phenotypic readout, are early indicators of chondrogenic differentiation capacity of MSC lines. By establishing linear relationships between metrics of MSC aggregate morphology and multiple quantitative chondrogenic outcomes, we show that emergent 3D aggregate phenotypes correlate strongly with long‐term chondrogenesis assessed through cell content, gene expression, and protein expression. We also identified a number of genes from the transcriptomic analysis of undifferentiated MSCs that correlated with their chondrogenic differentiation capacity. These findings lay a foundation for better understanding MSC functional heterogeneity in the context of chondrogenesis and represent a generalizable analytical approach for quickly assessing cell quality and guiding cell manufacturing strategies that could be applied to additional cell‐types and therapeutic applications.

Materials and Methods

MSC Culture and Expansion

Human bone marrow‐derived MSCs were obtained from eight different donors purchased from RoosterBio (RB9, RB12, RB14, RB16; Frederick, MD) at passage 1 (P1), Lonza (8F560, 167696; Walkersville, MD) at passage 2 (P2), and AllCells (PCBM1632, PCBM1641; Alameda, CA) at P2 (see Fig. 1 for donor specifications). MSCs were cultured and expanded following previously established procedures 20. P2 RoosterBio MSCs and P3 Lonza/AllCells MSCs were used to represent early passage (EP) groups as the MSCs were derived after one passage from the indicated passage at purchase. P5 MSCs for all groups were used to represent late passage (LP).

Figure 1

A schematic of the experimental workflow for the analysis of MSC chondrogenic capacity is shown. Using MSC lines derived from eight different donors at two passages each (n = 10–15), we aim to primarily investigate the effects of cell passage on the morphological dynamics of formed MSC aggregates, how these morphological characteristics correlate with MSC chondrogenic capacity, and whether these effects are conserved across MSCs derived from multiple donors. MSCs at P2 or P3 and P5, respectively designated as early and late passage, are culture expanded for 72 hours and then formed into aggregates of defined size for 3D culture over 21 days. Images of aggregates are acquired and tracked during the culture period at various time points for morphological analysis using quantitative computational methods (CellProfiler). MSC aggregate morphological features (14 total) are then correlated with their respective chondrogenic outcomes, which are measured at day 21. Abbreviation: MSCs, multipotent stromal cells.

MSC Spheroid Culture and Chondrogenic Differentiation

To generate MSC aggregates for spheroid culture, MSCs derived from each donor were trypsinized and then resuspended in chondrogenic induction medium (high glucose Dulbecco's MEM supplemented with 10 ng/ml TGF‐β3 (Peprotech, Rocky Hill, NJ), ITS + Premix (6.25 μg/ml insulin, 6.25 μg/ml transferrin, 6.25 μg/ml selenious acid, 5.35 μg/ml linoleic acid, and 1.25 μg/ml bovine serum albumin) (Corning, Corning, NY), 50 mg/l ascorbic acid‐2‐phosphate, 10−7 M dexamethasone, 1% l‐glutamine, and 1% penicillin–streptomycin) at a concentration of 105 cells/150 μl. 150 μl of cell suspension was added to each well of ultralow‐attachment 96‐well plates (Costar; Corning) and then centrifuged at 250g for 5 minutes. MSC aggregates were incubated in parallel in a humidified chamber at 37°C and 5% CO2 and cultured for 3 weeks. The chondrogenic induction medium was changed at 24 hours and then at every other day following aggregate formation. At day 21, samples were retrieved for analyses. A minimum of 10 individual spheroids were generated for each experimental group.

Quantitative Morphological Analysis of MSC Aggregates

MSC aggregates (Fig. 1) were imaged at days 1, 4, 7, 11, 14, 18, and 21 during aggregate culture via light microscopy (Olympus CKX41 with Nikon DS‐Fi2 camera and Nikon DS‐L3 attachment). Automated analysis and quantification of morphological features were performed using the open‐source image analysis software CellProfiler 2.2.0 to provide high‐dimensional morphological characterization of MSC aggregates comprising 14 shape and size features (Supporting Information Table S1). The CellProfiler algorithm (pipeline) used to perform the morphological characterization can be viewed in Supporting Information Table S2. Principal component analysis (PCA) was performed to reduce these large multivariate morphological datasets into a smaller number of artificial (dimensionless) principal components (PC), which represent an uncorrelated linear combination of optimally weighted observed variables that capture maximal variance in the original dataset. From this PCA, only the first two PC were selected as they accounted for at least 70% of the cumulative explained variance of the original data.

Quantitative Evaluation of Chondrogenic Capacity

At day 21, MSC aggregates were collected from spheroid culture and frozen at −20°C until use for biochemical analysis. More detailed descriptions of this protocol are found in the Supporting Information.

DNA content, used as an indirect indicator of aggregate cellularity, was quantified using the Quant‐iT PicoGreen dsDNA Assay Kit (Molecular Probes, Eugene, OR) according to the manufacturer's protocols. Sulfated glycosaminoglycan (GAG) content, used as a direct indicator of extracellular matrix deposition and synthetic activity, was measured using the established 1,9‐dimethylmethylene blue dye calorimetric assay as previously described 24. To measure the synthetic activity of MSC aggregates, the resulting total GAG amounts from each sample were normalized to the total DNA amount from that same sample.

MSC aggregates were processed for real‐time reverse transcriptase polymerase chain reaction (RT‐PCR) to quantify the expression of selected genes as gauges for chondrogenic differentiation. Total RNA were isolated from samples using the RNeasy Mini Kit (Qiagen, Valencia, CA) following the manufacturer's protocols. The isolated RNA samples were reverse‐transcribed into cDNA using the Verso cDNA Synthesis Kit (ThermoFisher Scientific, Waltham, MA), where the final cDNA transcripts were used for the quantification of selected genes listed in Supporting Information Table S3 25, 26 by RT‐PCR using established primers and Power SYBR Green Master Mix (Applied Biosystems, Foster City, CA).

PCA was performed on supervised datasets measuring chondrogenic outcomes to generate composite scores principle component 1 (PC1) that captures maximal variance in the data representing matrix synthetic activity (SynthAct; DNA, GAG, GAG/DNA) and chondrogenic gene expression (ChondroGene; type IIα1 collagen, aggrecan, and type IIα1/type I collagen fold ratio) of the MSCs. Unsupervised PCA was also performed on all available chondrogenic outcome measures to generate a composite metric that can be used as an overall chondrogenic index (ChondroInd).

Histology

More details regarding the processing of samples for histology are found in the Supporting Information. Raw histology images were then quantified using the CellProfiler algorithm shown in Supporting Information Table S4. Briefly, the total intensity quantified from the histological stains were normalized to the corresponding MSC spheroid area to account for differences in aggregate size.

Gene Expression Microarray Hybridization, Data Acquisition, and Statistical Analysis

All gene expression microarray data files have been uploaded into the public repository, Gene Expression Omnibus (GSE110755). A more detailed description of the experimental and analytical methods used in this study are found in the Supporting Information.

Results

Chondrogenic Synthetic Activity and Gene Expression Exhibit Cell Line‐ and Passage‐Dependence

To assess the chondrogenic differentiation potential of each MSC line, we quantified cellularity (DNA content), matrix production, and chondrogenic gene expression after 21 days of chondrogenic stimulation. Regarding the DNA content, EP MSC aggregates generally maintained higher cellularity when compared to LP MSC aggregates (Fig. 2A). Differences between MSC cell lines were more evident in terms of extracellular matrix production (as measured by GAG content). Despite exhibiting greater matrix production overall versus LP aggregates, not all EP MSC aggregates (8F3560, 167696, and PCBM1632) displayed high amounts of GAG deposition by D21 (Fig. 2B). Among LP aggregates, only RB9 and RB16 showed moderate amounts of GAG deposition. This pattern in GAG deposition became much more apparent when normalized to DNA. Of note, the normalized matrix production (GAG/DNA) of LP RB9 aggregates reached levels comparable to EP groups with high GAG production (Fig. 2C).

Figure 2

Following 21 days of 3D culture, the (A) DNA, (B) GAG, and (C) GAG/DNA were quantified using biochemical assays to provide an indirect measure of cellularity, matrix deposition, and matrix synthetic activity, respectively. Additionally, the expression levels of canonical chondrogenic genes, which include (D) type IIα1 collagen, (E) aggrecan, (F) type I collagen, (G) type X collagen, and the (H) type IIα1/I collagen fold ratio were also quantified. * Indicates a statistically significant difference between early passage and late passage samples (*p < .05, and ****p < .0001, two‐way ANOVA). Data are presented as the mean ± SD values of the biological replicates included per group (n = 4) for these analyses. (I): Histological sections of multipotent stromal cells aggregates illustrate GAG (blue) and collagen (red) distribution after 21 days in 3D culture. Scale bar = 500 μm.

Gene expression profiling of well‐known chondrogenic marker genes revealed trends between cell‐lines and passages that deviated from the DNA and GAG content (Fig. 2D–2G). For instance, at EP, both 8F3560 and PCBM1641 aggregates showed the highest levels of type IIα1 collagen expression despite the former having produced the least amount of GAG (Fig. 2D). Other such discrepancies between synthetic activity and gene expression also appeared for the expression of aggrecan (Fig. 2E). For type I collagen expression, a cartilage dedifferentiation marker, LP aggregates generally maintained greater levels of expression compared to EP aggregates (Fig. 2F). MSC aggregates with low matrix production from both EP and LP groups displayed significantly increased type X collagen gene expression, a hypertrophic marker, when compared to MSC lines of higher matrix production (Fig. 2G). When the type IIα1/I collagen fold ratio (Fig. 2H) was measured to indicate the extent of hyaline (as opposed to fibrous) cartilage differentiation, we observed that the pattern of this metric between groups generally paralleled that of their matrix production (Fig. 2I).

To further characterize MSC chondrogenesis, extracellular matrix deposition, and distribution was quantified using image analysis of histological sections of aggregates (Fig. 2I). At EP, six out of the eight MSC lines (RB9, RB12, RB14, RB16, 8F3560, and PCBM1641) displayed moderate‐to‐significant GAG deposition throughout the entire spheroid as indicated by the intense Alcian blue staining. Only two of these MSC lines, namely RB9 and RB16, showed histologically detectable GAG deposition at LP. Overall, quantified GAGs from the histology correlated very strongly with biochemically detected GAGs (Supporting Information Fig. S1A). Most MSC lines display moderate‐to‐high amounts of collagen deposition at both EP and LP, especially at the periphery of the chondrogenically stimulated aggregates. Although EP aggregates produced more total collagen than LP aggregates, the quantified collagen from the histology yielded no apparent patterns between MSC lines and did not correlate with biochemically measured GAG matrix production (Supporting Information Fig. S1B).

Generation of a Composite Quantitative Metric for Overall Chondrogenic Capacity

As singular measures from individual assays may not fully elucidate the extent of MSC chondrogenic differentiation, we used PCA to integrate multiple assay outcomes into single quantitative values representing both overall and key aspects of chondrogenic potential. Supervised PCA was performed as outlined in Figure 3A, where the PC accounting for the maximum observed variance (PC1) from each respective PCA was used to construct quantitative composite measures of collective chondrogenic synthetic activity (DNA, GAG, and GAG/DNA) and chondrogenic gene expression (type IIα1 collagen, aggrecan, and type IIα1/I collagen fold ratio), referred to as the SynthAct and ChondroGene scores, respectively. Type I collagen and type X collagen were not included in the derivation of the ChondroGene scores as their inclusion resulted in decreased captured variance in the first two PC scores for the dataset (i.e., 83.4% vs. 62.9% of captured variance), making the use of a single composite score representative of chondrogenic gene expression difficult (Fig. 3, Supporting Information Fig. S2A). Furthermore, increases in type I collagen and type X collagen expression, unlike the other chondrogenic gene expression metrics, do not reflect hyaline‐like chondrogenic differentiation but rather dedifferentiation and hypertrophy, respectively. Overall MSC chondrogenic potential, referred to as the ChondroInd, is then measured by PC1 of the unsupervised PCA performed on all outcome measures (including type I collagen and type X collagen) and embodies the greatest amount of variance in the unsupervised dataset at 45%. When experimental results representing various aspects of MSC chondrogenic differentiation were plotted in PC space, we observed clear separations between experimental groups, indicating the use of raw PC scores as a viable method for assessing MSC chondrogenic potential.

Figure 3

(A): PCA was performed using supervised datasets comprising metrics that include the greatest variability in synthetic activity (DNA, GAG, and GAG/DNA) and chondrogenic gene expression (type IIα1 collagen, aggrecan, and type IIα1/type I collagen fold ratio) as illustrated by the colored boxes. Unsupervised PCA was also performed on all outcome measurements to derive an overall quantitative composite score (ChondroInd) that captures aspects of both the actual synthetic activity as well as chondrogenic gene expression of multipotent stromal cell lines. Normalized PC1 scores from each PCA performed (supervised and unsupervised), which account for the greatest variance in SynthAct (81.5%), ChondroGene (48.2%), and ChondroInd (45%) datasets, are shown in (B), (C), and (D), respectively. * Indicates a statistically significant difference between early passage and late passage samples (**p < .005 and ****p < .0001, two‐way ANOVA). Data are presented as the mean ± SD values of the biological replicates included per group (n = 4) for these analyses.

For improved clarity, the derived chondrogenic composite scores were then rescaled via normalization to remove negative values (Fig. 3B–3D). By compartmentalizing MSC chondrogenic differentiation into individual metrics for chondrogenic synthetic activity and gene expression, we observed that while EP MSC aggregates maintained much higher SynthAct scores over LP aggregates (Fig. 3B), the same pattern was not observed for ChondroGene scores (Fig. 3C). For SynthAct, MSC groups exhibited scores that matched the observed matrix deposition from the histological analysis. However, there were notable inconsistencies in the ChondroGene scores as some LP groups (PCBM1641 and 8F3560) displayed increased ChondroGene scores despite possessing correspondingly low SynthAct scores. Indeed, the SynthAct scores were only weakly correlated with the ChondroGene scores (Supporting Information Fig. S2B). Additionally, both histomorphometrically derived metrics for GAGs and collagen either weakly correlated or failed to correlate with chondrogenic gene expression, respectively (Supporting Information Fig. S2C). Despite this mismatch, the overall ChondroInd, which incorporates both SynthAct and ChondroGene, provided differences between MSC groups that were more congruent with their expressed phenotype (Fig. 3D), highlighting the utility of such a composite score for quantifying the degree of overall MSC chondrogenic differentiation.

Morphological Profiling of MSC Aggregates Reveal Cell Line‐ and Passage‐Dependent Differences in Aggregate Size and Shape Change

In our study, culture‐expanded MSCs (from eight different donors at EP and LP) were formed into aggregates and stimulated with standard chondrogenic induction medium to assess MSC chondrogenic capacity in tandem with morphological profiling, where high‐dimensional morphological features were quantified at various time points. Unsupervised PCA was performed (using all 14 morphological features) to reduce the multivariate dataset into two PC scores, with PC1 and PC2 capturing 61.2% and 23.2% of the morphological data variance, respectively. From the factor loading plot, size‐based parameters were more highly correlated with PC1 while shape‐based parameters were more highly correlated with PC2 (Supporting Information Fig. S3). As a result, PC1 (Size) and PC2 (Shape) represent overall morphological scores that relate to MSC aggregate size and shape, respectively. Therefore, we assessed whether these scores could be used to determine 3D morphological dynamics of different cell‐lines and passages following chondrogenic stimulation (Fig. 4A).

Figure 4

(A): Multipotent stromal cell (MSC) aggregates exhibit different patterns in aggregate size changes, as represented by the PC1 score derived from the PCA performed on the entire morphological dataset, during the culture period depending on both the cell line and passage. The PC1 (Size) scores of the aggregates were compared between groups and shown at Day 1, Day 4, Day 7, and Day 21. * Indicates a statistically significant difference between early passage and late passage samples (****p < .0001, two‐way ANOVA). Data are presented as the mean ± SD values of the biological replicates included per group (n = 10–15) for the experiments. Bright‐field images depicting the MSC aggregates from each group at Day 7 and Day 21 and at both passages for each time point are shown in (B). Scale bar = 500 μm.

Several MSC lines at LP demonstrated incomplete aggregation into spheroids as indicated by high Day 1 PC1 (Size) scores (Fig. 4A). While most MSC aggregates decreased in size by Day 4, several MSC aggregates exhibited size increases after this time point where EP aggregates continuously maintained larger sizes when compared to LP aggregates. Comparing PC1 (Size) scores beginning at Day 1, differences in aggregate size due to donor variations became readily apparent, where different morphological trends also emerged depending on passage. By Day 21, the same five larger EP groups (RB9, RB12, RB14, RB16, and PCBM1641) had exhibited continued increases in size over time whereas the initially smaller EP groups (167696 and PCBM1632) had further compacted during culture. Although most LP MSC lines continued to decrease in size until Day 21, one LP group (RB9) grew significantly larger and was comparable in size to aggregates from other EP MSC lines. There were no discernible trends in PC2 (Shape) dynamics with respect to cell line or passage except for RB12 (Supporting Information Fig. S4). Indeed, shape differences between MSC lines were not apparent until Day 14. Representative aggregates from all MSC lines at Days 7 and 21 are shown in Figure 4B.

We next investigated cell line variations in the morphological dynamics of MSC aggregates by mapping the raw PC1 (Size) and PC2 (Shape) scores over the entire culture duration for each MSC line (Fig. 5A). We immediately observed from their morphological trajectories (Fig. 5B) that RB14, PCBM1641, and to an extent 8F3560 MSCs shared very similar changes in aggregate morphologies over time. EP and LP aggregates from these cell lines were all similarly sized at the beginning of culture but then significantly diverged in size after Day 4. Unlike RB14 and PCBM1641, however, EP 8F3560 aggregates did not surpass their original size at D1 over time. RB16 and RB12 shared similar morphological dynamics in that both their LP aggregates were initially larger and differently shaped (due to incomplete aggregation) compared to EP aggregates. Over time, LP RB16 and RB12 aggregates became more spherical but decreased significantly in size when compared to EP aggregates, which grew larger than their Day 1 counterparts. 167696 and PCBM1632 aggregates showed very comparable morphological profiles, where both EP and LP aggregates diminished significantly in size after Day 7. The RB9 cell line was unique in that both EP and LP aggregates regained size over time after an initial decrease.

Figure 5

(A): The dynamics of the multipotent stromal cell aggregate size and shape over the culture period are shown for each cell line, where (B) the mean trajectory in morphology PC scores for each cell line are also tracked over time for both early and late passage groups.

Morphological Phenotypes Emerge upon Chondrogenic Stimulation That Predict Overall Chondrogenic Capacity

To identify emergent MSC aggregate morphological phenotypes predictive of chondrogenic capacity, correlations were made between MSC morphological scores at each time point and their measures of chondrogenic differentiation (Fig. 6). Raw scores were rescaled via normalization prior to plotting to remove negative values for improved clarity. We observed that PC1 (Size) scores became strongly correlated with the SynthAct and ChondroInd scores as early as Day 7 (Fig. 6A), where R values exceeded 0.8 after 1 week of culture for both SynthAct and ChondroInd measures. However, variability in MSC aggregate size‐based morphological features did not correlate with ChondroGene scores at any time point. In contrast to PC1 (Size) scores, the PC2 (Shape) scores achieved statistically significant correlations with SynthAct and ChondroInd scores much later at Day 14 and beyond (Fig. 6B). Furthermore, correlations between PC2 (Shape) scores and the chondrogenic outcomes were not as strong as those for PC1 (Size) scores (lower R values). As was the case for PC1 (Size) scores, PC2 (Shape) scores did not correlate with ChondroGene scores at any time point.

Figure 6

Normalized PC1 scores representing MSC synthetic activity, chondrogenic gene expression, and overall chondrogenic index were correlated with (A) PC1 (Size) and (B) PC2 (Shape) morphology scores at various time points measured during the 21‐day aggregate culture. The correlation coefficient R measures the strength of the linear relationship between correlated variables, where * indicates a statistically significant correlation (*p < .05, **p < .005, ***p < .0005, and ****p < .0001, Pearson product–moment correlation).

Individual morphological features measuring both MSC aggregate size and shape were also independently correlated with the multiple chondrogenic metrics to elucidate which morphological feature may be most indicative of MSC chondrogenic potential. While most of the size‐based morphological features became significantly correlated with SynthAct and the ChondroInd by Day 7, the three measures of aggregate radii (MaxRadius, MeanRadius, and MedianRadius) were already strongly correlated by Day 4 (Supporting Information Table S5). Correlations between some shape‐based features and the chondrogenic metrics including the ChondroGene scores were also observed, namely at Day 11. Of the multiple shape‐based morphological features measured, only the form factor (i.e., circularity) maintained significant correlations with the SynthAct, ChondroGene, and ChondroInd scores at Day 11, where less circular spheroids tended to have higher chondrogenic potential (Supporting Information Table S6).

Gene Expression of Undifferentiated MSCs Correlate with Measures of Chondrogenic Differentiation

Global gene microarray analysis (62,976 probes) was performed on the MSC lines without chondrogenic stimulation to identify inherent differences in gene expression that may correlate with their chondrogenic capacity. 750 probes were found to be significantly different when comparing EP versus LP MSCs (Supporting Information Table S7), where only 32 genes exhibited absolute fold changes greater than 1.5 (Supporting Information Fig. S5A). When MSC lines were grouped and compared according to higher (RB9EP/LP, RB12EP, RB14EP, RB16EP/LP, PCBM1641EP, 8F3560EP) or lower (RB12LP, RB14LP, PCBM1641LP, 8F3560LP, PCBM1632EP/LP, 167696EP/LP) chondrogenic potential (based on ChondroInd; cutoff value set by RB16LP), 66 probes were identified (Supporting Information Table S8) to be either significantly upregulated or downregulated (Supporting Information Fig. S5B). When correlating the microarray data of the undifferentiated MSC lines with the corresponding ChondroInd scores following chondrogenic stimulation, a total of 107 probes were observed to be statistically significant (Supporting Information Fig. S5C; Supporting Information Table S9). A preliminary pathway analysis was also performed to identify associated biofunctions (Supporting Information Table S10). A hierarchical clustering heat map was generated to illustrate differences in gene expression, where MSC lines are ordered from high to low chondrogenic potential based on ChondroInd scores (Fig. 7).

Figure 7

Hierarchical clustering heat map illustrates gradual increases and decreases in the expression of genes in undifferentiated multipotent stromal cells (MSCs) found to be significantly correlated with their ChondroInd scores at Day 21. MSC lines are ordered from highest to lowest chondrogenic potential based on ChondroInd scores.

Discussion

MSC heterogeneity represents a major obstacle to establishing critical quality attributes that can effectively predict the functional capacity of MSCs relevant to therapeutic applications such as chondrogenesis. MSC functional heterogeneity (differences in therapeutic function 17) across donors and even within a single population can derive from both intrinsic and extrinsic factors 19, the latter of which also include various manufacturing parameters (i.e., different cell isolation/expansion protocols, general media composition, and other operator effects). New innovative approaches have begun to correlate the morphology of MSC subpopulations within larger populations to their differentiation (adipogenic and osteogenic) 20, 21 and immunosuppressive 23 potential. While several culture methods and assays for characterizing MSC chondrogenic differentiation do exist 27, 28, there still lacks well‐defined approaches that seek to identify cell critical quality attributes predictive of MSC chondrogenic potential. Toward this goal, we characterized the morphology of MSC aggregates under chondrogenic stimulation from multiple donors and demonstrated that MSC aggregates in 3D culture exhibited unique morphological dynamics that were both passage and cell line dependent. We showed that MSC chondrogenic differentiation becomes highly correlated with size‐based features of formed aggregates much earlier during 3D culture versus shape‐based features, suggesting the ability to use quantitative aggregate morphology as a simple tool to estimate the chondrogenic capacity of an MSC line.

MSC heterogeneity is still often overlooked as many research groups continue to generalize the regenerative potential of selected cell lines to all MSCs 18. However, research efforts have been more recently focused on the characterization of MSC heterogeneity with the hope of developing adaptable tools using certain quality attributes for predicting cell behavior. Comparing eight different MSC lines, it was shown that undifferentiated MSCs isolated from different donors exhibited varying colony‐forming unit (CFU) capacities 29, indicating varying degrees of stemness. Yet, stemness levels as measured by both %CFU and flow cytometry (for MSC surface markers) did not reflect their ability to undergo adipogenic differentiation, highlighting inadequacies in the use of cell surface markers alone as a biomarker for MSCs (a common practice). Recent efforts characterizing the gene expression profile of proliferating MSCs via microarray analysis also showed that gene expression changes associated with MSC aging during passage failed to correspond with any of the markers established by the International Society for Cellular Therapy 30, 31. However, the identification of new epigenetic 32 and genetic markers 33 correlated with cell function could provide quality tools for assaying MSC potency.

One significant advantage of our investigation involves the use of multiple MSC lines that were also previously evaluated for differences in their adipogenic 29, osteogenic 20, proliferative 30, and immunosuppressive 23 capacities. It was noted previously that while the proliferative patterns (time to confluency, cell size, and morphology) of MSCs were similar at earlier passages 29, 33, 34, they began to diverge at later passages as the MSCs derived from various donor source exhibited differential gene expression of markers indicative of senescence 30. However, MSCs manufactured from different donor sources exhibited differences in genetic stability and varied in adipogenic, osteogenic, and immunosuppressive potential regardless of passage. Unlike the referred studies, we evaluated MSCs in 3D culture where variations in the observed readouts may potentially reflect dissimilarities in cellular processes originating from mesenchymal condensation. Previous studies have shown that MSCs cultured as 3D spheroids versus 2D monolayers improved overall survival 35 and increased their functional tropism 36. Although not directly tested in our case, it is a possibility that the transition from 2D monolayers to 3D spheroids did not elicit the same universal response in all MSCs. We noted that under similar chondrogenic conditions, six out of eight MSC lines produced significant amounts of GAG at EP whereas only two out of eight MSC lines displayed any GAG production at LP. Interestingly, histological analysis performed on MSC groups with high GAG production (at both passages) revealed not only larger aggregates overall but also extracellular matrix distributions similar to native articular cartilage; GAGs were deposited more densely at the aggregate core while collagens were more heavily aligned along the spheroid periphery. Such was not the case for MSC groups presenting minimal extracellular matrix production where the GAG deposition overlapped with that of collagen. This perhaps indicates an impaired ability of these MSCs to rearrange their aggregate packing density for metabolic adaptation 35, which merits further investigation.

As it was previously shown that the mRNA expression of the canonical gene aggrecan failed to correlate with MSC matrix production at the single‐cell level 37, we developed an approach using PCA to segregate metrics of chondrogenic differentiation into individual composite scores representative of chondrogenic synthetic activity (SynthAct) and gene expression (ChondroGene). We evaluated the correlation between functional chondrogenic differentiation and MSC aggregate morphology, a phenotypic readout that can be readily measured over time using nondestructive imaging techniques. Although some studies have peripherally explored the association between MSC spheroid morphology and differentiation states 9, 38, few attempts have since been made at establishing quantitative relationships. While our study used 2D projections of multicellular spheroids for high‐dimensional morphological analysis, the evaluation of volumetric spheroid morphology 39 could potentially provide further biological insight. Compared to larger MSC aggregates, smaller MSC aggregates may consist of MSC populations containing more subpopulations of highly contractile cells that led to overall greater aggregate compaction, which requires further study. Indeed, histological examination of smaller aggregates generally revealed more collagen deposition at the aggregate cores, indicating perhaps the increased cortical tension and cell‐matrix interactions characteristic of more contractile cells 40. These results suggest the strong possibility of adapting this phenotypic readout for quickly estimating functional chondrogenic matrix accumulation in lieu of other more costly and arduous processes. Delay in the time to correlation, as well as lower correlation strengths, between aggregate shape and measures of chondrogenic functional capacity point to a disconnect between cell spheroid size and spheroid shape that may be of biological interest. Additionally, the lack of any correlation between ChondroGene and aggregate morphology at any time point further highlights the need to establish a phenotypic basis in chondrogenic gene expression.

Despite significant discrepancies between the chondrogenic gene expression data and phenotypic outcomes (i.e., synthetic activity and histology), we further explored potential genetic signatures of MSC lines prechondrogenic induction that may be indicative of chondrogenic potential. Of the 107 probes (out of 62,976) found to be significantly correlated with the ChondroInd of the MSC cohort, a number of prominent genes were identified in the literature to be implicated in several biological functions related to cell–cell signaling, tissue and organ morphology, and musculoskeletal development. Given the varying dynamics observed in MSC aggregate morphology, which embodies a form of organoid culture, it is expected that such biological functions were found to be associated with chondrogenic potential. Genes found to be positively (JAM2, COL12A1, LUM, and FGF7) and negatively (DNER, DKK1, ADAMTS14, and SHB) correlated with ChondroInd suggest that MSCs of higher chondrogenic capacity may be inherently more migratory, more efficient at matrix reorganization, and more metabolically adaptable 41, 42, 43, 44, 45, 46, which prompted further investigation. Comparing EP and LP MSCs, genes revealed to be differentially expressed largely corroborated previous data comparing the effects of passage in other MSC lines 30, 33. While we provide an exploratory analysis of the relationship between MSC chondrogenic phenotype and their source gene expression, a more comprehensive genomic analysis could reveal gene signatures predictive of MSC chondrogenic potential, which is the topic of our ongoing work.

Conclusion

In summary, we demonstrated that multiple MSC lines exhibited both cell line‐ and passage‐dependent aggregate morphologies during spheroid culture that correlated strongly with chondrogenic capacity. Using temporal high‐dimensional morphological analysis of MSC aggregates, we created composite scores for morphological data based on size and shape morphological features, which were used to show that aggregates formed by MSCs from different donors diverged more in size than shape between cell lines following chondrogenic stimulation. Applying similar methodology for the analysis of chondrogenic outcomes, we further derived quantitative scores to represent the individual as well as combined matrix accumulation and gene expression aspects of chondrogenic differentiation. We showed via correlation analysis that functional matrix accumulation (SynthAct) but not chondrogenic gene expression (ChondroGene) correlated strongly with aggregate morphology as early as Day 7, highlighting a discrepancy between chondrogenic phenotype and gene expression. Furthermore, MSC aggregate size correlated with chondrogenic synthetic activity much earlier than MSC aggregate shape, which became significantly correlated by Day 14. Correlation of these phenotypic data with base gene expression of unstimulated MSCs featured several biological pathways and markers that provide key starting points for future research endeavors toward assaying MSCs. Finally, by developing a simple nondestructive approach demonstrating a strong correlation between early MSC aggregate morphology and chondrogenic potential, we provide a well‐defined tool for the early estimation of MSC chondrogenic differentiation capacity. Indeed, the simple measurement of MSC aggregate morphology represents a useful resource for the improvement of quality for manufactured cell products.

Author Contributions

J.L. and K.E.S.: conceived and designed the research; J.L., I.H.B., and R.A.M.: performed the experiments and data analysis; J.L., I.H.B., R.A.M., S.R.B., R.K.P., and K.E.S.: contributed intellectual input and critically reviewed and wrote the manuscript.

Disclosure of Potential Conflicts of Interest

The authors indicated no potential conflicts of interest.

Appendix S1: Supporting Information

Click here for additional data file.

Appendix S2: Supporting Information

Click here for additional data file.

Supplemental Table S1: Morphological size and shape features measured by the CellProfiler image analysis software. Formal definitions for each feature were acquired from the CellProfiler manual (http://cellprofiler.org/manuals.shtml)) under the “MeasureObjectSizeShape” section.

Click here for additional data file.

Supplemental Table S2: CellProfiler algorithm (pipeline) used to quantify high‐dimensional morphological size and shape features of MSC aggregates

Click here for additional data file.

Supplemental Table S3: Primer sequences used for RT‐PCR

Click here for additional data file.

Supplemental Table S4: CellProfiler algorithm (pipeline) used to quantify histology images of MSC aggregates.

Click here for additional data file.

Supplemental Table S5: Correlations between MSC aggregate size‐based morphological features and composite metrics of chondrogenic differentiation. The correlation coefficient measures the strength of the linear relationship between correlated variables, where * indicates a statistically significant correlation (*p < 0.05, Pearson product–moment correlation).

Click here for additional data file.

Supplemental Table S6: Correlations between MSC aggregate shape‐based morphological features and composite metrics of chondrogenic differentiation. The correlation coefficient measures the strength of the linear relationship between correlated variables, where * indicates a statistically significant correlation (*p < 0.05, Pearson product–moment correlation).

Click here for additional data file.

Supplemental Table S7: A list of statistically significant genes with adjusted P values and indicated fold changes comparing early versus late passage MSC lines

Click here for additional data file.

Supplemental Table S8: A list of statistically significant genes with adjusted P values and indicated fold changes comparing undifferentiated MSC lines identified to have high versus low chondrogenic potential

Click here for additional data file.

Supplemental Table S9: A list of genes exhibiting statistically significant correlations and slopes with the Chondrogenic Index metric

Click here for additional data file.

Supplemental Table S10:

Click here for additional data file.

Supplemental Figure S1: Quantification of (A) GAG (Alcian Blue) and (B) Collagen (Picrosirius Red) histology show donor‐ and passage‐dependent differences. * indicates a statistically significant difference between early passage and late passage samples (****p < 0.0001, two‐way ANOVA). Data are presented as the mean ±SD values of the biological replicates included per group (n = 3) for these analyses. Quantified measures of GAG and Collagen histology were correlated with biochemically detected matrix deposition (GAG/DNA). The correlation coefficient R measures the strength of the linear relationship between correlated variables, where * indicates a statistically significant correlation (****p < 0.0001, Pearson product–moment correlation).

Click here for additional data file.

Supplemental Figure S2: (A) PCA was performed using supervised datasets comprising all metrics of chondrogenic gene expression (type X collagen, type I collagen, type IIα1 collagen, aggrecan, and type IIα1/type I collagen fold ratio) as illustrated by the colored box. The ▪ and ▲ represent an early passage or late passage sample, respectively. The normalized PC1 scores representing MSC SynthAct were correlated with (B) the normalized ChondroGene scores (supervised for type IIα1 collagen, aggrecan, and type IIα1/type I collagen fold ratio only), where correlation coefficient R measures the strength of the linear relationship between correlated variables. (C) ChondroGene scores were also correlated with quantified measures of GAG and Collagen histology. * indicates a statistically significant correlation (*p < 0.05, Pearson product–moment correlation).

Click here for additional data file.

Supplemental Figure S3: The factor loading plot of the performed PCA is shown, where “Size” PC1 captures most of the variance in the morphological dataset related to size‐based parameters and “Shape” PC2 captures most of the variance in the morphological dataset related to shape‐based parameters. The associated table of correlations for the PCA performed on the unsupervised morphological dataset (percentages indicate amount of variance captured by each PC) is also shown.

Click here for additional data file.

Supplemental Figure S4: Normalized “Shape” PC2 scores derived from the PCA performed on the entire morphological dataset were compared between groups over time. Comparisons are also shown at Day 1, Day 4, Day 14, and Day 21. * indicates a statistically significant difference between early passage and late passage samples (*p < 0.05, and ****p < 0.0001, two‐way ANOVA). Data are presented as the mean ±SD values of the biological replicates included per group (n = 12–15) for the experiments.

Click here for additional data file.

Supplemental Figure S5: Hierarchical clustering heat map comparing the differential gene expression between (A) early versus late passage MSCs and between (B) MSC lines identified to have higher versus lower chondrogenic potential. (C) Using a regression model, a slope and correlation coefficient were calculated per probe for the ChondroInd scores. A cutoff slope was determined per probe using the percent coefficient of variation and technical variability per probe. Probes whose calculated slope from the regression model was greater than the cutoff slope were considered to have passed this filter and subjected to further statistical analysis. Probes that passed this filter typically had the largest calculated slopes and correlation coefficients closest to 1 or − 1.

Click here for additional data file.