Single-Cell RNA Sequencing Reveals CXCLs Enriched Fibroblasts Within Odontogenic Keratocysts

Introduction

Odontogenic keratocyst (OKC) is a locally invasive jaw cystic lesion characterized by a relatively high growth potential and a tendency towards recurrence.¹ The histology of OKC usually shows a thin parakeratinized epithelium, with columnar cells in the basal layer possessing focal reverse polarization.¹ Genetic and molecular researches had revealed the presence of PTCH1 mutations and the activation of the Hedgehog signaling pathway related to OKC pathogenesis.^1,2 The unusual behavior of OKC, grow along the long axis of the jaw and relatively high recurrence, suggests that aberrant molecular events behind its local invasiveness.³ Angiogenesis and cellular invasion are complex multistage processes that begin with the degradation of the extracellular matrix caused by proteolytic enzymes, and, especially, tumor-associated matrix metalloproteinases.⁴

Chemokines are important transmitter for leukocyte migration to inflammatory sites which play significant roles in inflammation.⁵ Inflammation had been frequently detected in OKCs.^6–8 Previous study had illustrated that a focal increase of proliferating epithelial cells were detected adjacent to moderate to severe inflammation in OKC, but without a significant effect on the overall proliferative activity of the cysts.⁶ Compared to classical OKC and dentigerous cyst (DC), inflamed OKC and inflamed DC showed significantly higher proliferating cell nuclear antigen expression.⁷ The thicker lining epithelia was also found in OKC with focal inflammation sites than those without inflammation.⁸ Therefore, inflammation might be closely related with the focal proliferation of epithelia. Recent years, chemokines are thought to be involved in endothelial cells functions especially angiogenesis.⁹ Chemokine receptors are now proved to be expressed on endothelial cells and chemokines could bind to the receptors and lead to the activation of these receptors which finally lead to various cellular responses including endothelial cells migration, proliferation, and tube formation.⁹

Mechanical stress had been vital in musculoskeletal health and disease. In addition to the genetic factors, intracystic fluid pressure had been long present in the odontogenic cystic lesions and affected the development of odontogenic cystic lesions.^10,11 Marsupialization, which reduced intracystic fluid pressure, had shown obvious clinical benefits in OKC.¹² New bone formations within OKC could be detected after marsupialization.¹² Intracystic fluid pressure had been reported to induced IL-1α expression and relating to the development of OKC.¹¹ Mechanical stress is now also emerging as a possible trigger for some forms of chronic inflammatory lesions, such as periodontitis and rheumatoid arthritis.^13,14 However, the potential relation between mechanical stress and chemokine release within OKC remains unclear.

Body tissue comprises a heterogeneous population of various cells, secreted factors and extracellular matrix proteins, together forming the tissue microenvironment.¹⁵ Tissue was influenced and regulated by various extracellular signals including cytokines, chemokines which are conserved between mice and humans at transcriptional levels.¹⁵ Heterogeneity is a fundamental and intrinsic characteristic underlining the presence of various cell clusters, but the discrepancy is covered when unlabelled cells are used for omic analysis.¹⁶ Single-cell RNA sequencing (scRNA-seq) represents an innovative method to discover the transcriptional heterogeneity in health and diseases, revealing new insights into tissue composition, and disease progression,¹⁷ and actually, the discovery of the cell populations that contribute to OKC pathogenesis is in urgent need.

In this study, single-cell RNA sequencing was performed using fresh OKC tissue samples. The result revealed that CXCLs were highly enriched in the subset of fibroblasts and showed close interactions with endothelial cells and myeloid cells which were partially regulated by mechanical stress. We investigated the pathological mechanism of OKC and providing a new insight into the pathogenesis of OKC.

Materials and Methods

Tissue Dissociation and Preparation

Three OKC tissue samples were collected from the Department of Oral and Maxillofacial Head Neck Surgery, School & Hospital of Stomatology, Wuhan University. The patients agreed and signed the informed consent form, and ethical approval was acquired from the Ethics Committee of the School of Stomatology, Wuhan University (Wuhan, PR China). The processing of clinical samples was conducted in accordance with the Declaration of Helsinki. Specimens of OKC were carefully washed with Hanks solution for 5 times and cutted into 2–3 mm pieces. Then, 2 mL GEXSCOPE^® Tissue Dissociation Solution was used for the tissue pieces digestion at 37°C for 15 mins in 15 mL tubes. The solutions were filtered using 40-micron sterile strainers and then centrifuged at 200 × g for 3 mins. After centrifugation, the supernatant was carefully removed, and the deposit was gently suspended in 2 mL phosphate buffered saline (PBS) (Servicebio, Wuhan, China). Red blood cells were removed using 3 mL red blood cell lysis buffer (Servicebio, Wuhan, China) which was incubated with the collected cells at 37°C for 5 mins. Then, the buffer was removed by centrifugation with speed of 200 × g for 3 mins. The pellet was suspended in 2 mL PBS and stained by trypan blue (Servicebio, Wuhan, China). The stained suspension was observed and evaluated in microscope.

Single-Cell RNA Sequencing

Single-cell suspension was diluted to the concentration of 1 × 10⁵ cells/mL in PBS (Servicebio, Wuhan, China). Microfluidic chips were used and the single-cell suspension was added to the microfluidic chips. The scRNA-seq libraries were established using GEXSCOPE^® Single-Cell RNA Library Kit (Singleron Biotechnologies, Nanjing, China). Individual libraries were prepared using 4nM and then pooled for the sequencing. Pools were later sequenced on Illumina HiSeq X using 150 bp paired end reads.

Analysis of Raw Read Data

FastQC was used for the process of raw reads and fastp was used to remove low quality reads.¹⁸ Cutadapt was employed to remove Poly-A tails and adaptor sequences. After quality control, STAR (version 2.5.2b) was used to map reads to the reference genome GRCh38 (ensembl version 92 annotation). Feature Counts software was used for the acquirement of Gene counts and Unique Molecular Identifier (UMI) counts. Based on the gene counts and UMI counts, files of expression matrix were obtained for subsequent analyses. The sequencing data of the present study was acquired and deposited in the Gene Expression Omnibus (GEO) database (accession number GSE176351).

Quality Control, Dimension-Reduction and Clustering

Cells were screened by gene counts between 300 and 4000. In addition, cells were also filtered by UMI counts below 30,000. Cells which had over 30% mitochondrial content were dropped out. After filtering, OKC cells were reserved for the subsequent analyses. Seurat R package (v3.0.2) was used for dimension-reduction and clustering.¹⁹ NormalizeData() and ScaleData() was used for the normalization of gene expressions. FindVariableFeautres() was used for top 2000 variable genes selection for principal component analysis. As results, basal cells were characterized into five subsets by FindClusters, with the top 20 principle components and resolution parameter at 1.0. For subclustering of basal cells and fibroblasts, the resolution was set to 1.2.

Differentially Expressed Genes (DEGs) Analysis and Cell Type Annotation

Seurat FindMarkers() was employed for DEG analysis. Specifically, genes which were expressed in greater than 10% of the total cells within a cluster and with average log (Fold Change) of larger than 0.25 were regarded as DEGs. The cell type identity of each cluster was determined with the expression of canonical markers found in the DEGs. Heatmaps, dot plots and violin plots showing the expression of markers used to identify each cell type were generated by Seurat DotPlot.

Pathway Enrichment and Gene Set Enrichment Analysis

To investigate the potential functions of subtypes of basal cells, the Gene Ontology (GO) analysis was used with the “clusterProfiler” R package version.²⁰ The P_adj values smaller than 0.05 were regarded as significantly enriched pathways. The activity of selected genes in cells were evaluated by AUCell version package (version 1.2.4). To determine cells with the given active gene set, the assignment threshold was set to 0.2 for each gene set in each cell.

Trajectory Analysis

To map differentiation/conversion of cell subtypes of fibroblasts during shift from fibroblasts #1 to fibroblasts #3, pseudo time trajectory analysis was performed with Monocle R package (version 2.10.1).²¹ For constructing the trajectory, top highly variable genes were selected from different clusters by Seurat. The trajectory was visualized by plot cell trajectory.

Cell–Cell Interaction Analysis

Cell–cell interactions between fibroblasts #3 and other cell types were predicted based on known ligand–receptor pairs by Cellphone DB version 2.0.²² Permutation number for calculating the null distribution of average ligand-receptor pair expression in randomized cell identities was set to 1000. Individual ligand or receptor expression was thresholded by a cutoff based on the average log gene expression distribution for all genes across each cell type. Predicted interaction pairs with p value < 0.05 and of average log expression > 0.1 were considered as significant and visualized by circlize R package (v 0.4.10).

Primary Culture of OKC Fibroblast Culture and Hydrostatic Pressure Treatment in Culture

The primary fibroblast culture was performed as described previously.²³ The tissue of OKC was immediately rinsed by PBS (Servicebio, Wuhan, China) and cutted into 5 × 5 mm cubes which were placed on the cell culture dishes. DMEM (Servicebio, Wuhan, China) medium supplemented with 10% FBS and 1% penicillin and streptomycin (Servicebio, Wuhan, China). The medium was changed every 3 days. After 7 days, spindle shape cells were present in the periphery of the tissue of OKC. The cells were cultured and seeded at 1 × 10⁵ cells/well in six-well culture plates. Press 200 System (Shanghai Naturethink Life & Scientific CO., Ltd) was used for the hydrostatic pressure treatment. This computer-regulated bioreactor applies hydrostatic pressure to the cultured cells (80 mmHg, 6 hours). Through setting atmospheric pressure, the compressive forces were applied to the cultured cells within the incubator.

Cytokine Antibody Array

The samples for cytokines and chemokines detection were stored as previously described.²³ Supernatant of cultured fibroblasts (5 × 10⁵/well in six-well culture plates) was collected and stored in −80°C. Following the manufacturer’s protocols, detection of cytokines and chemokines was performed by Bio-Plex Pro Human Cytokine Screening 40-plex panel using Bio-PlexMAGPIX System (Bio-Rad) (Wayen Biotechnologies (Shanghai), Inc).

Immunohistochemical Staining of CXCL12 and CD31 in OKC

Immunohistochemical staining was performed as described previously.²³ Tissue microarray (TMA) of 30 cases of primary OKC was collected (Servicebio, Wuhan, China). Briefly, tissue slides were dewaxed, hydrated, washed and endogenous peroxidase neutralized. Then, tissue sections with 4 μm thick were incubated with goat serum at 37°C for 20 mins. Anti-CXCL12 (1:100; Proteintech, Wuhan, China) and CD31 (1:200; Proteintech, Wuhan, China) were incubated with slides overnight with at 4°C. The second antibodies were incubated and followed by DAB staining (Zhongshan golden bridge, Beijing, China). The nuclei were counterstained with hematoxylin (Servicebio, Wuhan, China).

Scoring System and Analysis

Aperio Image Scope CS2 scanner (CA, USA) was used for TMA slides scan. Aperio Quantification software (Version 9.1) was applied for the quantification. The stromal area was manually selected for quantification. Then, the formula (1 × the percentage of weakly positive staining) + (2 × the percentage of moderately positive staining) + (3 × the percentage of strongly positive staining) was applied to count H-scores of stromal area. The scaled values (H-scores) of CXCL12 and CD31 expressions were subsequently converted to the dataset (Microsoft Excel). Then, the stromal expressions (H-scores) of CXCL12 and CD31 were assessed by the Spearman correlation.

Statistical Analysis

The data were represented as mean ± SD. Student’s t-test was applied to analyze the data between control group and experimental group. Correlations between two variables were assessed by the Spearman correlation. P < 0.05 was considered as statistically significant differences. GraphPad Prism 7 were used for all statistical analyses.

Results

Single-Cell Sequencing and Cell Type Identification

Fresh OKC tissue samples were divided into two parts: one was used for hematoxylin eosin staining, the other one for scRNA-seq. Clinical characteristics of OKC are listed in Table 1. As regard the scRNA-seq, fresh tissue was digested until obtaining a single-cell suspension containing 1.5 × 10⁵, 1.0 × 10⁵ and 1.3 × 10⁵ cells and the viability of cells was 81.60%, 83.50% and 84.3%, respectively. After quality control, 14,072 single cells were clustered into eight major clusters. The cell types were annotated by cluster-specific genes combined with classic markers: basal cells (KRT5+; 8019 cells); endothelial cells (CDH5+; 1450 cells); fibroblasts (LUM+; 2242 cells); pericytes (RGS5+; 258 cells); myeloid cells (LYZ+; 1737 cells); T cells (CD3D+; 186 cells); melanocytes (MLANA+; 37 cells) and plasma cells (MZB1+; 143 cells) (Figure 1A and B). Then, the gene-expression heatmap was shown to reveal the top differentially expressed genes in each group (Figure 1C). Cell components of OKC showed high heterogeneity among the donors (Figure 1D). The KRT5+ basal cells were clustered again into 5 groups (Figure 1E). Fibroblasts (LUM+) were clustered into three different types (Figure 1F).

Table 1 Clinical Features of Three OKC Patients

Figure 1 Identification of cell types in odontogenic keratocyst (OKC) tissue samples. (A) UMAP plot of single cells profiled in the presenting work colored by major cell types (basal cells, fibroblasts, myeloid cells, endothelial cells, T cells, plasma cells, pericytes and melanocytes); (B) Cell number of major cell types in OKC; (C) Heatmap of the top expressed genes in major cell types. (D) Proportion of major cell types in 3 patients with OKC (cases 1, 2, and 3). (E) Re-clustering of basal cells (KRT5+) into 5 subtypes represented by UMAP plot. (F) Re-clustering of fibroblasts (LUM+) into 3 subtypes represented by UMAP plot.

Fibroblasts Heterogeneity in the Odontogenic Keratocyst Tissue Samples

Fibroblasts (LUM+) were clustered into three clusters. SFRP4 and SPP1 were highly enriched in fibroblast #1 (Figure 2A and Table 2). Fibroblast #3 exhibited a strong enrichment in the genes of various chemokines, including CXCL1, CXCL13, CXCL12, CXCL6 and CXCL8 (Figure 2A, Table 3). Fibroblast #2 exhibited a strong enrichment in the genes of cell proliferation, including H2AFZ, TOP2A, and MKI67 (Table 4). The differentiation trajectory revealed that the fibroblast #3 was distributed at the beginning (Figure 2B). As regard to distribution of the clusters among patients, we found that fibroblast subsets varied in the cases, but fibroblast #3 account for more than 50% and was the most common type (Figure 2C). GO analysis was performed on the DEGs of fibroblast #3 and others to investigate the function of each subgroup. Fibroblast #1 was related to extracellular matrix constituent and glycosaminoglycan binding (Figure 2D). Fibroblast #2 was related to tubulin binding and microtubule binding (Figure 2E). Fibroblast #3 showed enriched genes in cytokine activity and chemokine activity (Figure 2F). These results indicated that fibroblast #3 was the most common cell types for OKC fibroblasts and were enriched in chemokine and cytokine genes.

Table 2 Differentially Expressed Genes (Top20) in Fibroblast 1#

Table 3 Differentially Expressed Genes (Top20) in Fibroblast 3#

Table 4 Differentially Expressed Genes (Top20) in Fibroblast 2#

Figure 2 Three different fibroblast subtypes identified in odontogenic keratocyst (OKC). (A) Differential gene expression panel from fibroblasts #1 into fibroblasts #3. (B) Trajectory analysis of the differentiation from fibroblasts #1 into fibroblasts #3 in 3 OKC samples (cases 1, 2, and 3). (C) Proportion of major clusters of fibroblasts in 3 patients with OKC (cases 1, 2, and 3). (D–F) Enriched GO analysis of the up-regulated genes in fibroblasts #1 to fibroblasts #3. Colored circle indicates P values of differential genes among the fibroblast clusters.

Mechanical Stress Regulated Genes Were Enriched in Fibroblast #3

To further explore the potential characteristics of fibroblasts subsets, the cytoskeleton-related genes including TAGLN2 and integrin (ITGA8) were explored and they were highly expressed in the fibroblast #3 indicating the potential mechanical regulation in the fibroblast #3 (Figure 3A). Odontogenic keratocysts are usually semi-fluids contained lesions and the intracystic fluid pressure had been examined and showed about 80 mmHg at the initial stage. To mimic the potential effects of hydrostatic pressure, we employed the compressive loading system and applied 80 mmHg pressure on the fibroblasts for 6 h (Figure 3B). The cytokine antibody array was used for the test of cytokines, especially chemokines, within the supernatant of the cultured fibroblasts. The results showed that CXCL12 were significantly upregulated in the hydrostatic pressure treated fibroblasts compared with control group (P<0.01, Figure 3C). But other chemokines (CCL21, CXCL5, CX3CL1, CXCL1, CXCL8, CCL2, CCL25) showed no significant changes (Figure 3C). The CXCL5 chemokines are the most high in concentrations (over 2000 pg/mL) among the chemokines (Figure 3C).

Figure 3 Mechanical stress regulated CXCL12 secretion in fibroblasts. (A) Dot plot showing the expression of TAGLN2and ITGA8 in 3 fibroblast clusters (fibroblasts #1, #2, and #3). The red dotted rectangle represent average relative expressions >0.5 in 3 fibroblast clusters. (B) Diagram of hydrostatic pressure device used to apply pressure on the cultured fibroblasts of odontogenic keratocyst (OKC). (C) Chemokine profiles (CCL21, CXCL5, CX3CL1, CXCL1, CXCL8, CCL2, CXCL12, and CCL25) from the supernatant of the hydrostatic pressure (80 mmHg, 6h) treated fibroblasts of OKC. Graphs in C represent the Mean ± SD. The green histogram represent values in control cultures. Differences were analyzed using Student’s t-test. **P < 0.01.

Fibroblast #3 is Closely Related to Angiogenesis of OKC

To further evaluate the receptor of CXCLs, we evidenced that atypical Chemokine Receptor 1 (ACKR1) was enriched in endothelial cells (Figure 4A). The cell–cell interaction network among the cell types identified in this study was investigated using CellphoneDB v2.0. Notably, fibroblast #3 showed specific strong interactions with endothelial cells (Figure 4B). To further validate the CXCLs expression levels, we constructed tissue microarray of OKC (n = 30), we evaluated the CXCL12 expressions using immunohistochemistry, and we found CXCL12 were mainly expressed in the fibroblasts and immune cells (Figure 4C). By statistical analysis of H-scores of CXCL12 and CD31, the stromal expressions of CXCL12 which were closely correlated with CD31 expressions (R²= 0.4853, P < 0.001) (Figure 4D). These results indicated that CXCL12 expressions within the stroma might be closely related to the angiogenesis of OKC.

Figure 4 Chemokine enriched fibroblasts were closely related to endothelial cells. (A) Violin plot showing the expression of ACKR1 in major cell types. (B) Heatmap showing the number of potential ligand-receptor pairs among cell groups predicted by CellphoneDB v2.0. (C) CXCL12 and CD31 staining in tissue microarray of odontogenic keratocyst (OKC). (D) Stromal CXCL12 was closely related to CD31 staining in OKC (R2= 0.4853; P < 0.001). Scale bar for (C) = 200 µm.

Discussion

Our scRNA-seq data identified heterogeneous cell populations within OKC. The CXCLs were enriched in the subset of fibroblasts (fibroblast #3, chemokine enriched fibroblasts). We evidenced that CXCLs enriched fibroblasts were closely related with the endothelial cells and immune cells. Furthermore, hydrostatic pressure, which was frequently present in odontogenic lesions could induce the CXCL12 secretion in OKC fibroblasts. A regulatory network was also constructed based on fibroblast #3 in a bioinformatic way. Due to lack of the single-cell analysis of OKC, the results were firstly used to elucidate the comprehensive cellular regulatory network of OKC.

Heterogeneity is a fundamental and intrinsic characteristic underlining the presence of various cell populations.²⁴ However, due to the lack of detection method at single levels, the heterogeneity of OKC had long been unexplored. scRNA-seq enables the exploration of the transcriptome of single cells which provide the opportunities for cell types identification and cellular heterogeneity recognition. Presently, we had identified 8 major cell types including basal cells, fibroblasts, myeloid cells, endothelial cells, pericytes, T cells, plasma cells and melanocytes. The most common cell types were basal cells accounting for about 50% of total cell populations. Fibroblast proportions varied in individuals. But fibroblast #3 (chemokine enriched fibroblasts) was the most common subset of fibroblasts populations indicating the indispensable roles of chemokine enriched fibroblasts in OKC.

OKCs are frequently local inflamed being caused by microbial infection or other unknown factors.^3,4 Even more, mechanical stress, including compressive and tensile force, triggers a local immunological response, inflammation and bone resorption.^13,14 In addition, marsupialization which reduced the hydrostatic pressure within OKC significantly boost the new bone formation and lower the proliferation of basal cells.¹² The chemokine enriched fibroblasts were also highly expressed mechanical stress related genes (TAGLN2 and ITGA8) indicating the potential close relation between mechanical stress and chemokine release. Presently, it was shown that fibroblasts produced chemokine (CXCL12) in response to hydrostatic pressure as mechanical stress. In a study with 14 radicular cysts, the found that the initial intracystic pressure was 47 mmHg and the odontogenic keratocyst was 80 mmHg.^11,25 In another study, the intracystic pressure of odontogenic keratocyst ranged from about 20 mmHg to 120 mmHg in 9 OKC patients.¹⁰ Oka et al found that the pressure of 80 mmHg caused IL-1α mRNA elevations and increased the secretion of MMP-1, MMP-2, MMP-3, and PGE2 in the fibroblasts.¹¹ We employed intracystic pressure (80mmHg) to mimic the physical environment in vitro for OKC study.

Cytokine profiles had been investigated in the cyst fluids and tissue of odontogenic lesions. Kolokythas et al reported that aspirate contained IL-1 and TNF-α were relatively enriched in the OKC group (n = 6) indicating the potential association with epithelial proliferation.²⁶ Dezerega et al reported that monocyte chemotactic protein-3 (MCP-3) was expressed in apical lesions relating to clinical symptoms. MCP-3 might be involved in stimulating mononuclear chemotaxis.²⁷ Cavalla et al reported apical lesions showed elevated CXCL12 expressions and were related to persistent recruitment of immune cells, particularly of mast cells.²⁸ However, studies regarding chemokines in tissue of odontogenic keratocyst were rare.

Although many studies had reported the functions of CXCL12-CXCR4 in various diseases,²⁹ only a few reports had described CXCL12 secretion regulated by mechanical forces.³⁰ CXCL12, also named stromal cell derived factor 1α (SDF-1α), is a vital player in angiogenesis, inflammation, and embryonic development.³¹ Study showed that autocrine CXCL12 could increase migration ability and tubule formation ability of human umbilical vein endothelial cells, without VEGF, through activating the CXCR7 receptor.³² Importantly, several chemokines such as CXCL1, CXCL8 and CXCL12 induce angiogenesis directly through binding to G-protein-coupled chemokine receptors located on endothelial cells, resulting in elevated angiogenesis ability.³³ In murine cortical osteocytes and periosteal cells, mechanical loading could significantly upregulated the SDF-1 expressions.³⁴ Others have revealed that mechanical stretch induces expression of CXCL1 in liver sinusoidal endothelial cells to enroll neutrophils.³⁵ Static stretch significantly increased mRNA expressions of iNOS, IL-6, and MCP-1 in the culture of colonic circular smooth muscle cells.³⁶ Until now, no studies had reported CXCL12 expressions in OKC and CXCL12 secretions under the hydrostatic pressure. Herein, we treated primary OKC fibroblasts using 80 mmHg and found significant elevations of CXCL12 in the supernatant of OKC fibroblasts. By CellPhoneDB analysis, the communication relationship between chemokine enriched fibroblasts, myeloid cells and endothelial cells was confirmed. This refers to the angiogenic roles of chemokine-related fibroblasts which was also evidenced by close relations between tissue expression of CXCL12 and vascular endothelial marker CD31.

However, there are some limitations in this study. Firstly, small cases of OKC are not enough to infer about molecular subtypes of OKC. Especially the PTCH gene mutations should be considered in the further studies. Secondly, the post-translational modifications were not explored in this study. The transcriptional level and the protein level should be studied in broad spectrum at the same time. Thirdly, short-term and long-term hydrostatic pressure effects should be explored in the future to elucidate the mechanism of OKC pathogenesis.

Conclusion

In conclusion, this work is the first study that fully characterized the landscape of OKC cells at a genomic level by scRNA-seq. Our work illustrated the role of one fibroblast subset, named chemokine enriched fibroblasts, and the close relationship with angiogenesis and immune infiltrations which was partially regulated by hydrostatic pressure. Our results revealed new insight on OKC cell subtypes and the relationship between physical microenvironment and angiogenesis of OKC. Thus, these results facilitated the understanding of the heterogeneity within OKC lesions and provide new pathogenesis for OKC.