scRNA-Seq: First Atlas and Cellular Landscape of Lacrimal Sac: Implications in Primary Acquired Nasolacrimal Duct Obstruction Pathogenesis

Abstract

Purpose: The aim of this study was to describe the transcriptional changes of individual cellular components in the lacrimal sac in patients with primary acquired nasolacrimal duct obstruction (PANDO) and attempt to construct the first lacrimal sac cellular atlas to elucidate the potential mechanisms that may drive the disease pathogenesis.

Methods: Lacrimal sac samples were obtained intra-operatively during the endoscopic dacryocystorhinostomy (EnDCR) procedure from five patients. Single-cell RNA sequencing was performed to analyze each individual cell population including epithelial and immune cells during the early inflammatory and late inflammatory phases of the disease.

Results: Eleven cell types were identified among 25,791 cells. T cells and B cells were the cell populations with the greatest variation in cell numbers between the two phases and were involved in immune response and epithelium migration-related pathways. The present study showed that epithelial cells highly expressed the genes of senescence-associated secretory phenotype (SASP) and were involved in influencing the inflammation, neutrophil chemotaxis, and migration during the late inflammatory stage. Enhanced activity of CXCLs-CXCRs between the epithelial cells and neutrophils was noted by the cell-cell communication analysis and is suspected to play a role in inflammation by recruiting more neutrophils.

Conclusions: The study presents a comprehensive single-cell landscape of the lacrimal sac cells in different phases of PANDO. The contribution of T cells, B cells, and epithelial cells to the inflammatory response, and construction of the intercellular signaling networks between the cells within the lacrimal sac has further enhanced the present understanding of the PANDO pathogenesis.

Primary acquired nasolacrimal duct obstruction (PANDO) is a common disease of the lacrimal drainage system with a predilection in perimenopausal female patients.¹ It presents as epiphora with or without discharge, and is usually managed by dacryocystorhinostomy.^2,3 The etiopathogenesis of PANDO is elusive but several factors have been implicated, including anatomic factors, vascular factors, autonomic dysregulation, microbial factors, local hormonal imbalance, changes in the tear composition, lysosomal dysregulation, surfactant dysregulation, and the overall change in the Lacriome.^4,5 One significant component of the Lacriome would also be to assess the cellular landscape of the diseased tissue and decode the cell to cell interactions in different phases of the disease. 

Based on pathological findings, the whole progression of PANDO consists of three distinct pathologic patterns^6,7: active chronic inflammation, proliferative chronic fibrosis, and complete lumina and subepithelial fibrosis. The first stage shows an active inflammatory pattern with signs of hypersecretion of goblet cells and subepithelial seromucous glands, with many lymphocytes infiltrating the epithelium and sub-epithelial tissues. The second stage shows transitional forms of chronic fibrosis, loss of intraepithelial goblet cells and epithelial cells to varying degrees, squamous metaplasia, and increased fibroblasts. The complete fibrotic stage was characterized by an obliterating mature fibrous tissue, with a total loss of goblet cells. It has been suggested that the development of PANDO may be regulated by bacteria, immune factors, anatomic factors, hormone, and cavernous body surrounding the nasolacrimal duct.^4,8–12 However, the changes at the cellular and molecular levels and the intercellular coordination and communications during the development of the disease remains unknown. 

Recently, single-cell RNA sequencing (scRNA-seq) analysis has created unprecedented opportunities to depict the transcriptome atlas of various tissues and average out underlying differences and interaction between cells. Various ocular tissues have been explored by scRNA-seq, including cornea,^13,14 sclera,¹⁵ retina,^16,17 and lacrimal gland,^18,19 but the lacrimal drainage system has never been explored before. Here, we aimed to establish the single-cell landscape of lacrimal sac cells with patients’ samples suffering from PANDO with different stages, and to further characterize the important changes and role of the epithelial cells and other cell types in the development of the disease. The present study characterized epithelial cells, T and B-lymphocytes, and neutrophils which are closely associated with the inflammatory response. Subsequently, the analysis was centered around epithelial and other cellular interactions to decode the cell-to-cell communication during the disease process. 

Materials and Methods

Ethical Approval

This study followed the tenets of the Declaration of Helsinki and was approved by the Ethics Committee of Shanghai Ninth People's Hospital. Written informed consents were obtained from all the patients. 

Patient and Specimen Details

A part of the lacrimal sac flap was obtained from five patients (all female patients, age range = 35–60 years) during the procedure of endoscopic dacryocystorhinostomy (EnDCR). Two patients were in the early inflammatory phase of the disease (<6 months of epiphora onset with dacryoendoscopy showing minimal edema with pinkish mucosa of the lacrimal sac) and three patients were classified as having late-phase inflammation (>6 months of onset, mucopurulent discharge, flaky hemorrhages of the lacrimal sac, submucosal scarring, and intra-sac synechiae). Immediately after its removal, the lacrimal sac tissues were kept in a tissue storage solution until further processing. 

Single-Cell Dissociation

Single-cell RNA-seq experiment was performed in the laboratory of Novel Bio Bio-Pharm Technology Co., Ltd. (Shanghai, China). The tissue samples were first washed with phosphate-buffered saline (PBS), minced into small pieces (approximately 1 mm³) on ice, and enzymatically digested with 0.1% Collagenase II (GIBCO, Cat #17101-015) and DNase I (Sigma-Aldrich, Cat #DN25) for 30 minutes at 37°C with gentle shaking. Then, the samples were sieved through a 70 µm cell strainer, and centrifuged at 300 g for 5 minutes. After the supernatant was removed, the pelleted cells were suspended in the red blood cell lysis buffer (Miltenyi Biotec, Cat #130-094-183) to lyse red blood cells and washed with PBS containing 0.04% BSA. The cell pellets were re-suspended in PBS containing 0.04% BSA and re-filtered through a 35 µm cell strainer. Dissociated single cells were then stained with AO/PI for viability assessment using Countstar Fluorescence Cell Analyzer. 

Single-Cell RNA Sequencing

The BD Rhapsody system was used to capture the transcriptomic information of the (5 sample-derived) single cells. Single-cell capture was achieved by random distribution of a single-cell suspension across >200,000 microwells through a limited dilution approach. Beads with oligonucleotide barcodes were added to saturation so that a bead was paired with a cell in a microwell. The cells were lysed in the microwell to hybridize mRNA molecules to barcoded capture oligos on the beads. Beads were collected into a single tube for reverse transcription and ExoI digestion. Upon cDNA synthesis, each cDNA molecule was tagged on the 5′ end (that is, the 3′ end of an mRNA transcript) with a unique molecular identifier (UMI) and cell barcode indicating its cell of origin. Whole transcriptome libraries were prepared using the BD Rhapsody single-cell whole-transcriptome amplification (WTA) workflow including random priming and extension (RPE), the RPE amplification PCR, and WTA index PCR. The libraries were quantified using a High Sensitivity DNA chip (Agilent) on a Bioanalyzer 2200 and the Qubit High Sensitivity DNA assay (Thermo Fisher Scientific). Sequencing was performed by Novaseq 6000 sequencer on a 150 bp paired-end run. 

Single-Cell RNA Statistical Analysis

The scRNA-seq data analysis was performed at the NovelBio Bio-Pharm Technology Co., Ltd. (Shanghai, China) on the NovelBrain Cloud Analysis Platform. Reads containing adapter and low-quality reads (q < 20) were trimmed using fasq with default parameter. UMI tools were applied for single cell transcriptome analysis to identify the cell barcode whitelist. The UMI-based clean data were mapped to human genome (Ensemble version 104) utilizing STAR mapping with customized parameter from UMI tools standard pipeline to obtain the UMIs counts of each sample. After that, we used the R (version 4.2.2) and Seurat package²⁰ (version 4.3.0) to load the gene count matrices and created a Seurat object for each sample with parameter min.cell 3 and then “merge” function was used to combine data as one Seurat object. Cells with detected genes above 200 and below 4000, the percentage of mitochondrial genes below 40 and the percentage hemoglobin genes below 10, were used for further analysis. Next, doublet cells were removed, which were detected by R package DoubletFinder²¹ (version 2.0.3). Mitochondria genes and hemoglobin genes were removed from the expression table. 

The filtered data were then normalized and scaled. The principal component analysis (PCA) was constructed based on the scaled data with the top 2000 high variable genes. To prevent batch effects from distorting subsequent data analysis, we used R package Harmony²² (version 0.1.1) for performing integration of each single cell datasets and the top 38 Harmony were used for Nearest-neighbor graph and UMAP construction. We performed the “FindClusters” function (resolution = 0.2, 0.5, 0.3, 0.2, and 0.8 for total cells, T cells, B cells, epithelial cell, and neutrophils, respectively) to cluster cells using the Louvain algorithm. 

Gene Ontology Enrichment Analysis

“FindAllMarkers” or “FindMarkers” function (test.use=“wilcox”, min.pct = 0.25, logfc.threshold = 0.25) was used to identify differentially expressed genes (DEGs) based on normalized data. Then, the DEGs with adjusted P values below 0.01 and at least 1.2-fold difference were used for Gene function analysis using R package “clusterProfiler”²³ (version 4.6.2). The selected GO terms among the top 100 ones were graphed with ggplot2 package.²⁴ 

Cell Communication Analysis

To enable a systematic analysis of cell–cell communication molecules, we proceeded with a cell-cell communication analysis between epithelial cells and other cell types in two periods by CellChat (version 1.6.1) with “CellChatDB.human,” an R package inference and analysis of cell-cell communication from single-cell data. Then, the cell-cell communications between the late inflammatory stage and early inflammatory stage were compared including interaction strength and dysfunctional signaling ligand-receptor pairs based on the differential gene expression analysis with 0.1 log fold change of ligands. 

Immunohistochemical Staining

All lacrimal sac tissues were acquired from DCR and fixed with 4% paraformaldehyde for 24 hours. After being dehydrated, the specimens were embedded in paraffin. Then, the 4 µm serial sections were cut and blocked with 3% BSA, followed by incubating with P16 Antibody (Abcam, Cat #ab108349), P27 KIP1 (Abcam, Cat #ab32034) overnight at 4°C. After that, the sections were washed by PBS 3 times and incubated with secondary antibody for 1 hour at room temperature. Diaminobenzidine was then used to detect the immunoreactivity and the sections were washed with flowing water. Finally, after hematoxylin re-staining and washing, the sections were observed under a microscope. 

Quantitative Measurement of Tear Cytokines

Tears of all the patients were collected before surgery by capillary tubes. Then, 20 µl of tears were obtained in diseased eyes and then preserved at −80°C before the test. Tear samples were diluted with diluent (RayBiotech, CAT #QA-SDB). Then, they added to the sample diluent into each well and incubated to block slides for 30 minutes at room temperature. After that, standard protein cocktail and biotinylated antibody cocktail were incubated for 1 to 2 hours, respectively, at room temperature. Cy3 equivalent dye-conjugated streptavidin was added after briefly spinning down and incubate at room temperature for 1 hour. Fluorescence detection was the last step and then we analyzed the data. 

Quantitative RT-PCR

For lacrimal sac tissues, total RNA was isolated by TRIzol reagent (Takara, CAT #9108) and then converted to DNA by PrimeScript RT reagent Kit (Takara, CAT #RR037Q). Quantitative PCR was performed using LightCycler 480 II (Roche, Switzerland) by UNICON qPCR SYBR Green Master Mix (Yeasen, CAT #11199ES03). Gapdh was used as endogenous control and gene expression level was evaluated by the 2 ^−∆∆Ct method. The primers used in this study are as follows: 

 

 

 

 

Results

Single-Cell Analysis of Lacrimal Sac in PANDO and Cell Type Identification

To generate a comprehensive single-cell atlas of the human lacrimal sac with PANDO, samples were obtained from five patients (Supplementary Figs. S1A, S1B). After preparing single-cell suspensions and removing dead cells, we performed scRNA-seq on these samples, resulting in transcriptome data from 25,791 cells (Fig. 1A, Supplementary Fig. S2A). Unsupervised Uniform Manifold Approximation and Projection (UMAP)-clustering analysis revealed the presence of 11 distinct cell types based on canonical lineage markers (Figs. 1B, 1C): T cells (CD3D⁺ and CD3E⁺), B cells (CD79A⁺, MS4A1⁺, and CD19⁺), plasma cells (CD79A⁺ and IGHG1⁺), monocytes (LYZ⁺ and CSG1R⁺), neutrophils (CDF3R⁺ and FCGR3B⁺), plasmacytoid dendritic cells (GZMB⁺), mast cells (KIT⁺ and CPA⁺), smooth muscle cells (MYH11⁺ and ACTA2⁺), fibroblast cells (COL1A1⁺), endothelial cells (PECAM1⁺ and ACKR1⁺), and epithelial cells (KRT5⁺, KRT8⁺, and EPCAM⁺). T cells, B cells, and neutrophils accounting for over 50% of all cells, increased in PANDO suggested their pivotal immune role in the pathogenic process (Figs. 1D-1F). Mast cells, smooth muscle cells, endothelial cells, and fibroblasts were decreased during the late inflammatory stage, and this was not unexpected. Epithelial cells, which are the first line of defense of the lacrimal mucosa, showed a decline in the late inflammatory phase (Supplementary Figs. S2B, S2C). 

Figure 1.

 

Single-cell RNA sequencing of lacrimal sacs cellular lineages reveal aberrant immune cell landscape in patients with PANDO. (A) Experimental design for single-cell sequencing of lacrimal sacs cells from patients with PANDO in the early inflammatory and in the late inflammatory stage. (B) UMAP plot showing the distribution of 25,791 cells from the early inflammatory (n = 2) and the late inflammatory stage (n = 3) lacrimal sacs. (C) Dot plot showing the expression of canonical markers for major cell types in the PANDO. (D) Stacked bar chart showing the relative proportion of major cell types in lacrimal sac from the early inflammatory stage and the late inflammatory stage patients with PANDO derived from scRNA-seq data. (E) Bar plot of the percentage change of T and B cells in two periods. (F) Density map showing the change of T cells, B cells, and neutrophils in two periods.

T Cells With Great Variation in Cell Number Actively Participate in Inflammatory Response in PANDO

T cells are the mainstay of the body's resistance to infections and tumors and their biological functions cannot be ignored, such as directly killing target cells, and assisting, or inhibiting antibody production by B cells. In the present study, the proportion of NK and T cells (NK&TC) in the early inflammatory stage was 24.76%, increasing to 44.23% during the late inflammatory stage (see Fig. 1E). We detected 9366 NK&TC as the largest cell population, which were divided into 8 clusters (Fig. 2A), including CD8⁺ memory T (Tm) cells (CD8⁺ and IL7R⁺), CD4⁺ central memory T (Tcm) cells (CD4⁺ and CCR7⁺), CD4⁺ effector memory T (Tem) cells (CD4⁺, IL7R⁺, and CCR7⁻), CD4⁺ exhausted T (Tex) cells (CD4⁺, PDCD1⁺, and TOX⁺), Treg cells (CD4⁺ and FOXP3⁺), proliferative T cells (CDC25A⁺ and CCNA2⁺), NKT cells (FCGR3A⁺ and GZMB⁺), and NK cells (NKG7⁺; Fig. 2B). CD4⁺ and CD8⁺ T cells were detected to increase in immunostaining experiments (Fig. 2C). Differential expressions of marker genes were used to annotate the various cell types and states. CD4⁺ Tcm exhibited the most significant increase in the late inflammatory groups (Supplementary Fig. S3A). This type of cell increases the immune response to re-stimulation by antigen and can proliferate and differentiate into effector T cells once stimulated by antigen.²⁵ 

Figure 2.

 

Characteristics of NK&TC subclusters in PANDO. (A) UMAP plot showing subsets of NK&TC. (B) Violin plot showing the expression of discriminative gene sets for subclusters in NK&TC. (C) Immunostaining for CD3, CD4, and CD8 in the lacrimal sac mucosa during two stages. (D) Volcano plot showing fold change and P value adjusted for the comparison of gene expression in the late inflammatory stage versus the early inflammatory stages in CD4⁺ Tcm cells. The genes of interest are indicated in volcano plots. (E) Representative GO terms enriched by enhanced DEGs of CD4⁺ Tcm, CD4⁺ Tem, and CD8⁺ Tm.

We then explored the functional changes of T cells and their subsets in PANDO by conducting a DEG analysis and Gene Ontology (GO) analysis. According to the top 10 DEGs of T cells in the late inflammatory/early inflammatory stages comparison, genes of IGHG1, IGHG3, IGHG4, IGLC3, and HLA-DRB1 were highly expressed, which were related to antigen receptor-mediated signaling pathway, defense response to bacterium and humoral immune response (Supplementary Figs. S3B, S3C, Supplementary Tables S1, S2). At the late inflammatory stage in PANDO, CD4⁺ Tcm exhibited upregulated IGHG, FCER1G, IGHA1, IGHA2, and IL6R genes, indicating the function of defense response to bacterium (Fig. 2D, see Supplementary Table S1). Adaptive immune responses were activated by CD4⁺ Tcm and CD8⁺ Tm, such as antigen receptor-mediated signaling pathway, B cell mediated immunity, and the immunoglobulin mediated immune response (Fig. 2E, see Supplementary Table S2). 

B Cells Regulated Migration and Proliferation of Epithelial Cells During Late Inflammatory Stage

The B cells were observed to be one of the main cell types affected in patients with PANDO. There was a significant increase in the proportion of B cells from 3.17% to 21.09% during the pathogenesis of PANDO (see Fig. 1E), which was verified in immunohistochemical results (Fig. 3A). To gain a deeper understanding of the heterogeneity of B cells and their immunological properties, we conducted a more detailed clustering analysis of B cells. This analysis led to the identification of five subsets of MS4A1⁺ B cells (Figs. 3B, 3C). Four clusters were categorized as memory B cells, and genes TMX4, FOSB, EGR, and RGS were highly expressed, respectively. One cluster was characterized by the expression of IGHD and IGHM, pointing to a B follicular cell cluster. Plasma cells (PCs), the terminally differentiated, non-dividing effector cells of the B-cell lineage, were also subjected to further subclustering analyses, which revealed two subpopulations (Figs. 3D, 3E, Supplementary Fig. S4A): IgA⁺ PCs (IGHA2⁺) and IgG⁺ PCs (IGHG1⁺). 

Figure 3.

 

Characteristics of B cells subclusters in PANDO. (A) Immunostaining for CD20 in the lacrimal sac mucosa during 2 stages. (B) UMAP plot showing subsets of MS4A1⁺ B cells. (C) Violin plot showing the expression of discriminative gene sets for subclusters in MS4A1⁺ B cells. (D) UMAP plot showing subsets of plasma cells. (E) Dot plot showing the expression of discriminative gene sets for subclusters in plasma cells. (F) Representative GO terms enriched by down-regulated DEGs of TMX4⁺ Bm cells and follicular B cells. (G, H) Volcano plot showing fold change and P value adjusted for the comparison of gene expression in the inflammatory stage versus early inflammatory stages in TMX4⁺ Bm cell (G) and follicular B cells (H).

An increase in subpopulations of TMX4⁺ and RGS⁺ memory B cells and follicular B cells was observed during the late inflammatory stage. Subsequently, we compared the gene expression in the late inflammatory stage versus the early inflammatory stages in TMX4⁺ memory B cells and follicular B. Interestingly, in TMX4⁺ memory B cells, with the biggest change of cell proportion (Supplementary Fig. S4B), genes related to epithelial cell proliferation (CCL2, EGR3, HMOX1, IGFBP5, TNF, and ZFP36) were downregulated specifically during the late inflammatory stage, suggesting a potential effect of memory B cells on epithelial cell proliferation (see Figs. 3F, 3G, see Supplementary Tables S3, S4). In the follicular B cells, the downregulated genes were associated with epithelial cell migration (DCN, DUSP10, EGR3, HMOX1, SERPINF1, and SOX9; see Figs. 3F, 3H; Supplementary Tables S3, S4). These results suggest that epithelial migration and proliferation during the late inflammatory stage was influenced by MS4A1⁺ B cells, at least to some extent. 

Epithelial Cells Expressed SASP and Upregulated Inflammatory Genes During the Late Inflammatory Stage

Epithelial cells are a vital group of parenchymal cells covered with a mucin layer possessing adhesive properties and antimicrobial peptides that play a crucial role in information exchange and transmission.^26,27 They are one of the first line of structural defenses against the microorganisms and antigens arriving in the lacrimal duct lumen. Hence, it is necessary to focus on the changes in epithelial cells. The present study identified four clusters of epithelial cells (Figs. 4A, 4B). One cluster exhibited high expression of secretory protein-related genes such as TFF3, TFF1, MUC5AC, and MUC5B, protecting mucosa from pathogen invasion and injury.^27–29 As a result, we defined this cluster as the “secretory” type, with a significant increase in cell number during the late inflammatory stage (Fig. 4C), which corroborates with the clinical phenomenon that the lacrimal sac epithelium hypersecretes during acute inflammations. Cluster 2 demonstrated higher levels of MALAT1 expression compared to other clusters, leading us to label them as “MALAT1⁺ epithelial cells.” MALAT1, a non-coding RNA, was first identified in non-small cell lung cancer. It has been reported that MALAT1 is related to cancer cell proliferation, inflammation, metastasis, fibroblast proliferation, and endothelial cell proliferation.³⁰ MALAT1 aggravates inflammation through modulating NF-κB and regulate airway epithelial cells to secrete IL-8 in pneumonia.³¹ In our study, genes of NFKB1 and CXCL8 were expressed in MALAT1⁺ epithelial cells (Supplementary Fig. S5A), but whether they were regulated by MALAT1 needed to be further researched. It was also noted that the cluster 3, whose ratio at 2 phases had no significant change, highly expressed the basal cell markers KRT5 and KRT8 as well as the progenitor marker S100A2,¹³ hence this type was labeled as basal epithelial cells. Basal epithelial cells are located adjacent to the basement membrane and have capacity of proliferation associated with tissue stem cells. Additionally, cluster 4 (CDH19⁺ epithelial cells) expressed S100B, CDH19, and PLP1, which are typically associated with myelination or growth regulation, and downregulated in the late inflammatory stage (Supplementary Fig. S5B). This could be due to the presence of nerve cells within the epithelial cell population, influencing their biological functions, or because of interactions between epithelial cells and nerves. 

Figure 4.

 

The expression of SASP in epithelial cells in PANDO. (A) UMAP plot showing subsets of epithelial cells. (B) The expression levels of selected genes in four subclusters of epithelial cells. (C) Stacked bar chart showing the relative proportion of subsets of epithelial cells from the early inflammatory and late inflammatory stages in PANDO. (D, E) The expression of genes related to SASP and cell cycle arrest in epithelial cells. (F) Quantitative real-time PCR analysis for CDKN2A, CDKN1B, and TP53. (G) Immunostaining of lacrimal sac epithelial cells for P27^KIP1(CDKN1B) and P16^INK4(CDKN2A). (H) Tear protein microarray model. (I) Heatmap of enhanced SASP-related cytokines in tear from 12 patients with PANDO during the early inflammatory stage and 13 during the late inflammatory stage.

The present study also observed that epithelial cells expressed a wide range of chemokines, particularly SASP (Fig. 4D, Supplementary Fig. S5C). Secretory epithelial cells in the late inflammatory phase enhanced the expression of IL-8, TNF, CXCL1, CXCL3, MMP2, and MMP10, and basal cells also upregulated various SASP. This finding suggested that epithelial cells of lacrimal sac possessed characteristics associated with SASP. Subsequently, we conducted further analysis on genes related to cell cycle arrest. It was found that gene P27^KIP1(CDKN1B) and P16^INK4(CDKN2A) were significantly upregulated during the late inflammatory stage (Fig. 4E). These genes are known to be commonly associated with senescence and hold great importance.^32,33 Moreover, AKT1 and AKT2, related to premature senescence and shortened cellular lifespan, exhibited similar upregulation during the late inflammatory phase.^34,35 To support these findings, we performed immunohistochemical staining and PCR on lacrimal sac samples obtained from patients with PANDO (Figs. 4F, 4G). Besides, we detected various chemokines in the tears from patients with PANDO during two stages. We first compared 21 chemokines between patients with PANDO and healthy people and we found that IL-6, TNF, IL-1β, and IL-10 were upregulated in tears. Then, we divided the patients into two groups, the early inflammation group and the late inflammation group, and the expression of IL-6, TNF, and IL-1β were apparently enhanced, whereas IL-10 had no difference between the two stages (Figs. 4H, 4I, Supplementary Table S5). Compared to the early inflammatory stage, the results demonstrated higher expression levels of SASP-related genes and cell cycle arrest genes (CDKN1B or CDKN2A) in the epithelial cells during the late inflammation stage. The concept of SASP in epithelial cells of PANDO has never been addressed in previous studies. It was noted in the present study that as the inflammation aggravated, aging genes were activated and epithelial cells of lacrimal sac secreted more SASP. 

To assess the biological function of epithelial cells during the inflammatory stage, DEG and GO analysis was conducted between the two stages of PANDO in different epithelial subtypes (Figs. 5A–5C, Supplementary Tables S6, S7). We observed an upregulation of mucin-related genes (MUC5AC, MUC5B, TFF1, TFF3, FCGBP, and ZG16) in secretory epithelial cells (see Figs. 5A, 5D). The production of mucus by epithelial cells was partly attributed to the immune response.³⁶ TFF-peptides like TFF1 and TFF3, were reported to play a role in anti-apoptotic and mitogenic processes as well as epithelial regeneration.^27,37 FCGBP, one of the early response proteins after microbial infection, is reported to be secreted by goblet cells and classified as mucin-like protein.^38–40 ZG16 is a small lectin-like protein that binds to Gram⁺ bacteria, preventing their proximity to the epithelium.^41,42 Therefore, these proteins form a mucus layer covering the wall of lacrimal duct. Secretory epithelial cells also exhibited high expression of genes, such as SAA1, SAA2, C3, S100A8, and S100A9 (Fig. 5E), which are associated with inflammation during the late inflammatory stage. Antibacterial-related genes including PI3, SLPI, and BPIFA1 were also detected (Fig. 5F). Secretory epithelial cells were predominantly involved in inflammation and anti-bacterial pathways. This suggested that secretory epithelial cells actively responded to inflammation and pathogen invasion during the late inflammatory stage. GO terms of defense response to bacterium, antimicrobial humoral response, neutrophil chemotaxis, neutrophil migration, and leukocyte chemotaxis were enriched in secretory epithelial cells (Fig. 5G, Supplementary Table S7). As shown in Figures 5A–5C, 5E, and Supplementary Figures S5D, as secretory epithelial cells, basal epithelial cells, and MALAT1⁺ epithelial cells also exhibited similar genes (SAA1, S100A8, S100A9, and C3). Subsequent GO analysis also exhibited enriched pathways related to acute inflammatory response and acute-phase response and other immune response pathways in basal epithelial cells, whereas immune response and inflammatory response-related pathways were also upregulated in MALAT1⁺ epithelial cells. These findings suggested that the expression of inflammatory genes was enhanced in the three subclusters mentioned above. 

Figure 5.

 

Characteristics of epithelial cells subclusters in PANDO. (A–C) Volcano plots showing fold change and P value for the comparison of gene expression in the late inflammatory stage versus the early inflammatory stages in secretory epithelial cells (A), MALAT1⁺ epithelial cells (B), and basal epithelial cells subtypes (C). (D) Split violin plots showing mucin-related genes (MUC5AC, MUC5B, TFF1, TFF3, FCGBP, and ZG16). (E) Split violin plots showing inflammatory genes in secretory epithelial cells, MALAT1⁺ epithelial cells and basal epithelial cells subtypes. (F) Split violin plots showing antibacterial-related genes (PI3, SLPI, and BPIFA1). (G) Representative GO terms enriched by enhanced DEGs of secretory epithelial cells, MALAT1⁺ epithelial cells, and basal epithelial cells.

Neutrophils are Recruited by Epithelial Cells Through CXCLs-CXCRs

To elucidate the impact of specific cell populations on epithelial cells, as well as the reciprocal effects on those cell populations, the cellular interaction network was investigated based on immune cell clusters (Figs. 6A, 6B). The present analysis revealed that the highest number of interactions occurred between epithelial cells and neutrophils, followed by interactions between epithelial cells and monocytes, with the epithelial cells exhibiting predominant effects on both. To find which subcluster of neutrophils interacted with epithelial cells, neutrophils were categorized into four distinct clusters (Figs. 6C, 6D). Cluster 1 exhibited TNFAIP6 expression, whereas cluster 3 showed high levels of S100A4, and cluster 4 expressed GBP1. However, cluster 2 did not demonstrate any characteristic gene expression but highly expressed IL1B. We found the interactions of chemokines (CXCL1, CXCL3, CXCL5, CXCL6, and CXCL8) and their receptor (CXCR1 and CXCR2) between epithelial cells and neutrophils were increased during the late-inflammatory stage (Fig. 6E, Supplementary Fig. S6A). CXCLs have been reported to be secreted by epithelial cells, including CXCL8 secreted by lung epithelium,⁴³ CXCL1/2 secreted by intestinal epithelium,⁴⁴ aggravating inflammation by recruiting neutrophils via CXCR1 or CXCR2 of neutrophils. Besides, CXCLs/CXCRs signaling axis was reported to be relevant to cellular senescence. Mitochondrial dysfunction in renal tubular epithelial cells through activating P16^INK4 and β-catenin by CXCR2.⁴⁵ Not only epithelial cells, but neutrophils also can secrete CXCLs to recruit more neutrophils. The interactions involving CXCL1-CXCR1/2 and CXCL8-CXCR1/2 were particularly pronounced in neutrophils during the late inflammatory stage (Fig. 6F), which suggested that they might recruit additional neutrophils through these pathways after being recruited to the inflammatory sites by epithelial cells. These results demonstrate that neutrophils may actively participate in the inflammatory response in PANDO in the form of chemotaxis by CXCLs secreted by epithelial cells and themselves to produce a cascade amplification effect. 

Figure 6.

 

Characteristics of neutrophils subclusters in PANDO. (A, B) Visualized network graphs of number of interactions showing the cell-type-specific cell-cell interactions among the early inflammatory stage (A) and the late inflammatory stage (B). (C) UMAP plot showing subsets of neutrophils. (D) Violin plot showing the expression of discriminative gene sets for four subsets in neutrophils. (E) The enhanced interaction between various cell types in individuals with PANDO during the late inflammatory stage. (F) The interaction between neutrophils in individuals with PANDO.

Furthermore, the interaction involving CSF3-CSF3R, associated with proliferation of neutrophils and essential for granulocytic maturation was detected. It was also found that the expression of SAA1-FPR2 was upregulated during the late inflammatory stage. It was proven that SAA1 could be chemotactic for leukocytes via activation of FPR2.⁴⁶ 

The Network Between Epithelial Cells and Other Cell Types

Next, the cell-cell communication between epithelial cells and other cell types in PANDO were analyzed (Figs. 7A, 7B). We discovered an upregulated MDK-NCL signaling pathway between epithelial cells and T cells during the late-inflammatory stage (Fig. 7C). MDK, a heparin-binding growth factor, is known to regulate various processes including inflammatory response, cell proliferation, cell adhesion, cell growth, cell survival, cell differentiation, and cell migration.^47,48 NCL has been reported to interact with multiple ligands involved in cell proliferation and apoptosis at the cell membrane.⁴⁹ We observed MDK expression in epithelial cells, whereas NCL exhibited high expression in CD4⁺ Tcm, CD4⁺ Tem, and CD8⁺ Tm cells. This finding suggests that epithelial cells may be involved in recruiting T cells through the secretion of MDK. Furthermore, we found that the interaction of MDK-NCL was the main pathway between epithelial cells and B memory cells, and this interaction was enhanced during the late inflammatory stage (Figs. 7D, Supplementary Fig. S7A). These findings further supported the notion that epithelial cells may recruit B memory cells to the inflammatory sites. 

Figure 7.

 

An epithelial cell-based cell-cell communication network in PANDO. (A, B) Network graphs of number (A) and weights (B) of interactions showing the interactions among all cell types. (C) The interaction of epithelial cells with NK&TC via MDK-NCL in patients with PANDO of two stages. (D) The interaction of epithelial cells with BCs via MDK-NCL in patients with PANDO of two stage. (E) The interaction of epithelial cells with fibroblasts, smooth muscle cells, monocytes, and endothelial cells in patients with PANDO of two stages.

As fibroblasts, smooth muscle cells, and endothelial cells with large changes in cell proportions between the two periods, there is a non-negligible interplay between them and epithelial cells (Fig. 7E). CEACAM1-SELE was the most significant ligand-receptor between epithelial cells and endothelial cells. CEACAM1 was proven as a pivotal factor in inflammation,⁵⁰ whereas SELE plays the role of immuno-adhesion in endothelial cells. We found that the expression of platelet-derived growth factor (PDGF) in epithelial cells, acting on fibroblasts and promoting their proliferation,⁵¹ was enhanced. 

In conclusion, given that epithelial cells are the first to be influenced by pathogen in the lacrimal lumen, they take a pivotal role in PANDO. We speculate that epithelial cells recruited neutrophils via SASP on the one hand and directly enhanced the inflammatory response via upregulating inflammatory genes on the other (Fig. 8). In addition, T cells and B cells, which showed the greatest change in number, were also recruited by epithelial cells. The increase in immune cells undoubtedly worsened the inflammatory environment. 

Figure 8.

 

The hypothesis that epithelial cells influence other cell types and finally lead to inflammation.

Discussion

In this study, the first comprehensive single-cell atlas aiming to elucidate the cellular diversity during the early inflammatory and late inflammatory stages of PANDO was constructed. Our analysis identified 11 major cell types, each exhibiting distinct functions in the context of PANDO. To the best of the authors’ knowledge, this is the first study to provide a comprehensive overview of the dynamic cellular profile of PANDO. This is one more step at unravelling yet another dimension of Lacriome looking specifically at the molecular and cellular origins of inflammation in PANDO. 

Although we were unable to conduct a detailed functional analysis of each cell type, there were several noteworthy findings. First, it was observed that epithelial cells highly expressed cytokines, such as IL-8, TNF, CXCL1, CXCL3, MMP2, and MMP10, all of which belonged to SASP. Further analysis found that the expression of senescence genes was also upregulated in the late inflammatory period. We additionally took the lacrimal sacs of patients in both periods for immunohistochemistry and verified the above results, suggesting that epithelial senescence might be a possible mechanism for the development of PANDO. It was the first time to explore SASP in epithelial cells and its function in promoting immune imbalance of lacrimal duct, resulting in PANDO. We also found that neutrophils were recruited by epithelial cells and themselves through CXCLs/CXCRs signal axis. It might be evidence that SASP secreted by epithelial cells resulted in chronic inflammation. It is reported that senescence of airway epithelial cells leads to inflammation in the lungs and accumulation of leukocytes.^52,53 The senescent cells in skin triggers the elevated production of IL-6.⁵⁴ In most populations, inflammaging increases with age. The overall inflammatory state is associated not only with an increase in individual inflammatory markers, but also with complex interactions between various inflammatory mediators.⁵⁵ For another, inflammation would aggravate senescence.⁵⁶ Inhibiting IL-6 has been reported to improve erythroid progenitors function in aged mice.⁵⁷ Besides, a low-grade inflammatory environment can promote the senescence of neutrophils.⁵⁸ But whether SASP was the first step to induce a series of inflammation in the lacrimal sac was still unknown. Previous studies in terms of senescence elsewhere in the body focused more on fibrosis.⁵⁹ The interaction of epithelial cells and fibroblasts might be a potential point to reveal the progression of PANDO. 

Furthermore, it was noted that a significant increase in immune cells, particularly TCs and BCs, during the late inflammatory stage. The proportions of TCs and BCs rose from 24.76% and 3.17% in the early inflammatory phase to 44.23% and 21.09% in the late inflammatory phase respectively, accompanied by a moderate increase in neutrophils and monocytes. Epithelial cells of the lacrimal sac in this research accounted for only 2.81% in the early inflammatory stage and 2.41% in the late inflammatory stage, which corresponds with the previous study showing that epithelial cells make up a small share of cell numbers.³⁶ As inflammation in the lacrimal sac progresses, epithelial cells decrease even disappear. Interestingly, we noticed a slight decrease in the percentage of epithelial cells, although their function appeared to undergo changes. Consistent with previous studies,⁶ we observed an increased secretion of mucins by epithelial cells during this phase. Both mucins and trefoil factor (TFF) play crucial roles in protecting the lacrimal drainage mucosa, and the elevated levels observed during the late-inflammatory phase suggest the initiation of self-protective mechanisms in the lacrimal tissue. This response aims to reduce epithelial cell damage caused by foreign pathogens, possibly triggering their own apoptosis. Given their role as the first line of defense against mucosal invasion within the lacrimal duct, changes in epithelial cells during the progression of PANDO could potentially induce cascading alterations in other cell types. Our findings indicated the presence of multiple ligand-receptor relationships between epithelial cells and immune cells, suggesting potential pathways through which epithelial injury may trigger immune cell migration. The dysregulation of immune cells, particularly T cell subsets, may prove to be a critical factor in the development of PANDO.⁹ Previous work by Yang et al.⁹ has reported an imbalance between Th1 and Th2 cells, with Th1 cells prevailing in PANDO. Unfortunately, in our study, we did not isolate Th cells, and some marker genes appeared to be challenging to detect. 

The smaller sample size of patients in the present study is also one of its limitations. We collected two early and three late inflammatory samples in this research and there is no doubt that it would be better to include more to reduce the random errors brought by patient heterogeneity. We had conducted immunohistochemical staining and qRT-PCR, but it is still necessary to do cell and animal experiments in the future to verify our conjectures, such as further cell-cell interactions and epithelial cell senescence during the late inflammatory stage. This study is the starting point towards enhanced and next-level molecular understanding of the etiopathogenesis of PANDO. 

Acknowledgments

Supported by Shanghai Key Laboratory of Orbital Diseases and Ocular Oncology (2022SKLE-KFKT001), Shanghai Key Clinical Specialty, Shanghai Eye Disease Research Center (2022ZZ01003), Shanghai Sailing Program (21YF1423000) and National Natural Science Fund of China (82388101). The authors thank NovelBio Co., Ltd. for the support of bioinformatics analysis with their NovelBrain Cloud Analysis Platform (www.novelbrain.com). 

Data Availability Statements: The raw sequence data reported in this paper can be available from the GEO (https://www.ncbi.nlm.nih.gov/geo/) databases, under the accession number: GSE252058. The scripts used to generate all analyses can be accessed from https://github.com/augsdihh/PANDO.git. 

Author Contributions: C.X., H.S., M.J.A., and Y.X. designed the study. W.Z. and H.H. wrote the manuscript and analyzed the scRNA-seq data. X.L. acquired the samples. W.Z., L.Z., and L.L. performed the experiments. Y.D., M.J.A., and X.L. analyzed the data. All authors contributed to read and approved the submitted version. C.X., H.S., and Y.X. supervised the project. 

Disclosure: W. Zhang, None; H. Huang, None; X. Liu, None; L. Zhang, None; L. Li, None; Y. Ding, None; Y. Xiao, None; M.J. Ali, None; H. Sun, None; C. Xiao, None