Cellular features of localized microenvironments in human meniscal degeneration: a single-cell transcriptomic study

Abstract
Editor's evaluation
Introduction
Methods
Results
Discussion
Data availability
References
Article and author information
Metrics

Abstract

Background:

Musculoskeletal tissue degeneration impairs the life quality and function of many people. Meniscus degeneration is a major origin of knee osteoarthritis and a common threat to athletic ability, but its cellular mechanism remains elusive.

Methods:

We built a cell atlas of 12 healthy or degenerated human meniscus samples from the inner and outer meniscal zones of 8 patients using scRNA-seq to investigate meniscal microenvironment homeostasis and its changes in the degeneration process and verified findings with immunofluorescent imaging.

Results:

We identified and localized cell types in inner and outer meniscus and found new chondrocyte subtypes associated with degeneration. The observations suggested understandings on how cellular compositions, functions, and interactions participated in degeneration, and on the possible loop-like interactions among extracellular matrix disassembly, angiogenesis, and inflammation in driving the degeneration.

Conclusions:

The study provided a rich resource reflecting variations in the meniscal microenvironment during degeneration and suggested new cell subtypes as potential therapeutic targets. The hypothesized mechanism could also be a general model for other joint degenerations.

Funding:

The National Natural Science Foundation of China (81972123, 82172508, 62050178, 61721003), the National Key Research and Development Program of China (2021YFF1200901), Fundamental Research Funds for the Central Universities (2015SCU04A40); The Innovative Spark Project of Sichuan University (2018SCUH0034); Sichuan Science and Technology Program (2020YFH0075); Chengdu Science and Technology Bureau Project (2019-YF05-00090-SN); 1.3.5 Project for Disciplines of Excellence of West China Hospital Sichuan University (ZYJC21030, ZY2017301); 1.3.5 Project for Disciplines of Excellence – Clinical Research Incubation Project, West China Hospital, Sichuan University (2019HXFH039).

Editor's evaluation

This paper will be of interest to researchers studying meniscus homeostasis and knee osteoarthritis. It uncovers distinct subtypes of cell populations in the inner and outer parts of the human meniscus using single-cell RNA sequencing. In particular, this work further identifies how alterations in meniscal cell populations may contribute to inflammation and osteoarthritis and thus serves as a resource paper for the field.

https://doi.org/10.7554/eLife.79585.sa0

Introduction

The knee is the largest and most complicated hinge joint in the human body. The meniscus is an important component of the knee joint with essential roles in the load-bearing, shock absorption, nutrition, and lubrication of the knee joint articular cartilage and is very susceptible to injury and has a limited reparative potential (Murphy et al., 2019; Newman et al., 1989; Rai and McNulty, 2017). The meniscus degeneration is a slow-developing disease of middle-aged or older people. Horizontal cleavage of the meniscus usually happens in the disease, and the lesion can be identified through knee MRI. Patients typically have no clear history of acute knee injury (Beaufils et al., 2017). Meniscal degeneration is one important risk factor for osteoarthritis (OA) and joint dysfunction, causing huge social and economic burdens (Englund et al., 2012). It has been found through clinical observations that nearly half of patients with meniscal degeneration will eventually develop OA over time (Lohmander et al., 2007). Studies have shown that inflammatory responses and biomineralization affect the meniscus cellular microenvironment and contribute to disease occurrence (Sun et al., 2010; Goldring and Goldring, 2016). However, the changes in the cell microenvironment in degenerated meniscus remain largely unclear. An adequate understanding of the changes in the degenerative meniscus is critical for preventing meniscal injury in young and middle-aged patients and for relieving symptoms in elderly patients with knee OA.

The meniscus can be anatomically divided into the inner region (the white–white zone) and the outer region (the red–red zone). Some researchers also use the ‘red–white’ zone to denote the in-between transition region. The inner and outer menisci differ in neurovascular distribution, structural composition, and recovery ability from injury. There are blood vessels and nerves in the outer area, while few blood vessels and nerves are observed in the inner part. The meniscal extracellular matrix (ECM) is the physical foundation of menisci’s biological roles, and alterations in the ECM may lead to meniscus degeneration and dysfunction. Collagens and proteoglycan (PG) are ECM’s two most important components. The outer meniscus is dominated by type I collagen, while the inner meniscus is dominated by type II collagen. The inner part of the meniscus has a relatively higher percentage of PGs compared with the outer. Aging and mechanical injuries may lead to meniscus degeneration. Different anatomical regions have varied recovery capabilities: the degeneration in the outer meniscus is more likely to heal, while the degeneration in the inner part tends to be irreversible (Makris et al., 2011; Danzig et al., 1983). The cellular and molecular basis behind anatomical regions remains to be explored. It is unclear what cell type contributes to the inner and outer variations and what cell types mediate the degeneration.

Building a detailed meniscal cell landscape is essential to understanding meniscus characteristics. Early profiling of meniscal cell heterogeneities dates back to the work of Ghadially et al., 1983. They studied injured and uninjured human menisci with the electron microscope and stated that chondrocytes, a few fibroblasts, myofibroblasts, and intermediate state cells between chondrocytes and fibroblast existed in menisci. Scotti et al. also discussed that the cells in the superficial zone were fusiform, while cells laid deeper were polygonal (Scotti et al., 2013). These classifications were mainly based on the shape and the surrounding matrix content. The recent development of single-cell omics has facilitated researchers’ understanding of various cell types across multiple cartilage tissues, including knee joint cartilage (Ji et al., 2019), intervertebral disc (Gan et al., 2021), meniscus (Sun et al., 2020), etc. These researches initiated a rudiment of the chondrocyte cell types and states, covering parts of the heterogeneous cell subpopulations in the meniscus. Sun et al., 2020 sketched the outline of the meniscal cell populations following the cell type definition convention of a knee articular cartilage study (Ji et al., 2019). They reported cartilage progenitor cells (CPC), regulatory chondrocytes (RegC), prehypertrophic chondrocytes (PreHTC), hypertrophic chondrocytes (HTC), fibrochondrocytes (FC), fibrochondrocytes progenitors (FCP), proliferating fibrochondrocytes (ProFC), degenerated progenitor (DegP), endothelial cells, and other immune cells in the meniscus, and showed links between meniscus progenitors and the progression of meniscal degeneration (Sun et al., 2020).

We aimed to build a refined atlas complementary to the existing data by systematically revealing cell heterogeneities in inner and outer meniscal zones under different health states. We created a single-cell transcriptomic atlas of 45,744 cells from 12 healthy/degenerated meniscus samples of 8 patients and provided an online cell browser at http://meni.singlecell.info:3000/ to support flexible explorations of the data. We performed multiple rounds of cell-level quality control steps to remove and adjust noisy data such as empty droplets, doublets, dissociation-affected cells, and ambient RNA-polluted cells. Among the current single-cell sequencing studies involving articular cartilage and meniscus, we have a larger number of sequencing samples and cells. Our dataset is the first single-cell sequencing study data characterizing menisci’s inner/outer zonation differences. We built a hierarchical cell type classification framework and identified five chondrocyte subtypes and two pericyte-like cell subtypes. Functional enrichment analyses suggested these subtypes may have specialized ECM construction and remodeling duties. We profiled various immune cell types and observed damage-induced inflammatory responses in degenerated samples. Taken together, we inferred a hypothetical model of ECM disassembly, inflammation, and angiogenesis in the meniscal tissues. The model may also fit degenerative situations in other cartilage tissues, such as articular cartilage and intervertebral disc, because they contain similar cell types and ECM contents (Chen et al., 2017b). Our discoveries may suggest new clues for treating various types of joint degeneration in the human body.

Methods

Human meniscus cell sample preparation

Four degenerative meniscus specimens were obtained from patients with severe OA who underwent total knee arthroplasty, and four normal meniscus specimens were obtained from patients with bone tumors or severe trauma who underwent amputation. The meniscus tissues were removed from the patient’s knee joints and divided into the inner (white–white zone) and outer (red–red zone) areas. The degeneration group is comprised of patients A, B, C, and D; and the normal group is comprised of patients E, F, G, and H. Patients A, B, E, and F contributed paired inner–outer samples. Patients C and H contributed unpaired inner samples, and patients D and G contributed unpaired outer samples because the other pairs had obtained in the surgery had poor sample states. Detailed designs of experiments can be found in Supplementary file 1.

After carefully dividing the meniscus tissue into the inner and outer parts, the menisci were cut into small pieces. We used 0.25% Trypsin (Gibco 25200072) to digest the pieces for 0.5 hr under 37℃. Next, we used 2 mg/ml collagenase IV to dissociate the pieces for 4–6 hr under 37℃ and finally obtained dissociated cells. The dissociated cells were resuspended at a concentration of 250–1200 cells/μl and with viability between 67–87% for microfluidics (Chromium Single-Cell Controller, 10× Genomics). According to the manufacturer’s instructions, cells were loaded into the chip and run using the Chromium Single Cell 3′ Reagent Kit v2 (10× Genomics). In brief, Gel beads with GemCode barcodes and primers and individual cells were encapsulated in oil droplets. Next, within each oil droplet, mRNA were released from lysed cells, barcoded, reverse transcribed to cDNA, and sequenced. Sequencing-ready single-cell transcriptome libraries were mapped to the human reference genome GRCh38-3.0.0 and quantified by CellRanger 3.1.0.

Multiplex immunofluorescence (OPAL) staining

Human meniscus normal and degenerated samples were fixed in 4% paraformaldehyde and embedded in paraffin. We sliced the embedded paraffin samples in series, and each slice was 4 µm thick. Firstly, water-bath heating was used for antigen retrieval. Then various cell marker primary antibodies were incubated with the paraffin slide of meniscus tissue to conduct the continuous staining with the Opal Polaris Multiple-Color Manual IHC Kit (NEL861001KT). We used different primary antibodies to simultaneously label Ch.1 (SERPINA1), Ch.2 (MMP14), Ch.4 (CYP1B1), and PCL (ACTA2) on the same tissue slide. The automated staining system (BOND-RX, Leica Microsystems, Vista, CA) was used to perform the chromogen-based multiplex immunohistochemistry labeling. SERPINA1 (ab207303), MMP14 (ab51074), CYP1B1 (ab33586), and ACTA2 (ab5694) were all purchased from Abcam and diluted at a concentration of 1:100. The Opal Polaris dyes were used to pair with these antibodies containing fluorophores for tyramide signal amplification to enhance sensitivity. The sequence of labeling for detecting each marker was optimized: CYP1B1 (Opal 570), MMP14 (Opal 620), SERPINA1 (Opal 690), ACTA2 (Opal 780), and DAPI (4',6-diamidino-2-phenylindole). Including an autofluorescence section, the staining process is the same as above, but no primary antibody is added. Multiplex analysis was operated to analyze the results of simultaneously stained tissue slides. In addition, common immunofluorescence staining was performed using CDON (ab227056), CD31 (ab9498), and CD45 (ab40763) labeled Ch.3, endothelial cells, and immune cells. All experiments have three biological replicates.

Single-cell data quality control steps

We adjusted the ambient RNA expression with SoupX (Young and Behjati, 2020) v1.2.1 and generated read counts free of background noises for all downstream analyses. We used scanpy (Wolf et al., 2018) and Seurat (Stuart et al., 2019; Butler et al., 2018) to conduct the basic filtering out outliers cells with high n_genes metrics (the number of observed genes per cell), high n_counts metrics (the number of UMI per cell), and high percent_mito metrics (the mitochondrial gene fraction). Next, we assigned draft cell type identities to the cells using SingleR (Aran et al., 2019). Note that the SingleR-generated cell type labels were only used to perform quality control analysis to avoid removing biologically meaningful cells with poor data qualities. We calculated doublet scores using DoubletFinder (McGinnis et al., 2019) removed the cells with high doublet scores and outlier clusters with high doublet scores. We observed that DoubletFinder tends to identify small clusters as doublets incorrectly, so we assert these small populations as singlets if they had well-defined identities given by SingleR or appear clear biological functions such as cell cycles. To remove the confounding factor introduced by the dissociation-induced gene expression (van den Brink et al., 2017; Denisenko et al., 2020), we decided to remove the heavily affected cells by the dissociation steps. We calculated dissociation scores using methods described in van den Brink et al., 2017.

Raw single-cell RNA-seq datasets contain numerous low-quality droplets. Hence, we should ensure the majority of the data items we analyzed correspond to viable cells. It is known that one simple quality control filtering on the aforementioned metrics (n_genes, n_counts, percent_mt) can't remove all of them. We then conducted comprehensive quality control steps that can be summarized as (1) top-down refinements and (2) divide and conquer. First, we divide cells into several rough populations with distinct biological identities, for example, chondrocytes, immune cells, endothelial cells, etc. We zoomed into these populations one by one and performed subclustering within them. The within-population subclustering often works like a centrifuge – the low-quality cells mixed in the large rough population usually form tiny outlier groups when subclustering is performed on a single population. When a tiny group is mainly comprised of cells with high doublet scores, cells with high mitochondrial gene fractions, cells with high dissociation scores, we consider removing it. After manually inspecting and dropping these tiny outlier groups, we finally get much cleaner datasets.

Single-cell data clustering steps

To eliminate the genetic background variations and technical noises of the single-cell RNA-seq data, we used Harmony (Korsunsky et al., 2019) to perform integrative clustering across different samples. Major clusters corresponded to chondral cells, ACTA2⁺ cells, endothelial cells, and immune cells were identified across samples and validated with cluster correlation analysis (Figure 1—figure supplement 1C and Figure 2—figure supplement 1A).

To understand large cell populations with higher resolutions, we subset the chondral and ACTA2⁺ majors clusters, performed reclustering on harmony-derived embeddings, and obtained multiple fine clusters. To alleviate the scRNA-seq data noise and promote data consistency across samples (Hua and Zhang, 2019), we partitioned original cell data points into homogenous meta-cells with MetaCell (Baran et al., 2019). We then used HGC (Zou et al., 2021) to build hierarchical relationships based on the meta-cells.

Differentially expressed genes analysis

Differentially expressed genes (DEGs) were identified by the ‘FindMarkers’ function and the ‘FindAllMarkers’ function provided by the Seurat package. The DEGs of each chondrocyte cluster were identified by comparing the cluster with all the other clusters. A full list of DEGs in shown in Figure 2 is given in Supplementary file 2. The health group (degeneration/normal) markers were calculated in an overall comparison way and a per-cluster comparison way. We first extracted the total chondrocytes (from Ch.1 to Ch.5) and compared them between the degenerated and normal groups. We made the same comparison to the total PCL populations. The health group markers are reported in Figure 3—figure supplement 5A.

Gene set enrichment analysis and gene set variation analysis

Enrichment analyses were performed with top fold-change DEGs and a BH-adjusted p-value threshold of 0.05, using R package clusterProfiler (Wu et al., 2021). Cluster-wise enriched terms (Ch.1–5 and PCL.1/2, Figure 2—figure supplement 2) were obtained using cluster marker DEGs (one cluster vs. the others). Health state-specific terms were obtained using status DEGs (all chondrocytes in normal samples vs. all those in the degenerated) and displayed in Figure 3—figure supplement 5B, C.

We assigned gene set activity scores to individual cells using gene set variation analysis (GSVA) (Hänzelmann et al., 2013), with gene sets obtained from MSigDB (Subramanian et al., 2005; Liberzon et al., 2015). The GSVA-derived gene set scores were visualized in Figures 2G and 3C–F. The GSVA analyses used the following MSigDB-curated gene ontology terms: GO0050900 (leukocyte migration), GO1905563 (negative regulation of vascular endothelial cell proliferation), GO1904018 (positive regulation of vasculature development), GO0071772 (response to BMP), GO0030199 (collagen fibril organization), GO0070555 (response to interleukin 1), GO0071559 (response to transforming growth factor beta), GO:0034612 (response to tumor necrosis factor), and GO0030239 (myofibril assembly).

Gene regulation network inference

We used the python package pySCENIC (Van de Sande et al., 2020; Aibar et al., 2017) to infer gene regulatory networks (GRN) from single-cell expressional profiles. The workflow includes three substeps: identifying coexpression modules, trimming modules with prior knowledge, and calculating activating scores.

The workflow builds coexpression modules by choosing TFs that predict target gene expression well. First, the algorithm builds gradient boosting machine models predicting each target gene’s expression from TF expressions and keeps the models’ feature weights as a measure of TF–target regulatory scores. Second, the algorithm creates modules by adding/filtering weighted links between TFs and their targets. To reduce the stochasticity of the results, the algorithm independently uses six weight filtering rules to build modules. The six rules are: (1) keep TF–target pairs with weights > the 75% percentile of the weights; (2) keep TF–target pairs with weights > the 90% percentile of the weights; (3) for each TF, keep the top 50 TF–target pairs with the highest weights; (4) for each target, keep top 5 TF–target pairs with highest weights; (5) for each target, keep top 10 TF–target pairs with highest weights; (6) for each target, keep top 50 TF–target pairs with highest weights. After this step, many TF–target modules are created. The algorithm only keeps TF–target links that have positive regulation relationships (the Pearson correlation of TF and target expressions >0.03) and drops modules with <20 genes inside (TF itself included). Finally, we get a list of primary coexpression modules, each containing a TF and its positively regulated targets.

The module trimming process combines the regulatory motif information with the coexpression analyses. Given the coexpression-derived modules, we then keep target genes that are overlapped with TF’s regulatory motif region. We adopt the 10 kb up-/downstream regions of a target gene’s TSS as the gene’s distal enhancer occurrence region and the 500 bp upstream/100 bp downstream region as the gene’s proximal promoter occurrence region. After trimming, more reliable sets of modules are generated.

Finally, we quantified and binarized the activation levels of the TF–target modules with AUCell (Aibar et al., 2017) and visualized them with ComplexHeatmap (Gu et al., 2016) and cytoscape (Shannon et al., 2003).

Cell crosstalk analysis

Given the averaged normalized expression value of a ligand gene in a sender cluster $i$ as $E_{L, i}$ , and that of a receptor gene in receiver cluster $j$ as $E_{R, j}$ , we defined the crosstalk between cluster $i$ and cluster $j$ via the ligand–receptor pair L–R is their product $E_{L, i} E_{R, j}$ . We enumerated all ligand and receptor genes provided by CellTalkDB (Shao et al., 2021) and calculated crosstalk products for sender–receiver cluster pairs. We empirically evaluated the statistical significance of this product by counting the more extreme null products after shuffling the cluster labels. The overall crosstalk level between the sender cluster $i$ and the receiver cluster $j$ is defined as $\sum_{L} \sum_{R} E_{L, i} E_{R, j}$ if $E_{L, i} E_{R, j} > θ$ , and $θ = 1.5$ in our cases.

We compared the crosstalk variations between the degenerated and the normal group by subtracting the products in two conditions (visualized Figure 5—figure supplement 2). We identified ligand–receptor interactions that were high in the normal/degenerated group and grouped the differential interactions into ECM, TNF, TGFβ (Transforming growth factor-β), chemokine, and cytokine antigen-presentation categories (Figure 5—figure supplement 2A–L).

Online cell browser construction

We used the python package cellxgene to provide online cell browser service. The matrices along with other metadata were saved in scanpy h5ad format for cellxgene.

Results

The cellular landscape of the inner and outer parts of human menisci

To understand the composition and molecular profiles of the cells from different regions and degeneration states, we sampled four healthy menisci and four degenerated menisci from a patient cohort and performed single-cell RNA sequencing. Figure 1A illustrates the morphological patterns of the samples: Meniscal specimens in the normal group are smooth and complete with a clear structure. Specimens from the degenerated group are swollen, irregular in shape, damaged in structure, and vascularized. Each specimen’s inner part (white–white zone) and the outer part (red–red zone) were strictly separated and collected. These samples are dissociated and prepared for single-cell RNA sequencing (Methods), as shown in the overall workflow (Figure 1B).

Figure 1 with 2 supplements see all

Download asset Open asset

Single-cell RNA-seq reveals major cell classes in human menisci.

(A) Photographs of typical normal and degenerative meniscus specimens. (B) The overall workflow of the single-cell sequencing. Inner and outer meniscal parts are collected separately from patients with normal/degenerated menisci. (**C, D**) Distributions of cells in the inner meniscus samples. Cells of the degeneration batches and the normal batches are integrated together. (C) Uniform manifold approximation and projection (UMAP) visualization of donor origins (inner parts). (D) UMAP visualization of the cell identity compartments (inner parts). (E) The percentages of the identified cell classes in six samples of the inner meniscus. (**F, G**) Distributions of cells in the outer meniscus samples. Cells of the degeneration batches and the normal batches were integrated together. UMAP visualization of the donor origins and cell identity compartments (outer parts). (H) The percentages of the classes in six samples of the outer meniscus. (**I, J**) The expression levels of the major class markers of the inner cells (upper) and the outer cells (lower). Darker colors indicate higher expression levels.

After comprehensive quality control steps removing low-quality cells and inferred doublets (Methods), we recovered 45,744 cells from 12 samples of 8 patients (Supplementary file 1). The number of counts observed per cell and the number of genes observed per cell are visualized in Figure 1—figure supplement 1A, B. We can infer from the quality control steps that the general sequencing quality is good except for the E2 outer and H12 inner samples.

We conducted a survey on the inner and outer menisci’s major cell types. An unsupervised clustering algorithm partitioned the cells into separated clusters (Methods), which were visualized using uniform manifold approximation and projection (UMAP). UMAP and stacked bar plots in Figure 1C–E described the distributions of cells in the inner menisci, and Figure 1F–H described them in the outer menisci. Figure 1C, F visualized the origin of the cells, indicating that the batch effects of the samples were well eliminated, and the biological variations were preserved. The cell type correlation plots in Figure 1—figure supplement 1C also confirmed that cells of the same type, not the same batch, had higher affinities. We observed five major cell types in the inner menisci: chondrocytes, lymphocytes, myeloid cells, endothelial cells, and a group of ACTA2⁺ cells. We also observed the five types in the outer menisci, plus a group of Schwann cells, indicating there are nerve tissues in the outer menisci. The percentages of cells in Figure 1E, H showed that, in general, the outer menisci had higher cell type diversities. In the more diverse outer menisci, the degenerated group (A1, B2, and D2) had higher immune cell ratios and lower ACTA2⁺ cell ratios.

The cell types were determined by combinations of markers in Figure 1I, J and Figure 1—figure supplement 2A, B. We identified cell clusters as chondrocytes if they expressed a high level of cartilage-related genes such as DCN, PRG4, COL1A2, COL3A1, ACAN, COMP, FN1, CHI3L1, and CHI3L2. We identified cell clusters with high expression levels of PLVAP, PECAM1 as endothelial cells. Immune cells were identified using high PTPRC, myeloid markers CD68, C1QA, HLA-DQA1, and T cell markers CD3D, CD3E. In addition, there are a group of cells separated from other cells in both inner and outer areas of the meniscus, which has high expression levels of ACTA2. Besides high expressions of ACTA2, this group of cells also expressed some chondrocyte-associated gene markers, such as COL1A2, COL3A1, DCN, and FN1. We tentatively named this ACTA2⁺ group ‘pericyte-like cells’ or PCL (Figure 1I, J and Figure 1—figure supplement 2A, B) and reanalyzed them together with the chondrocytes in Figure 2 and associated texts.

Figure 2 with 4 supplements see all

Download asset Open asset

Identification of chondrocyte and PCL subclusters in human meniscus.

(A) Uniform manifold approximation and projection (UMAP) visualization of the chondrocyte and PCL class cells. (**B, C**) UMAP visualization of the distribution of chondrocytes at different anatomical sites (up) and different sample statuses (down). (D) A schematic diagram of the hierarchical classification of chondrocyte subgroups. The classification criteria are along the tree path, and each group’s highly expressed marker genes are given below each cluster. (E) A cluster-level heatmap shows extracellular matrix (ECM)-related gene expressions in each chondrocyte/PCL subpopulation. (F) A cell-level heatmap reveals the normalized expression of differentially expressed genes for each cluster defined in (D). (G) Violin plots score nine molecular themes of each cell type (leukocyte migration, suppression of angiogenesis, angiogenesis, response to BMP, collagen fibril, response to IL1, response to TGFβ, response to TNF, and myofibril). We calculated the scores using gene set variation analysis (GSVA) with gene sets picked from MSigDB.

Functional coordination of the heterogeneous chondrocytes and PCL cells

We built a comprehensive cell identity framework by charting the meniscal chondrocytes and PCL cells into a hierarchical coordination system. Unlike the other single-cell cartilage studies that assigned cell types with several marker genes, we defined cell types in a top-down way to simultaneously reveal subtle granular-view variations and exhibit general similarities.

To develop such a cell identity framework, we took multiple factors into account: clustering labels, clustering granularities, marker gene expression specificity, generalization abilities across samples and conditions, supporting evidences in other studies, matrix compositions, and biological functions (Methods). We finally identified five chondral clusters named chondrocyte-1 to chondrocyte-5 (Ch.1–Ch.5 in Figure 2A) and two pericyte-like cell clusters (PCL.1 and PCL.2 in Figure 2A) named ACTA2⁺ cells in Figure 1D, G. The distribution shown in Figure 2B, C shows cells from different batches were mixed well after data integration. The correlations illustrated in Figure 2—figure supplement 1A and the marker gene heatmap in Figure 2—figure supplement 1B indicated that inner and outer samples had consistent expression profiles. From these analyses, we can infer that batch effects are resolved well, and biological variations were retained.

We visualized the cell type identification framework in a tree in Figure 2D, which assigned labels to cells with hierarchical marker combinations. Common fibroblast markers like DCN and LUM are highly expressed in the chondral class. Smooth muscle cell and pericyte markers like ACTA2, TAGLN, and MYL9 are highly expressed in the PCL class. More detailed marker expressions in chondrocyte and PCL subpopulations were visualized in a heatmap in Figure 2F and Figure 2—figure supplement 1C. The marker genes differentially expressed in each subpopulation shown in the heatmap are given in Supplementary file 2.

For more granular cell subtypes, we analyzed their biological functions using GSVA (Hänzelmann et al., 2013) in Figure 3C–F, Figure 2—figure supplement 2H and gene set overrepresentation enrichment analysis (GSEA) (Wu et al., 2021) in Figure 2—figure supplement 2A–G. The gene expressions and function predictions suggested subtypes dominating homeostasis states, subtypes associated with the pro-/anti-angiogenesis process, as well as subtypes associated with the construction/disassembly of the ECM.

Figure 3 with 7 supplements see all

Download asset Open asset

Degeneration molecular patterns in chondrocytes and PCLs.

(A) A volcano plot shows differentially expressed genes with high fold-change values. The comparison was made between (1) degenerated chondrocytes/PCL cells and (2) normal chondrocytes/PCL cells. (B) The proportion of each subcluster of chondrocytes in the degenerated and normal meniscus. (**C–F**) Gene set variation analysis (GSVA) of each cluster under degenerated/normal conditions. Gene sets supported functions are evaluated using GSVA (Methods), and scaled values are visualized in the heatmaps. We observed Ch.2 and Ch.3 had high scores for matrix disassembly activities, angiogenesis activities, chemotaxis activities, and chondrocyte hypertrophy regardless of the normal/degeneration status. The cycling chondrocytes Ch.5 have high scores in the cell-cycle-related terms.

We constructed GRN using pySCENIC (Aibar et al., 2017; Van de Sande et al., 2020) to decipher the common and unique programs behind the chondrocyte and the PCL populations. In the inferred GRN, links originate from the upstream TFs to the downstream targets. Each TF forms a module of nodes that contains its targets. We calculated the module activation scores and showed the binarized states of the gene regulatory programs in Figure 2—figure supplement 3A. Figure 2—figure supplement 3B visualizes a core part of TFs and targets in the GRN. Several polygons were added to annotate subcluster-specific GRN modules. We found that PCLs shared common gene regulatory modules represented by STAT4, NR2F2, and chondrocyte subcluster Ch 2–5 shared TRPS1, FOXC1, HOXD10, and each subcluster possessed some specific modules. For example, Ch.2 had KLF6, FOXN3, AKR1A1, Ch.3 had SIX3, Ch.4 had CEBPA, and Ch.5 had MYBL2. This GRN analysis also explained intra-PCL heterogeneities and indicated PCL cells were a distinct population with divergent gene regulatory programs, although they have some collagen-producing transcriptomic features like chondrocytes.

The characteristics and function predictions of PCL cells in the meniscal microenvironment

The PCL cells present strong muscle contractile gene expression signatures (Figure 2F), including ACTA2, MYL9, TAGLN, etc. ACTA2, which encodes a smooth muscle actin involved in vascular contractility, as well as TAGLN (transgelin) and MYH/MYL (myosin heavy/light chains) genes. These genes have been used for identifying pericytes or smooth muscle cells across multiple human organs, including intervertebral disc (Gan et al., 2021), intestinal tract (Elmentaite et al., 2021), heart (Litvinukova et al., 2020), brain (Smyth et al., 2018), etc. The gene set enrichment analyses of these signature genes suggested the PCL population has a discrepant molecular identity from the chondrocyte lineage because GO terms such as ‘muscle contraction’, ‘muscle organ development’, and ‘smooth muscle cell proliferation’ were enriched (Figure 2—figure supplement 2F).

In the work of Sun et al., 2020, they named a similar ACTA2⁺ MYLK⁺ MCAM⁺ MYL9⁺cluster as the FCP and claimed FCP could differentiate into degenerated meniscus progenitors (DegP) in degenerated menisci, or into adipogenic/osteogenic lineages when cultured ex vivo. We did not name the population we found as FCP for three reasons: Firstly, enrichment analyses showed its similarity to pericytes/smooth muscle cells. Secondly, our PCL population showed weak differentiation potential into DegP because this group of cells in our data did not show concentrated expressions of DegP markers GREM1, CDCP1, and DNER. Thirdly, the only expressions of GREM1, CDCP1, and DNER were sporadically distributed on the UMAP plot (Figure 2—figure supplement 4B) and topologically separated from our ACTA2⁺ MYLK⁺ MCAM⁺ MYL9⁺PCL cluster (Figure 2A, Figure 2—figure supplement 4A).

We named this population as PCL considering multiple factors. The PCL cells are more abundant in the vascularized outer meniscus region than in the avascular inner part (Figure 1E, H). Our immunofluorescence imaging found that this population of cells lay around the vascular endothelial cells and formed a tube shape outside the blood vessels (Figure 4A, B). Like the meniscus, the intervertebral disc also has cartilage tissue and blood vessel structures. Gan et al., 2021 named a group of ACTA2⁺ MCAM⁺ TAGLN⁺ in the human intervertebral disc as pericytes, which can be a reference for defining meniscal cell clusters. For these reasons, we inferred that the ACTA2⁺ MYL9⁺ TAGLN ⁺ population in our meniscus samples belong to the mural structure around the meniscal vasculature and chose to name this cluster as the ‘pericyte-like cell’ considering its signature gene expressions, functional analyses, and spatial distributions. We admit that it is difficult to fully distinguish pericyte from smooth muscle cells purely based on the transcriptomes because the two cell types share numerous markers, and their markers appear to be dynamic and affected by the tissue types and pathogenic states (van Dijk et al., 2015).

Figure 4

Download asset Open asset

Immunofluorescent staining of human menisci demonstrating resident cell subtypes.

Representative immunofluorescence staining images. Color settings: MMP14 for Ch.2 cells (orange), SERPINA1 for Ch.1 cells (magenta), ACTA2 for PCL cells (white), CDON for Ch.3 cells (cyan), CYP1B1 for Ch.4 cells (yellow), CD31 for endothelial cells (green), and CD45 for immune cells (green). Nuclei are stained blue (DAPI). Scale bar 20 μm. (A) Multiplexed IF staining of the normal meniscus. (B) Multiplex IF staining of the degenerated meniscus. (C) IF staining of CDON in a normal sample. (D) IF staining of CDON in a degenerated sample. (E) IF staining of CD31 in a normal sample. (F) IF staining of CD31 in a degenerated sample. (G) IF staining of CD45 in a normal sample. (H) IF staining of CD45 in a degenerated sample (IF: immunofluorescence). The inner and outer zones of menisci were annotated with arrows.

It is known that pericytes attach to the surface of vascular endothelial cells and maintain vasculature stability via the endothelium–pericyte crosstalk (Eilken et al., 2017; Sugihara et al., 2020). We therefore inferred that these cells in the meniscal microenvironments could maintain vascular stability like pericytes according to the GSVA results and the decreased percentages of PCLs in the degeneration group. Previous studies have shown that pericytes have pro-angiogenesis capabilities (Eilken et al., 2017; Armulik et al., 2011). In our datasets, we observed expressions of multiple potential angiogenesis inhibition genes. For example, in Figure 3—figure supplement 1, we showed that both PCL.1 and PCL.2 expressed ECM components Collagen IVα1 (COL4A1), Collagen IVα2 (COL4A2), and Collagen XVIII (COL18A1), the non-collagenous fragments of which (i.e., arresten, canstatin, and endostatin, respectively) were known to have anti-angiogenesis functions (de Castro Brás and Frangogiannis, 2020; Marneros and Olsen, 2001). We also plotted multiple other pro-/anti-angiogenesis factors’ expressions in Figure 3—figure supplement 1B to have a full view of the pro-/anti-effects and found the pro-angiogenesis genes were seldom expressed in PCL either.

We identified two subtle subclusters within the PCL class. One subcluster PCL.1 expresses a higher level of FABP4, RGS5, and the other PCL.2 subcluster has relatively higher expressions of MYH11, MYH9, CNN1, and TAGLN. Kumar et al., 2017 analyzed these genes’ functions in the specification and diversification of pericyte/smooth muscle cells from mesenchymoangioblasts. The GSVA suggested that PCL.1 had more active responses to the pro-inflammatory cytokines IL-1 and TNF-α (Figure 2G) and expressed higher levels of chemokines genes CCL2, CCL19, CCL21, CXCL12 (Figure 3—figure supplement 1D), and COL4A1/2 (Figure 3—figure supplement 1E). The GSVA scores suggested that PCL.2 had higher myofibril gene set activation scores, indicating this population could have higher muscle contractility and may regulate the blood flow in the microvasculature (van Dijk et al., 2015).

Chondrocyte subpopulations and their potential roles in ECM remodeling

The ECMs in meniscal tissues undergo continuous remodeling processes to keep homeostasis of the microenvironment. We aimed to reveal meniscal chondrocyte heterogeneities and to find their associations with the ECM homeostasis and remodeling in normal and degeneration states. With the hierarchical cell type definition framework, we identified two major branches and a small cluster undergoing cell-cycle phase transitions in the chondrocyte class. These branches of cells were predicted to have varied roles in the construction and maintenance of the ECM, angiogenesis, and leukocyte chemotaxis.

The largest branch is the chondrocyte-1 (Ch.1) population that highly expresses CHAD, CHI3L1, CHI3L2, ECRG (C2orf40), and APOD. Some of these genes (e.g., CHI3L1/2) have been described as the markers of RegC (Ji et al., 2019). This cluster’s percentage decreased in the degeneration group. The gene set activation analysis indicated that this population of cells had anti-angiogenesis effects (Figure 2G). In Figure 2E, we visualized the gene expressions that play essential roles in the ECM according to the Matrisome project (Naba et al., 2012). ECM-related gene expressions indicated Ch.1 highly expressed ECM regulators TIMP1 and TIMP4, which were inhibitors of ECM destruction proteins. Hence, we inferred that the CHAD⁺ Ch.1 population is a cluster of chondrocytes that maintain homeostasis in the meniscal tissue. We named another part of this branch that expresses high levels of PRG4, HTRA1, FN1, CD55, TIMP3, CDON, AMTN as Ch.3 (PRG4). We inferred that Ch.3 cells had lubrication function because the protein proteoglycan 4 (PRG4) has a lubrification function on the articular cartilage surface, and these cells were located near the surface of the meniscus. Our immunofluorescence staining further confirmed that its marker CDON expresses along the surface region of the meniscus (Figure 4). This population highly expresses TIMP1/2/3 (Figure 3—figure supplement 1A) and may inhibit the matrix decomposition by opposing the MMP/ADAM/ADAMTS family enzymes.

One branch that highly expresses FNDC1 and THY1 has two subclusters. One TGFBI⁺ subcluster in this branch has particularly high expression levels of FNDC1, so we named it chondrocyte2-FNDC1 (Ch.2). The Ch.2 population highly expresses OA marker genes (Fisch et al., 2018) and chondrocyte hypertrophy genes (Ji et al., 2019) such as type I, III, and VI collagen genes (COL1A1/2, COL3A1, COL6A1/2/3), POSTN, MMP2, SPARC, and MXRA5 (Figure 3—figure supplement 2A, B). Among these upregulated genes, FNDC1 and MXRA5 have been recently reported as signatures of chondrocyte-mediated valve calcification in aortic stenosis (Bouchareb et al., 2021). It is also known that MXRA5 (matrix remodeling-associated protein 5) is a TGF-β1 regulated gene (Poveda et al., 2017), which is consistent with the high GSVA scores of ‘Response to TGFβ’ in Ch.2 shown in Figure 2G. Besides MMP2, this population also expresses other MMP, ADAM, and ADAMTS family genes, such as MMP11/13/14, ADAM9/12, and ADAMTS2, which can target collagen, aggrecans or other ECM components (Figure 3—figure supplement 1A). The increased expression of these genes indicated that Ch.2 (FNDC1) is a chondrocyte population associated with aberrant ECM degradation and remodeling, and this population may contribute to the progress of human meniscal degeneration. The other CFD⁺ subcluster in this branch expresses high levels of CFD, CFH, C1R, C1S, CXCL12, and VCAN. We name it Ch.4 (CFD). This subcluster expresses genes associated with alternative complement activation pathways (C1R, C1S, C2, C3, C6, and C7) and may promote inflammation by recruiting macrophages and neutrophils (Figure 3—figure supplement 1C).

Ch.5 represents chondrocytes in proliferation states. They highly expressed the S- and G2M-phase cell-cycle genes such as STMN1, CDK1, CENPF, HMGB2, CDKN3, TOP2A and were mostly found in the degenerated samples. It is generally thought that cells during mitosis condense their chromatins, downregulated genes that maintain cell identities, and reestablish the gene expression patterns when the mitosis exits (Palozola et al., 2017). Hence, it is hard to judge whether Ch.5 is a high-stemness chondrocyte progenitor cluster or a dividing differentiated cell cluster based on these highly expressed markers. We may assume that the cells, though in a proliferation state, retained some transcripts before mitosis, and we may infer their origins or fates via the reference mapping analysis. We re-assigned Ch.5 cell identities using cells in Ch.1–Ch.4 clusters as a reference with SingleR (Aran et al., 2019) and found 98.63% Ch.5 cells were mapped to Ch.2 and the rest were mapped to Ch.3 (Figure 3—figure supplement 7F). This identity mapping suggested the proliferative Ch.5 could be the origin of the Ch.2 cells increased in the degenerated group.

We did not directly adopt the chondrocyte terminology used in other studies such as Ji et al., 2019 and Sun et al., 2020, but instead used Ch.1–Ch.5 to represent the subclusters because our data spans chondrocytes of higher quantity and greater heterogeneity. We have examined the literature markers on our data in Figure 2—figure supplement 4A, B. Those markers did not show highly concentrated expressional patterns. We reanalyzed the data in Sun et al., 2020 and validated our markers (Figure 3—figure supplement 2C–E). We found Ch.1 markers CHI3L1/2 expressed on clusters 1, 2, 7, Ch.2 markers FNDC1 expressed on clusters 0 and 3, and cycling cell markers STMN1 expressed on clusters 5 and 6. The expressions of FNDC1, TGFBI, PRG4, CFD, and CFH were highly overlapped with each other, which is slightly different with our data (Figure 3—figure supplement 2E). We performed a comprehensive integration of our data and the public data in Sun et al., 2020 using Seurat (Stuart et al., 2019). Figure 3—figure supplement 3A, B shows the transferred UMAP coordinates after integration. From the comprehensive integration results, we inferred that our data spanned broader transcriptomic states and the public data mainly corresponded to a subset near Ch.2, Ch.3, and Ch.5 locations. We calculated the percentages of transferred labels in the public clusters (Figure 3—figure supplement 3C–E) and revealed their mutual correspondences. We found the ACTA2⁺ cluster 3 in public data partially mapped to the PCL cells in our data. We used a hierarchically clustered joint heatmap of our data and public data to visualize the cluster correspondences (Figure 3—figure supplement 3F). Considering these integration results, we believe our cluster labels and markers span a larger state space and could better represent cellular heterogeneities.

Degeneration and zonation-associated cell subpopulation changes

We compared the overall gene expression differences between the normal and degenerations conditions for chondral and PCL populations and derived a series of conditional DEGs. We identified genes upregulated in the degeneration group, such as COL1A1, COL1A2, COL3A1, COL6A1, POSTN, SPARC, CILP, TGFBI, DPT, and ATMN, and genes upregulated in the normal group, such as APOD, CHI3L1, CHI3L2, CEBPD, C11orf96, C2orf40 (ECRG4), ACTA2, TAGLN, BCAM, and MCAM (shown in Figure 3A and Figure 3—figure supplement 4, an extended list of genes in Figure 3—figure supplement 5A). The enrichment analysis interpreting the potential functions of these DEGs is available in Figure 3—figure supplement 5B, C. We believe the overall expression differences are mainly contributed by the cell type compositional changes because these genes are identified as markers of the chondrocyte subpopulations. For example, the expressions of Ch.2 marker genes COL1A1/2, COL3A1, COL6A1, and SPARC are higher in the degenerated group; and Ch.3 markers AMTN and CILP also have higher expressions in the degenerated group; whereas PCL cell markers ACTA2, TAGLN, BCAM, and MCAM are highly expressed in the normal group. All these expression differences are consistent with our compositional observations that Ch.2 and Ch.3 increased in degeneration and Ch.1 and PCL decreased in degeneration (Figure 3B and Figure 3—figure supplement 6A). Besides overall differences, we also calculated the DEGs in each cluster between the degenerated and normal conditions (Figure 3—figure supplement 6B) and found more subtle expressional shifts in the degeneration states. Considering all factors above, it appears that the overall expressional changes after degeneration could be explained by the increase of Ch.2 and Ch.3 as well as the decrease of PCL cells.

To explore each subcluster’s biofunction in degeneration, we performed GSVA analyses of each population in two conditions. We first analyzed the GSVA scores of ECM-related gene ontology terms in Figure 3C and found Ch.2 could have high ECM disassembly scores while Ch.3 could negatively regulate the disassembly regardless of the degeneration states. They could play essential roles in the anabolic/catabolic balance. Ch.2 and Ch.3 also got high scores in leukocyte chemotaxis (Figure 3D), angiogenesis (Figure 3E), and chondrocyte hypertrophy-associated terms (Figure 3F); while in contrast, Ch.1 had relatively lower scores on these terms. This is consistent with the fact that Ch.1 cells dominate the normal meniscus (Figure 3B), and leukocyte infiltration and angiogenesis mainly happen in the degenerated meniscus. Ch.5 got high scores on the cell-cycle-related terms (Figure 3F). PCL cells got higher pro-angiogenesis scores in degenerated samples than in normal samples. The GSVA analyses indicated the increased or decreased subclusters in degenerated samples had pro- or anti-degeneration functions and further supported that the compositional changes of subclusters primarily contributed to the degeneration process.

It is acknowledged that the inner and outer zones have varied recovery abilities, ECM compositions, cellular compositions, and mechanical characteristics. But it is still unknown which zone is the primary contributor to joint inflammation in OA. We analyzed the zonation-dependent cellular state shift and computationally evaluated inflammation response scores in two regions. From the cell density plots in Figure 3—figure supplement 7A, B, we found that the inner zone chondrocytes shift from the Ch.1 positions to the Ch.3 positions, while the outer zone chondrocytes shift to Ch.2 and Ch.4 positions. Both zones have cell state shifts after degeneration. We compared the cell densities instead of directly comparing the subpopulation percentages. Since the molecular profiles of subclusters vary on a continuum and the clustering boundaries could be vague, directly comparing the percentages with paired/unpaired statistical tests could miss the state shifts (Figure 3—figure supplement 7C, D). Compared to Ch.3, the Ch.2 cells have higher scores on the responses to inflammation cytokines TNF and IL1 (Figure 2—figure supplement 2H). We also evaluated the GSVA scores of three inflammation gene sets on the chondrocytes and found the outer-degenerated groups consistently got significantly higher inflammation scores than the outer-normal group. In contrast, the inner zone cells did not have this trend (Figure 3—figure supplement 7E). From the computational analyses, we inferred that outer menisci could contribute more to joint inflammation in OA.

Spatially resolve cell subtype distributions with in situ imaging

We used immunofluorescence staining to validate the subpopulation markers and study the spatial distribution of the subclusters in the human meniscus. Among the top 20 marker genes (according to log2 fold-change values) for each chondrocyte subpopulation, we chose SERPINA1, MMP14, CDON, CYP1B1, ACTA2, CD31, and CD45 for staining. Normal and degenerated meniscus sample slides were simultaneously stained with the following antibodies: anti-SERPINA1 for Ch1(CHAD), anti-MMP14 for Ch2(FNDC1), anti-CYPIB1 for Ch4(CFD), and anti-ACTA2 for PCL cells. We used anti-CD31, anti-CDON, and anti-CD45 to label vascular endothelial cells, Ch3(PRG4), and immune cells for immunofluorescence staining. From the staining results, we can clearly find that CDON + chondrocytes were mainly distributed on the surface of the meniscus tissue, and ACTA2⁺ pericytes were presented surrounding blood vessels. SERPINA1⁺, CYP1B1⁺, and MMP14⁺ chondrocytes were mixed and scattered inside the meniscus. Immune cells were mainly distributed around the blood vessels in the meniscus. Compared with the normal samples, the number of ACTA2⁺ chondrocytes was decreased, the number of MMP14⁺ chondrocytes was markedly increased, and the number of vascular endothelial cells was significantly increased in degenerative specimens (Figure 4A–H).

Crosstalk changes of chondrocytes, leukocytes, and endothelial cells in meniscus degeneration

Besides chondrocytes, various types of leukocytes, and endothelial cells also have indispensable roles in the occurrence and development of cartilage degeneration and OA (Hsueh et al., 2021; Alahdal et al., 2021).

For the leukocyte part, we identified macrophages/monocytes, neutrophils, mast cells, dendritic cells, T cells, and cycling immune cells in normal and degenerated meniscal samples (Figure 5A, B, Figure 5—figure supplement 1A). Macrophage/monocyte makes up the largest portion of the total immune population. The percentages of macrophage/monocyte, T cell, dendritic cell, and total immune cell increased in the degeneration group, which was consistent with the vascularization and the damage-induced inflammation (Figure 5B, C). We noticed some pro-inflammatory cytokine/chemokine genes upregulated in degeneration samples, for example, interleukin genes IL1A, IL1B, IL6, IL15, and IL18 upregulated in the cycling subpopulation (Figure 5A); chemokine genes CXCL1, CCL3, and CXCL3 upregulated in macrophages/monocytes (Figure 5—figure supplement 1D). We also noticed increased expressions of metallothionein genes in the degenerated group, for example, MT1X and MT1G in macrophages/monocytes (Figure 5—figure supplement 1B), MT1X, MT2A, and MT1E in T cells (Figure 5—figure supplement 1C).

Figure 5 with 2 supplements see all

Download asset Open asset

Immune, endothelial cells, and their crosstalk with chondrocytes.

(A) Uniform manifold approximation and projection (UMAP) visualization of the immune cell types in the meniscus. (B) UMAP visualization of the statuses of samples (degenerated/normal). (C) Per-sample bar plots visualize the immune cell percentage changes between the normal and degenerated group. (D) Changes in immune cell percent between normal and degenerated groups. Error bars show the standard deviations of the data. Wilcoxon p values: Mac.Mo 0.004, Cycling 0.10, DC 0.14, Neu 0.19, T 0.22, Mast 1.0. (E) UMAP visualization of the endothelial class cells’ subtypes. (F) UMAP visualization of sample status. (G) UMAP visualization of the anatomical regions. (H) Top differentially expressed genes between different health states (degenerated vs. normal). The heatmap shows z-score-scaled gene expression values. (I) General cell–cell crosstalk between large populations. The ligand–receptor pair crosstalk was evaluated at the large population level (Methods). (J) Representative crosstalk was significantly enhanced in the degeneration group.

We identified an arterial–capillary–venous gradient in the meniscal endothelial cell population (Figure 5E–G) with other organs’ endothelial markers (Schupp et al., 2021; Kalucka et al., 2020). The capillary–arterial endothelial cells highly express GJA5; the capillary endothelial cells highly expressed PLVAP; and the capillary–venous endothelial cells highly expressed VWF. These endothelial cells comprise the microvascular and capillary vessels. A tiny cluster of lymphatic endothelial cells was also observed, which may come from the lymphatic vessels in the meniscus. We compared the total endothelial cells between the normal/degenerated groups and reported a list of DEGs (Figure 5H). In the normal group, CST3, ACKR1, RBP7, and IGFBP4 were upregulated, while in the degenerated group, S100A4, DDIT4, CRIP1, and MGP were upregulated. CST3 is ranked first in the normal group DEGs, which was reported to have inhibition effects on the proliferation, migration, tube formation, and permeability of endothelial cells (Li et al., 2018). ACKR1 is also highly expressed in the normal group, which encodes a chemokine decoy receptor inhibiting the effectiveness of other chemokines (Rot, 2005). These molecules contributed to the blood vessel’s stability in the normal group.

With all major classes of meniscal cells being identified, we inferred the cell–cell interactions (CCIs) based on the ligand and receptor expressions and compared the interaction intensities between the normal/degenerated samples. The results suggested that CCI intensities generally increased in the degeneration group, especially for the endothelial cells and immune cells (Figure 5I). These enhanced CCIs are also consistent with the angiogenesis and leukocyte recruitment phenomena in the degeneration group (Figure 6). Figure 5J shows some instances of the significant CCIs originating from the sender cell type to the receiver cell type. For example, in the degenerated meniscus, we observed that a chemokine gene CXCL8 was upregulated by macrophages/monocytes, neutrophils, DCs; and its receptor gene SDC2 was expressed in chondrocytes and the pericyte-like cells. We also observed significant CALM1/2-INSR CCIs occurred in the degeneration in which angiogenesis happened – many cell types upregulated CALM1/2, and endothelial cells upregulated the insulin receptor gene INSR. This observation is consistent with a previous study showing insulin receptors’ pro-angiogenesis functions. (Walker et al., 2021). We also identified multiple other cell–cell crosstalks enhanced in the degeneration group that are involved in immunoregulations. They are mediated by TNF superfamily interactions (TNF, TNFSF10 → VSIR, TNFRSF1A/B), and TGF superfamily interactions (INHBA → ENG), and Class I HLA → APLP2 interactions. In Figure 5—figure supplement 2A–L, we provided a detailed resource of the predicted CCIs. We illustrated sender–receiver CCI pairs high in the normal or degeneration groups. Note that all the CCIs mentioned above were derived using gene expressions and curated lists of ligands and receptors. We did not examine these computational inferences with experiments.

Figure 6

Download asset Open asset

A schematic diagram of the microenvironment changes between the normal and the degenerated meniscus.

The left side visualizes the homeostasis meniscus, where Ch.1 is the dominant chondrocyte population. In this situation, the extracellular matrix (ECM) decomposition and synthesis reach a dynamic equilibrium, and the aberrant proliferation of blood vessels is inhibited. The right side visualizes the degenerated meniscus, where the orchestrated microenvironment balance is broken. The pericyte-like cells detach from the blood vessels and lower the stability of the blood vessels. The endothelial cells grew and formed new blood vessels in the degenerated area. The vascular permeability increased and allowed more immune cells to infiltrate the degenerated tissue. Ch.2 and Ch.3 increased in quantity and produced more matrix disassembly enzymes, damaging the ECM, releasing angiogenesis factors, and recruiting more immune cells. The immune cells, such as macrophage/monocytes, DCs, and T cells, also have chemotaxis functions. There are also cycling chondrocytes and nerve tissue cells in the meniscus.

Based on the above observations, we inferred that the immune cells and endothelial cells could be important regulators for chondrocytes and ECMs (Figure 6). It has been known that the degenerated ECM released damage-associated molecular patterns (DAMPs) that initiated inflammatory responses (Martel-Pelletier et al., 2016; Chen et al., 2017a). According to our inference, immune cells could sense the DAMP signals, migrate to the inflammatory meniscus and produce cytokines and chemokines. The inflammation could in turn enhance ECM catabolic enzymes like MMP/ADAM/ADAMTS, stimulate blood vessel proliferation, and further promote immune infiltration in the meniscal tissue. These factors work in a manner like a ‘positive feedback loop’. In common conditions, the DAMP-induced inflammation will reverse when the DAMP molecules are cleared (Zindel and Kubes, 2020). However, if the above factors work in the positive feedback loop in some perturbed conditions, they can mutually promote each other and drive the microenvironment away from homeostasis.

Discussion

Recent single-cell studies have shown the importance of meniscal progenitors in tissue engineering and given a rough picture of the meniscal cell types. However, the existing data were still insufficient to reveal degeneration participants in the microenvironment. In this work, we investigated meniscal chondrocytes, developed a hierarchical classification system, and classified the cells into five major populations (Ch.1–Ch.5). We identified endothelial cells enriched in both outer meniscal parts and inner parts, which are typically conceived as avascular ‘white–white zone’. We inferred that these unanticipated cells were the degeneration product because newly generated vascular endothelial cells usually occur in cartilage as OA progresses (Goldring and Goldring, 2016). Surrounding the microvessels are the pericyte-like cells possessing muscle contraction features. We also observed Schwann cells, which typically wrap around the neurons, indicating there are nerve tissues in the human meniscus. For the immune cells, we observed monocytes/macrophages, DCs, T cells, neutrophils, mast cells, and some cycling immune cells.

The meniscal microenvironment changes significantly from the homeostasis when the degeneration progresses (Figure 6). In normal meniscal tissues, the synthesis and disassembly of the matrix molecules were balanced in the orchestrated ECM metabolism. We inferred that the dominant Ch.1 population was essential to the balance maintenance because they expressed TIMP family molecules that inhibited the ECM decomposers (MMP/ADAM/ADAMTS). The PCL population that adhered to the surface of endothelial cells maintained vasculature stability. The high-level expressions of atypical chemokine receptors like ACKR1 neutralized the chemokines and attenuated the angiogenesis (Xu et al., 2007). The patrolling immune cells mainly pass by the tissue in the blood vessel without residing or being attracted.

However, the balanced microenvironment was altered significantly in the degenerated meniscal tissues. The expression of ECM decomposers (MMP/ADAM/ADAMTS family molecules) from chondral subpopulations Ch.2 and Ch.3 might surpass the inhibition forces (TIMP) and decompose the ECM by breaking large PG molecules into fragments. We believe the blood vessels have two sides in the homeostasis of the meniscus. On the one hand, the blood vessels can help to restore acute injuries to health; on the other hand, the proliferated blood vessels can build pathways for the immune cells so that they can infiltrate the meniscal matrix and cause inflammation. We observed multiple synergistic factors promoting angiogenesis and increased endothelial percentages in the degenerated samples. While the pericyte stabilizes the blood vessel in health, the detachment and decrease of the pericyte may lead to pathological states (Figueiredo et al., 2020; Armulik et al., 2011; Diéguez-Hurtado et al., 2019). On the vascular wall, the endothelial cells could downregulate ACKR1 expressions, respond to the pro-angiogenesis factors like VEGFA, FGF, IGF1, TGFβ produced by the Ch.2, Ch.3 Ch.4, grew new branches, increased the vascular permeability. Changing the chondrocyte subtype ratios and limiting blood vessel proliferation may delay or reverse the meniscus degeneration.

In addition to these findings from the analyses of the meniscal single-cell sequencing data, we have also made an interactive website to store and let readers to browse the processed data. Readers can directly access the website through the link http://meni.singlecell.info:3000/, and perform online analyses of the results, comparison of differential genes between two special clusters, or data downloading.

In the design of this study, we only partitioned samples into two health states: normal and degenerated. If we distinguish the disease severity of the meniscus samples, we could establish more fine-grained links between the molecular profiles and the disease. Another limitation of this study is that females make up the majority of our donors (Supplementary file 1), since OA is more common in middle-aged and older females. The cellular profiles observed in this study may better represent the female population. We observed a small number of nerve tissue cells, which also hindered detailed analysis of them. Our functional coordination of chondrocyte subtypes was confined to the meniscus environment, even though there were signs that chondrocytes across different cartilage tissues were not distinct and had many common features. We hope a future molecular coordination framework could integrate chondrocytes across multiple cartilage tissues, different health states, anatomical regions, and differentiation stages, guiding the precise taxonomy of chondrocytes when more single-cell cartilage studies are released.

In conclusion, we systematically profiled the cellular diversities in the inner and outer meniscus and reported the microenvironmental alterations in the healthy and degenerated states. Our study is an informative complement to the existing meniscal single-cell sequencing data and provides an important reference for the study of meniscal degeneration. The study suggested that we should pay attention to the blood vessel’s functional duality in acute and chronic inflammation. Preserving the pericyte-like cells wrapping around the vessels might be a potential strategy to alleviate chronic degeneration. Systematically coordinating of angiogenesis, inflammation response, ECM catabolism, and preventing their mutual reinforcement could be a strategy for delaying or reversing meniscus degeneration. The meniscus is an important and representative fibrocartilage in the human body. We speculated that the cellular and molecular mechanisms observed in this study on meniscal degeneration, such as angiogenesis and inflammation, could also apply to intervertebral discs and articular cartilage degeneration. It could shed new light on the diagnosis and treatment of other degeneration-associated musculoskeletal diseases like low back pain and joint pain.

Data availability

Data are available in a public, open access repository. The single-cell RNA-seq data, cluster annotations are available at GSA for human (https://ngdc.cncb.ac.cn/gsa-human/) with the accession number PRJCA008120.

The following data sets were generated

1. Weili Fu1, Sijie Chen, Runze Yang, Chen Li, Haoxiang Gao, Jian Li, Xuegong Zhang
(2022) BioProject
ID PRJCA008120. Single-cell transcriptomic profiling of healthy and degenerated human menisci.

https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA008120

References

1. Aibar S
2. González-Blas CB
3. Moerman T
4. Huynh-Thu VA
5. Imrichova H
6. Hulselmans G
7. Rambow F
8. Marine J-C
9. Geurts P
10. Aerts J
11. van den Oord J
12. Atak ZK
13. Wouters J
14. Aerts S
(2017) Scenic: single-cell regulatory network inference and clustering
Nature Methods 14:1083–1086.

https://doi.org/10.1038/nmeth.4463
- PubMed
- Google Scholar
1. Alahdal M
2. Zhang H
3. Huang R
4. Sun W
5. Deng Z
6. Duan L
7. Ouyang H
8. Wang D
(2021) Potential efficacy of dendritic cell immunomodulation in the treatment of osteoarthritis
Rheumatology 60:507–517.

https://doi.org/10.1093/rheumatology/keaa745
- PubMed
- Google Scholar
1. Aran D
2. Looney AP
3. Liu L
4. Wu E
5. Fong V
6. Hsu A
7. Chak S
8. Naikawadi RP
9. Wolters PJ
10. Abate AR
11. Butte AJ
12. Bhattacharya M
(2019) Reference-Based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage
Nature Immunology 20:163–172.

https://doi.org/10.1038/s41590-018-0276-y
- PubMed
- Google Scholar
(2011) Pericytes: developmental, physiological, and pathological perspectives, problems, and promises
Developmental Cell 21:193–215.

https://doi.org/10.1016/j.devcel.2011.07.001
- PubMed
- Google Scholar
1. Baran Y
2. Bercovich A
3. Sebe-Pedros A
4. Lubling Y
5. Giladi A
6. Chomsky E
7. Meir Z
8. Hoichman M
9. Lifshitz A
10. Tanay A
(2019) MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions
Genome Biology 20:206.

https://doi.org/10.1186/s13059-019-1812-2
- PubMed
- Google Scholar
1. Beaufils P
2. Becker R
3. Kopf S
4. Englund M
5. Verdonk R
6. Ollivier M
7. Seil R
(2017) Surgical management of degenerative meniscus lesions: the 2016 ESSKA meniscus consensus
Joints 5:59–69.

https://doi.org/10.1055/s-0037-1603813
- PubMed
- Google Scholar
(2021) Proteomic architecture of valvular extracellular matrix: FNDC1 and MXRA5 are new biomarkers of aortic stenosis
JACC. Basic to Translational Science 6:25–39.

https://doi.org/10.1016/j.jacbts.2020.11.008
- PubMed
- Google Scholar
1. Butler A
2. Hoffman P
3. Smibert P
4. Papalexi E
5. Satija R
(2018) Integrating single-cell transcriptomic data across different conditions, technologies, and species
Nature Biotechnology 36:411–420.

https://doi.org/10.1038/nbt.4096
- PubMed
- Google Scholar
1. Chen S
2. Fu P
3. Wu H
4. Pei M
(2017a) Meniscus, articular cartilage and nucleus pulposus: a comparative review of cartilage-like tissues in anatomy, development and function
Cell and Tissue Research 370:53–70.

https://doi.org/10.1007/s00441-017-2613-0
- PubMed
- Google Scholar
1. Chen D
2. Shen J
3. Zhao W
4. Wang T
5. Han L
6. Hamilton JL
7. Im HJ
(2017b) Osteoarthritis: toward a comprehensive understanding of pathological mechanism
Bone Research 5:16044.

https://doi.org/10.1038/boneres.2016.44
- PubMed
- Google Scholar
(1983)
Blood supply to the normal and abnormal menisci of the human knee

Clinical Orthopaedics and Related Research 1:271–276.
- PubMed
- Google Scholar
1. de Castro Brás LE
2. Frangogiannis NG
(2020) Extracellular matrix-derived peptides in tissue remodeling and fibrosis
Matrix Biology 91–92:176–187.

https://doi.org/10.1016/j.matbio.2020.04.006
- PubMed
- Google Scholar
1. Denisenko E
2. Guo BB
3. Jones M
4. Hou R
5. de Kock L
6. Lassmann T
7. Poppe D
8. Clément O
9. Simmons RK
10. Lister R
11. Forrest ARR
(2020) Systematic assessment of tissue dissociation and storage biases in single-cell and single-nucleus RNA-seq workflows
Genome Biology 21:130.

https://doi.org/10.1186/s13059-020-02048-6
- PubMed
- Google Scholar
1. Diéguez-Hurtado R
2. Kato K
3. Giaimo BD
4. Nieminen-Kelhä M
5. Arf H
6. Ferrante F
7. Bartkuhn M
8. Zimmermann T
9. Bixel MG
10. Eilken HM
11. Adams S
12. Borggrefe T
13. Vajkoczy P
14. Adams RH
(2019) Loss of the transcription factor RBPJ induces disease-promoting properties in brain pericytes
Nature Communications 10:2817.

https://doi.org/10.1038/s41467-019-10643-w
- PubMed
- Google Scholar
1. Eilken HM
2. Diéguez-Hurtado R
3. Schmidt I
4. Nakayama M
5. Jeong HW
6. Arf H
7. Adams S
8. Ferrara N
9. Adams RH
(2017) Pericytes regulate VEGF-induced endothelial sprouting through VEGFR1
Nature Communications 8:1574.

https://doi.org/10.1038/s41467-017-01738-3
- PubMed
- Google Scholar
1. Elmentaite R
2. Kumasaka N
3. Roberts K
4. Fleming A
5. Dann E
6. King HW
7. Kleshchevnikov V
8. Dabrowska M
9. Pritchard S
10. Bolt L
11. Vieira SF
12. Mamanova L
13. Huang N
14. Perrone F
15. Goh Kai’En I
16. Lisgo SN
17. Katan M
18. Leonard S
19. Oliver TRW
20. Hook CE
21. Nayak K
22. Campos LS
23. Domínguez Conde C
24. Stephenson E
25. Engelbert J
26. Botting RA
27. Polanski K
28. van Dongen S
29. Patel M
30. Morgan MD
31. Marioni JC
32. Bayraktar OA
33. Meyer KB
34. He X
35. Barker RA
36. Uhlig HH
37. Mahbubani KT
38. Saeb-Parsy K
39. Zilbauer M
40. Clatworthy MR
41. Haniffa M
42. James KR
43. Teichmann SA
(2021) Cells of the human intestinal tract mapped across space and time
Nature 597:250–255.

https://doi.org/10.1038/s41586-021-03852-1
- PubMed
- Google Scholar
1. Englund M
2. Roemer FW
3. Hayashi D
4. Crema MD
5. Guermazi A
(2012) Meniscus pathology, osteoarthritis and the treatment controversy
Nature Reviews. Rheumatology 8:412–419.

https://doi.org/10.1038/nrrheum.2012.69
- PubMed
- Google Scholar
(2020) Phosphoinositide 3-kinase-regulated pericyte maturation governs vascular remodeling
Circulation 142:688–704.

https://doi.org/10.1161/CIRCULATIONAHA.119.042354
- PubMed
- Google Scholar
1. Fisch KM
2. Gamini R
3. Alvarez-Garcia O
4. Akagi R
5. Saito M
6. Muramatsu Y
7. Sasho T
8. Koziol JA
9. Su AI
10. Lotz MK
(2018) Identification of transcription factors responsible for dysregulated networks in human osteoarthritis cartilage by global gene expression analysis
Osteoarthritis and Cartilage 26:1531–1538.

https://doi.org/10.1016/j.joca.2018.07.012
- PubMed
- Google Scholar
1. Gan Y
2. He J
3. Zhu J
4. Xu Z
5. Wang Z
6. Yan J
7. Hu O
8. Bai Z
9. Chen L
10. Xie Y
11. Jin M
12. Huang S
13. Liu B
14. Liu P
(2021) Spatially defined single-cell transcriptional profiling characterizes diverse chondrocyte subtypes and nucleus pulposus progenitors in human intervertebral discs
Bone Research 9:37.

https://doi.org/10.1038/s41413-021-00163-z
- PubMed
- Google Scholar
(1983)
Ultrastructure of normal and torn menisci of the human knee joint

Journal of Anatomy 136:773–791.
- PubMed
- Google Scholar
1. Goldring SR
2. Goldring MB
(2016) Changes in the osteochondral unit during osteoarthritis: structure, function and cartilage-bone crosstalk
Nature Reviews. Rheumatology 12:632–644.

https://doi.org/10.1038/nrrheum.2016.148
- PubMed
- Google Scholar
1. Gu Z
2. Eils R
3. Schlesner M
(2016) Complex heatmaps reveal patterns and correlations in multidimensional genomic data
Bioinformatics 32:2847–2849.

https://doi.org/10.1093/bioinformatics/btw313
- PubMed
- Google Scholar
(2013) GSVA: gene set variation analysis for microarray and RNA-Seq data
BMC Bioinformatics 14:7.

https://doi.org/10.1186/1471-2105-14-7
- PubMed
- Google Scholar
1. Hsueh MF
2. Zhang X
3. Wellman SS
4. Bolognesi MP
5. Kraus VB
(2021) Synergistic roles of macrophages and neutrophils in osteoarthritis progression
Arthritis & Rheumatology 73:89–99.

https://doi.org/10.1002/art.41486
- PubMed
- Google Scholar
1. Hua K
2. Zhang X
(2019) A case study on the detailed reproducibility of a human cell atlas project
Quantitative Biology 7:162–169.

https://doi.org/10.1007/s40484-018-0164-3
- Google Scholar
1. Ji Q
2. Zheng Y
3. Zhang G
4. Hu Y
5. Fan X
6. Hou Y
7. Wen L
8. Li L
9. Xu Y
10. Wang Y
11. Tang F
(2019) Single-Cell RNA-seq analysis reveals the progression of human osteoarthritis
Annals of the Rheumatic Diseases 78:100–110.

https://doi.org/10.1136/annrheumdis-2017-212863
- PubMed
- Google Scholar
1. Kalucka J
2. de Rooij LPMH
3. Goveia J
4. Rohlenova K
5. Dumas SJ
6. Meta E
7. Conchinha NV
8. Taverna F
9. Teuwen L-A
10. Veys K
11. García-Caballero M
12. Khan S
13. Geldhof V
14. Sokol L
15. Chen R
16. Treps L
17. Borri M
18. de Zeeuw P
19. Dubois C
20. Karakach TK
21. Falkenberg KD
22. Parys M
23. Yin X
24. Vinckier S
25. Du Y
26. Fenton RA
27. Schoonjans L
28. Dewerchin M
29. Eelen G
30. Thienpont B
31. Lin L
32. Bolund L
33. Li X
34. Luo Y
35. Carmeliet P
(2020) Single-cell transcriptome atlas of murine endothelial cells
Cell 180:764–779.

https://doi.org/10.1016/j.cell.2020.01.015
- PubMed
- Google Scholar
1. Korsunsky I
2. Millard N
3. Fan J
4. Slowikowski K
5. Zhang F
6. Wei K
7. Baglaenko Y
8. Brenner M
9. Loh PR
10. Raychaudhuri S
(2019) Fast, sensitive and accurate integration of single-cell data with harmony
Nature Methods 16:1289–1296.

https://doi.org/10.1038/s41592-019-0619-0
- PubMed
- Google Scholar
1. Kumar A
2. D’Souza SS
3. Moskvin OV
4. Toh H
5. Wang B
6. Zhang J
7. Swanson S
8. Guo L-W
9. Thomson JA
10. Slukvin II
(2017) Specification and diversification of pericytes and smooth muscle cells from mesenchymoangioblasts
Cell Reports 19:1902–1916.

https://doi.org/10.1016/j.celrep.2017.05.019
- PubMed
- Google Scholar
1. Li Z
2. Wang S
3. Huo X
4. Yu H
5. Lu J
6. Zhang S
7. Li X
8. Cao Q
9. Li C
10. Guo M
11. Lv J
12. Du X
13. Chen Z
(2018) Cystatin C expression is promoted by VEGFA blocking, with inhibitory effects on endothelial cell angiogenic functions including proliferation, migration, and chorioallantoic membrane angiogenesis
Journal of the American Heart Association 7:e009167.

https://doi.org/10.1161/JAHA.118.009167
- PubMed
- Google Scholar
(2015) The molecular signatures database (msigdb) hallmark gene set collection
Cell Systems 1:417–425.

https://doi.org/10.1016/j.cels.2015.12.004
- PubMed
- Google Scholar
1. Litvinukova M
2. Talavera-lopez C
3. Maatz H
4. Reichart D
5. Worth CL
6. Lindberg EL
7. Kanda M
8. Polanski K
9. Heinig M
10. Lee M
11. Nadelmann ER
12. Roberts K
13. Tuck L
14. Fasouli ES
15. Delaughter DM
16. Mcdonough B
17. Wakimoto H
18. Gorham JM
19. Samari S
20. Mahbubani KT
21. Saeb-parsy K
22. Patone G
23. Boyle JJ
24. Zhang H
25. Zhang H
(2020) Cells of the adult human heart
Nature 588:466–472.

https://doi.org/10.1038/s41586-020-2797-4
- Google Scholar
(2007) The long-term consequence of anterior cruciate ligament and meniscus injuries: osteoarthritis
The American Journal of Sports Medicine 35:1756–1769.

https://doi.org/10.1177/0363546507307396
- PubMed
- Google Scholar
(2011) The knee meniscus: structure-function, pathophysiology, current repair techniques, and prospects for regeneration
Biomaterials 32:7411–7431.

https://doi.org/10.1016/j.biomaterials.2011.06.037
- PubMed
- Google Scholar
1. Marneros AG
2. Olsen BR
(2001) The role of collagen-derived proteolytic fragments in angiogenesis
Matrix Biology 20:337–345.

https://doi.org/10.1016/s0945-053x(01)00151-2
- PubMed
- Google Scholar
(2016) Osteoarthritis
Nature Reviews. Disease Primers 2:16072.

https://doi.org/10.1038/nrdp.2016.72
- PubMed
- Google Scholar
(2019) DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors
Cell Systems 8:329–337.

https://doi.org/10.1016/j.cels.2019.03.003
- PubMed
- Google Scholar
(2019) The meniscus in normal and osteoarthritic tissues: facing the structure property challenges and current treatment trends
Annual Review of Biomedical Engineering 21:495–521.

https://doi.org/10.1146/annurev-bioeng-060418-052547
- PubMed
- Google Scholar
1. Naba A
2. Clauser KR
3. Hoersch S
4. Liu H
5. Carr SA
6. Hynes RO
(2012) The matrisome: in silico definition and in vivo characterization by proteomics of normal and tumor extracellular matrices
Molecular & Cellular Proteomics 11:M111.

https://doi.org/10.1074/mcp.M111.014647
- PubMed
- Google Scholar
(1989) Mechanics of the healed meniscus in a canine model
The American Journal of Sports Medicine 17:164–175.

https://doi.org/10.1177/036354658901700205
- PubMed
- Google Scholar
1. Palozola KC
2. Donahue G
3. Liu H
4. Grant GR
5. Becker JS
6. Cote A
7. Yu H
8. Raj A
9. Zaret KS
(2017) Mitotic transcription and waves of gene reactivation during mitotic exit
Science 358:119–122.

https://doi.org/10.1126/science.aal4671
- PubMed
- Google Scholar
(2017) MXRA5 is a TGF-β1-regulated human protein with anti-inflammatory and anti-fibrotic properties
Journal of Cellular and Molecular Medicine 21:154–164.

https://doi.org/10.1111/jcmm.12953
- PubMed
- Google Scholar
1. Rai MF
2. McNulty AL
(2017) Meniscus beyond mechanics: using biology to advance our understanding of meniscus injury and treatment
Connective Tissue Research 58:221–224.

https://doi.org/10.1080/03008207.2017.1312921
- PubMed
- Google Scholar
1. Rot A
(2005) Contribution of Duffy antigen to chemokine function
Cytokine & Growth Factor Reviews 16:687–694.

https://doi.org/10.1016/j.cytogfr.2005.05.011
- PubMed
- Google Scholar
1. Schupp JC
2. Adams TS
3. Cosme C
4. Raredon MSB
5. Yuan Y
6. Omote N
7. Poli S
8. Chioccioli M
9. Rose K-A
10. Manning EP
11. Sauler M
12. DeIuliis G
13. Ahangari F
14. Neumark N
15. Habermann AC
16. Gutierrez AJ
17. Bui LT
18. Lafyatis R
19. Pierce RW
20. Meyer KB
21. Nawijn MC
22. Teichmann SA
23. Banovich NE
24. Kropski JA
25. Niklason LE
26. Pe’er D
27. Yan X
28. Homer RJ
29. Rosas IO
30. Kaminski N
(2021) Integrated single-cell atlas of endothelial cells of the human lung
Circulation 144:286–302.

https://doi.org/10.1161/CIRCULATIONAHA.120.052318
- PubMed
- Google Scholar
(2013) Meniscus repair and regeneration: review on current methods and research potential
European Cells & Materials 26:150–170.

https://doi.org/10.22203/ecm.v026a11
- PubMed
- Google Scholar
1. Shannon P
2. Markiel A
3. Ozier O
4. Baliga NS
5. Wang JT
6. Ramage D
7. Amin N
8. Schwikowski B
9. Ideker T
(2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks
Genome Research 13:2498–2504.

https://doi.org/10.1101/gr.1239303
- PubMed
- Google Scholar
1. Shao X
2. Liao J
3. Li C
4. Lu X
5. Cheng J
6. Fan X
(2021) CellTalkDB: a manually curated database of ligand-receptor interactions in humans and mice
Briefings in Bioinformatics 22:bbaa269.

https://doi.org/10.1093/bib/bbaa269
- PubMed
- Google Scholar
(2018) Markers for human brain pericytes and smooth muscle cells
Journal of Chemical Neuroanatomy 92:48–60.

https://doi.org/10.1016/j.jchemneu.2018.06.001
- PubMed
- Google Scholar
1. Stuart T
2. Butler A
3. Hoffman P
4. Hafemeister C
5. Papalexi E
6. Mauck WM
7. Hao Y
8. Stoeckius M
9. Smibert P
10. Satija R
(2019) Comprehensive integration of single-cell data
Cell 177:1888–1902.

https://doi.org/10.1016/j.cell.2019.05.031
- PubMed
- Google Scholar
1. Subramanian A
2. Tamayo P
3. Mootha VK
4. Mukherjee S
5. Ebert BL
6. Gillette MA
7. Paulovich A
8. Pomeroy SL
9. Golub TR
10. Lander ES
11. Mesirov JP
(2005) Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles
PNAS 102:15545–15550.

https://doi.org/10.1073/pnas.0506580102
- PubMed
- Google Scholar
1. Sugihara K
2. Sasaki S
3. Uemura A
4. Kidoaki S
5. Miura T
(2020) Mechanisms of endothelial cell coverage by pericytes: computational modelling of cell wrapping and in vitro experiments
Journal of the Royal Society, Interface 17:20190739.

https://doi.org/10.1098/rsif.2019.0739
- PubMed
- Google Scholar
1. Sun Y
2. Mauerhan DR
3. Honeycutt PR
4. Kneisl JS
5. Norton JH
6. Hanley EN Jr
7. Gruber HE
(2010) Analysis of meniscal degeneration and meniscal gene expression
BMC Musculoskeletal Disorders 11:19.

https://doi.org/10.1186/1471-2474-11-19
- PubMed
- Google Scholar
1. Sun H
2. Wen X
3. Li H
4. Wu P
5. Gu M
6. Zhao X
7. Zhang Z
8. Hu S
9. Mao G
10. Ma R
11. Liao W
12. Zhang Z
(2020) Single-Cell RNA-seq analysis identifies meniscus progenitors and reveals the progression of meniscus degeneration
Annals of the Rheumatic Diseases 79:408–417.

https://doi.org/10.1136/annrheumdis-2019-215926
- PubMed
- Google Scholar
1. Van de Sande B
2. Flerin C
3. Davie K
4. De Waegeneer M
5. Hulselmans G
6. Aibar S
7. Seurinck R
8. Saelens W
9. Cannoodt R
10. Rouchon Q
11. Verbeiren T
12. De Maeyer D
13. Reumers J
14. Saeys Y
15. Aerts S
(2020) A scalable scenic workflow for single-cell gene regulatory network analysis
Nature Protocols 15:2247–2276.

https://doi.org/10.1038/s41596-020-0336-2
- PubMed
- Google Scholar
(2017) Single-Cell sequencing reveals dissociation-induced gene expression in tissue subpopulations
Nature Methods 14:935–936.

https://doi.org/10.1038/nmeth.4437
- PubMed
- Google Scholar
1. van Dijk CGM
2. Nieuweboer FE
3. Pei JY
4. Xu YJ
5. Burgisser P
6. van Mulligen E
7. el Azzouzi H
8. Duncker DJ
9. Verhaar MC
10. Cheng C
(2015) The complex mural cell: pericyte function in health and disease
International Journal of Cardiology 190:75–89.

https://doi.org/10.1016/j.ijcard.2015.03.258
- PubMed
- Google Scholar
1. Walker AMN
2. Warmke N
3. Mercer B
4. Watt NT
5. Mughal R
6. Smith J
7. Galloway S
8. Haywood NJ
9. Soomro T
10. Griffin KJ
11. Wheatcroft SB
12. Yuldasheva NY
13. Beech DJ
14. Carmeliet P
15. Kearney MT
16. Cubbon RM
(2021) Endothelial insulin receptors promote VEGF-A signaling via ERK1/2 and sprouting angiogenesis
Endocrinology 162:bqab104.

https://doi.org/10.1210/endocr/bqab104
- PubMed
- Google Scholar
(2018) SCANPY: large-scale single-cell gene expression data analysis
Genome Biology 19:15.

https://doi.org/10.1186/s13059-017-1382-0
- PubMed
- Google Scholar
1. Wu T
2. Hu E
3. Xu S
4. Chen M
5. Guo P
6. Dai Z
7. Feng T
8. Zhou L
9. Tang W
10. Zhan L
11. Fu X
12. Liu S
13. Bo X
14. Yu G
(2021) ClusterProfiler 4.0: A universal enrichment tool for interpreting omics data
Innovation 2:100141.

https://doi.org/10.1016/j.xinn.2021.100141
- PubMed
- Google Scholar
(2007) Duffy antigen/receptor for chemokines (DARC) attenuates angiogenesis by causing senescence in endothelial cells
Angiogenesis 10:307–318.

https://doi.org/10.1007/s10456-007-9084-y
- PubMed
- Google Scholar
1. Young MD
2. Behjati S
(2020) SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data
GigaScience 9:giaa151.

https://doi.org/10.1093/gigascience/giaa151
- PubMed
- Google Scholar
1. Zindel J
2. Kubes P
(2020) Damps, PAMPS, and lamps in immunity and sterile inflammation
Annual Review of Pathology 15:493–518.

https://doi.org/10.1146/annurev-pathmechdis-012419-032847
- PubMed
- Google Scholar
1. Zou Z
2. Hua K
3. Zhang X
(2021) HGC: fast hierarchical clustering for large-scale single-cell data
Bioinformatics 1:btab420.

https://doi.org/10.1093/bioinformatics/btab420
- PubMed
- Google Scholar

Article and author information

Author details

Weili Fu

Department of Orthopedics, Orthopedic Research Institute, West China Hospital, Sichuan University, Chengdu, China

Contribution
Conceptualization, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing

Contributed equally with
Sijie Chen

For correspondence
foxwin2008@163.com

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4438-2760
Sijie Chen

MOE Key Lab of Bioinformatics, Bioinformatics Division, BNRIST and Department of Automation, Tsinghua University, Beijing, China

Contribution
Conceptualization, Resources, Data curation, Software, Funding acquisition, Visualization, Methodology, Writing – original draft

Contributed equally with
Weili Fu

Competing interests
No competing interests declared
Runze Yang

Department of Orthopedics, Orthopedic Research Institute, West China Hospital, Sichuan University, Chengdu, China

Contribution
Formal analysis, Validation, Writing – original draft

Competing interests
No competing interests declared
Chen Li

MOE Key Lab of Bioinformatics, Bioinformatics Division, BNRIST and Department of Automation, Tsinghua University, Beijing, China

Contribution
Software, Visualization, Methodology

Competing interests
No competing interests declared
Haoxiang Gao

MOE Key Lab of Bioinformatics, Bioinformatics Division, BNRIST and Department of Automation, Tsinghua University, Beijing, China

Contribution
Software

Competing interests
No competing interests declared
Jian Li

Department of Orthopedics, Orthopedic Research Institute, West China Hospital, Sichuan University, Chengdu, China

Contribution
Conceptualization, Writing – review and editing

Competing interests
No competing interests declared
Xuegong Zhang
1. MOE Key Lab of Bioinformatics, Bioinformatics Division, BNRIST and Department of Automation, Tsinghua University, Beijing, China
2. School of Life Sciences and School of Medicine, Center for Synthetic and Systems Biology, Tsinghua University, Beijing, China
Contribution
Conceptualization, Data curation, Supervision, Funding acquisition, Project administration, Writing – review and editing

For correspondence
zhangxg@tsinghua.edu.cn

Competing interests
No competing interests declared

Funding

National Natural Science Foundation of China (81972123)

Weili Fu

National Natural Science Foundation of China (82172508)

Weili Fu

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors are grateful to Dr. Tao Wu from the Institute of Immunology, Tsinghua University, for his helpful discussions. The authors are grateful to Kaibo Zhang, Zhong Zhang, Beini Mao, Sike Lai, and Menglin Yao for the help of cell isolations and sample preparations. The authors are grateful to Yi Zhang, Yue Li, and Wanli Zhang for the help of multiplex immunofluorescence staining.

The experiments were approved by Ethics Committee on Biomedical Research, West China Hospital of Sichuan University with an identifier of No. 2020-921. All individuals have been informed of the experiments and analyses, and they have given consent to use their tissue samples for this research and publication.

Ethics

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.