Single-Nuclei Characterization of Lacrimal Gland in Scopolamine-Induced Dry Eye Disease

Abstract

Purpose: The lacrimal gland (LG) is the main organ responsible for tear secretion and an important pathogenic site for dry eye disease (DED). This study aimed to comprehensively characterize LG cellular heterogeneity under normal and DED conditions using single-nucleus RNA sequencing (snRNA-seq).

Methods: Single LG nuclei isolated from mice with or without DED induced by scopolamine (SCOP)/desiccating stress (DS) were subjected to snRNA-seq using the 10x Genomics platform. These cells were clustered and annotated using the t-distributed stochastic neighbor embedding (t-SNE) method and unbiased computational informatic analysis. Cluster identification and functional analysis were performed based on marker gene expression and bioinformatic data mining.

Results: The snRNA-seq analysis of 30,351 nuclei identified eight major cell types, with acinar cells (∼72.6%) being the most abundant cell type in the LG. Subclustering analysis revealed that the LG mainly contained two acinar cell subtypes, two ductal cell subclusters, three myoepithelial cell (MECs) subtypes, and four immunocyte subclusters. In the SCOP-induced DED model, three major LG parenchymal cell types were significantly altered, characterized by a reduced proportion of acinar cells with a lowered secretion potential and an augmented proportion of ductal cells and MECs. LG immunocytes in DED scenarios showed an intensified inflammatory response and dysregulated intercellular communication with three major LG parenchymal cells.

Conclusions: Overall, this study offers a systemic single-nucleus transcriptomic profile of LGs in both normal and DED conditions and an atlas of the complicated interactions of immunocytes with major LG parenchymal cells. The findings also facilitate understanding the pathogenesis of DED.

Dry eye disease (DED) is recognized as a multifactorial ocular disorder with a global incidence between 5% and 50% among adults.¹ DED is mainly classified into two distinct subtypes: aqueous deficient and evaporative dry eye.² According to the Tear Film and Ocular Surface Society Dry Eye Workshop II (DEWS II), DED is pathologically characterized by the loss of tear film stability and varying degrees of serious ocular discomfort, including eye dryness, burning, pain, and sensitivity to light.³ Emerging evidence has revealed that the pathogenesis of DED is closely associated with numerous causes, such as corneal nerve impairment, mucin dysregulation, meibomian gland dysfunction, and lacrimal secretion insufficiency.³ Moreover, various risk and lifestyle factors contribute to the development of DED, including mental, physical, and social health, as well as cosmetics and digital environments.⁴ Together, these causes form a vicious cycle leading to DED; however, the exact pathogenesis of DED remains unclear. 

Tear-fluid insufficiency is a major characteristic of DED. As an important tubular–acinar exocrine gland, the lacrimal gland (LG) is essential for maintaining ocular homeostasis through secretion of the aqueous components of the tear film.⁵ The secretion of tears is tightly controlled by the lacrimal functional unit, which is comprised of the ocular surface, LG, and nervous system.⁶ Growing evidence underscores the functional associations between LG damage and DED in multiple pathological conditions, including diabetes,⁷ systemic lupus erythematosus,⁸ and Sjögren's syndrome.⁹ Given the importance of LG in regulating tear secretion, the LG has recently received significant attention, including LG organoid^10,11 and DED¹² development. However, the systemic pathological alterations of LG under various DED conditions have not been fully characterized. 

To better understand the pathogenesis of DED, multiple representative mouse models have been established based on different pathogenic conditions to decipher its pathophysiology and underlying mechanisms. Widespread applied models include: (1) autoimmune and autoinflammatory disease models, such as Aire knockout mice, non-obese diabetic (NOD) mice, NOD.H-2b mice, and MRL/MPJ Fas^lpr/lpr mice¹³; (2) desiccating an environmental stress model, which can be used to investigate aqueous-deficient DED, such as scopolamine (SCOP) and desiccating stress (DS) dry eye^14,15; and (3) models based on the aging of mice applied for age-related chronic DED. Although these models are widely used, none of them can encompass all the symptoms of human DED. 

Although accumulating evidence has shown significant pathological alterations of the LG during DED,^16–18 there is still a lack of systemic comparison of the LG between normal and DED conditions. The techniques of single-cell RNA sequencing (scRNA-seq) and single-nucleus RNA sequencing (snRNA-seq) have made it available for comprehensively examining the complexity of human diseases, including the diversity of cell types, cellular heterogeneity, and intercellular communications.¹⁹ These two kinds of techniques have been widely applied to various diseases, including cancer,²⁰ autoimmune diseases,²¹ and eye diseases.²² SCOP/DS-induced DED is a well-established DED model using subcutaneous injection of SCOP in conjunction with a constant airflow chamber.^14,15 In this study, we comprehensively investigated the cellular and functional heterogeneity of the LG under normal and DED conditions using snRNA-seq. This study presents the systemic molecular profile of the LG under normal and DED conditions, aiding in subsequent pathogenesis and therapeutic studies. 

Methods

Animal Models

Adult C57BL/6 mice, 6 to 10 weeks of age and weighing 18 to 20 g, were procured from SPF Biotechnology Co., Ltd. (Beijing, China). The mice were housed in the animal center of the Shandong Eye Institute. Ethical approval for all experimental procedures was obtained from the Ethics Committee of the Shandong Eye Institute (Ethics number: SDSYKYJS No.20220413). All experimental procedures were conducted in accordance with the guidelines outlined in the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research. The dry eye mouse model was induced by subcutaneously administering SCOP hydrobromide (ab141080, 2.5 mg/mL in PBS; Abcam, Cambridge, UK) to the mice three times a day (at 9:00 AM, 3:00 PM, and 9:00 PM) and exposing them to an air-drying device continuously for 7 consecutive days.²³ The control group of mice received no treatment. The extraorbital LG was extracted after all experimental procedures. 

Tissue Processing

The freshly harvested LG was fixed in 4% paraformaldehyde for 4 hours and subsequently dehydrated in a 30% sucrose solution for 8 hours, then embedded in optimal cutting temperature (OCT) compound (Tissue-Tek, Torrance, CA, USA) and quickly frozen using dry ice. A cryostat was utilized to cut sections with a thickness of 7 µm, which were then placed on CITOGLAS slides (Citotest Scientific, Nanjing, China) and allowed to air dry overnight. These dried sections were subsequently stored at −80°C. 

Antibodies

The following antibodies were used for tissue immunohistochemistry: rabbit anti-MIST1/BHLHA15 (14896S; Cell Signaling Technologies, Danvers, MA, USA); rabbit anti-alpha smooth muscle actin (ab150301; Abcam); rabbit anti-annexin A1/ANXA1 (ab214486; Abcam); rabbit Sox6 polyclonal antibody (PA5-34616; Thermo Fisher Scientific, Waltham, MA, USA); rabbit anti-NDUBB11 polyclonal (YN0924; ImmunoWay Biotechnology, Newark, DE, USA); and mouse anti-FGFR-2 (ab10648; Abcam). 

Isolation of Nuclei From Mouse Tissue

For each experiment, LGs were obtained from C57BL/6 mice (age 8–10 weeks; weight 20–24 g). A modified DroNc-seq technique was used to separate nuclei.²⁴ Tissue was chopped to 1 mm³, homogenized with Nuclei EZ Lysis Buffer (NUC101; Sigma-Aldrich, St. Louis, MO, USA), and filtered through a 40- µm cell strainer. The pellet was then resuspended in PBS after centrifugation at 500g for 5 minutes at 4°C. Red fluorescent protein (RFP)⁺ nuclei were collected individually in a 384-well plate for Smart-seq2 or 10x Genomics sequencing after resuspension in the MoFlo Astrios EQ Cell Sorter (Beckman Coulter, Brea, CA, USA). 

Single-Nucleus RNA-Sequencing

A total of 12,000 nuclei were placed into a 10x Genomics Chromium controller (PN110203; 10x Genomics, Pleasanton, CA, USA) immediately after isolation. The Chromium Single Cell 3′ Library & Gel Bead Kit v2, 16 rxns (PN-120237; 10x Genomics); Chromium Single Cell A Chip Kit, 48 rxns (PN-120236; 10x Genomics); and Chromium i7 Multiplex Kit, 96 rxns (PN-120262; 10x Genomics) libraries were prepared according to the manufacturer's instructions for Experiments 1 and 2. Chromium Single Cell 3 GEM, Library & Gel Bead Kit v3, 16 rxns (PN-1000075; 10x Genomics); Chromium Single Cell B Chip Kit, 48 rxns (PN-1000073; 10x Genomics); and Chromium i7 Multiplex Kit, 96 rxns (PN-120262; 10x Genomics) were used for the third experiment. Libraries were sequenced using an Illumina NovaSeq 6000 Sequencing System (20012850; 10x Genomics) with 50,000 raw reads per nucleus. 

Sample Identity Distribution

Sample identification and distribution of integrated LG cells from nine control mice and three SCOP mice were performed using the souporcell algorithm, as previously described.²⁵ 

Data Processing and Downstream Analysis

Transcriptome data were aligned to the relevant reference genome using the 10x Genomics CellRanger 3.1.0 pipeline. Subsequently, read count matrices were generated for each sample using the CellRanger count. Upon importing the data into the Seurat R 3.2.0 package,²⁶ quality control measures were implemented involving the exclusion of cells with fewer than 500 detected genes or those with a mitochondrial gene ratio exceeding 10%. Genes expressed in fewer than five cells were filtered out. The DoubletFinder 2.0.2 package²⁷ was used to identify potential doublets by adjusting the bimodality coefficient and number of artificial doublets to 0.25. Concomitantly, cells expressing previously documented dissociation-induced gene signatures were progressively eliminated, unless alternative marker genes were expressed. The subsequent data analysis utilized the canonical correlation analysis method to mitigate batch effects across distinct experimental batches.²⁶ To normalize the data, the LogNormalize method was applied to normalize the data, accounting for inherent variation caused by mitochondrial gene expression. Principal component analysis was performed on highly variable genes to facilitate cell clustering. Utilizing the graph-based shared nearest neighbor method, clustering was executed at a resolution of 0.25 across the top 13 principal components. This process yields six unsupervised cell clusters. Clustering outcomes for individual or aggregated samples were visually represented using the Uniform Manifold Approximation and Projection (UMAP) method or t-distributed stochastic neighbor embedding (t-SNE). Cell types were determined by conducting differential expression analysis and identifying cluster-specific marker genes using the Seurat FindMarker function. 

Cell–Cell Communication Analysis

Systemic analysis of cell–cell communication based on ligand–receptor interactions was performed using CellPhoneDB 1.1.0 (https://github.com/Teichlab/cellphonedb).²⁸ Subsequently, we examined the receptors or ligands expressed in at least 10% of cells within a specific cell type, with a significance threshold of P < 0.05. For the selected significant receptor–ligand pairs, we utilized the R packages ggalluvial and circlize (R Foundation for Statistical Computing, Vienna, Austria) to create visual representations of the interaction links. 

Gene Enrichment Analysis

We conducted enrichment analysis of the significantly differentially expressed genes using ToppGene. We used a combined score to determine the ranking of pathway enrichment, ontologies, transcription factor networks, and protein networks. This combined score was calculated by taking the logarithm of the P value obtained from the Fisher's exact test and multiplying it by the z-score representing the deviation from the expected rank. 

Results

Single-Nuclei Transcriptomic Atlas of LGs From Control and SCOP/DS-Induced DED Mice

To comprehensively characterize the cellular diversity and molecular profiles of the LG under normal and DED conditions, LG samples from naïve control (Ctrl) and DED (SCOP) mice induced by SCOP/DS were collected for snRNA-seq (Fig. 1A). After quality control and removal of batch effects, 30,351 nuclei with an average of 1724 genes and 5525 transcripts per cell were obtained for downstream analyses (Supplementary Figs. S1A, S1B). 

Figure 1.

 

The snRNA-seq analysis of LGs from Ctrl and SCOP-induced DED mice. (A) Method flowchart. (B) The t-SNE clustering of mouse LG cells with colors based on eight distinct clusters. (C) Proportions of different cell types in normal LGs (top) and in SCOP and Ctrl groups (bottom). (D) Dot plot showing the expression of representative genes for each cell type. (E) The immunofluorescence of Bhlha15, Anxa1, and Acta2. (F) Heatmap for DEGs in each cell type (left) and representative gene ontology terms of all DEGs in each cell type (right).

We identified 11 clusters based on their specific transcriptional profiles (Supplementary Fig. S1C), with comparable distributions of all clusters across samples from different groups (Supplementary Fig. S1D). These clusters were classified into eight cellular types according to their specific cell markers: acinar cells, ductal cells, myoepithelial cells (MECs), mesenchymal cells, immune cells, endothelial cells, Kit⁺ cells, and a cluster of mitotic cells (Fig. 1B). Among these clusters, acinar cells (∼72.6%) were identified as the most prominent cell type in the LG based on the classic acinar cell markers Bhlha15²⁹ and Scgb2a2³⁰ (Figs. 1C, 1D). A cluster termed ductal cell (∼7.68%) was annotated according to the specific expression of Sftpd,³¹ Wfdc2,³² and Anxa1.³³ The clusters corresponding to MECs (∼3.74%) were identified by Acta2, Cnn1, and Actg2,^34,35 and mesenchymal cells (∼2.06%) were identified based on the classic markers of mesenchymal cells, including Dcn, Lum, and Gsn.^36,37 Immune cells were determined based on the classic markers Ptprc (CD45), Cd3e, and Cd83.³⁸ A cluster termed endothelial cells was determined according to classic markers, such as Cdn5, vWf, and Pecam1(Fig. 1D).³⁹ Then, we performed immunofluorescent staining using specific markers for acinar cells, ductal cells, and myoepithelial cells (Fig. 1E). We also identified a cluster of highly expressed Kit⁺ cells (∼5.73%) and a cluster of mitotic cells with unknown functions (Figs. 1B, 1C). Collectively, we demonstrated the cellular diversity of the LG. 

Subsequently, the biological functions of differentially expressed genes (DEGs) in each cell type were investigated using Gene Ontology (GO) (Fig. 1E). Specifically, the significantly enriched GO terms in acinar cells mainly include Golgi vesicle transport, protein exit from the endoplasmic reticulum, and Golgi vesicle budding. The enriched GO terms associated with ductal cells include urogenital system development, epithelial tube morphogenesis, and actin filament organization. The major functions of MECs include axonogenesis, muscular system processes, and muscle tissue development. Mesenchymal cells are functionally associated with extracellular organization, cell–substrate adhesion, and connective tissue development. The enriched pathways in immunocytes include T-cell activation, immune response–regulating signaling pathways, and positive regulation of cytokine production. The endothelial cells mainly regulate vasculature development, circulatory system processes, and cell junction organization. Moreover, when compared with the Ctrl group, a notable reduction in the proportion of acinar cells in the SCOP LGs was observed (Fig. 1C). Conversely, the ratios of ductal cells, MECs, mesenchymal cells, and immunocytes in SCOP-treated LG increased (Fig. 1C). These results suggest remarkable alterations in the proportion of cell types across the Ctrl and SCOP LG samples. 

Transcriptomic Alterations of LG Acinar Cells in Control and SCOP Mice

Unsupervised subclustering was performed to understand the heterogeneity of LG acinar cells. We found that the acinar cells could be divided into two subclusters: 50.17% C0 and 49.83% C1 (Fig. 2A). As shown in Supplementary Fig. S2A, these two clusters were uniformly distributed in the control and SCOP samples in different proportions. GO analysis of each subcluster was conducted based on the highly expressed genes (Supplementary Fig. S2B). The C0 subcluster is functionally associated with mRNA processing, covalent chromatin modification, GTPase activity regulation, and Golgi vesicle transport. Cytoplasmic translation and ribosomal small subunit biogenesis were enriched in the C1 subcluster (Supplementary Fig. S2B). The top 10 marker genes for each subtype are shown in Figure 2B. The C0 group contained Pded4d, Auts2, Naalad12, Zbtb20, Tmem117, Chm3, Thrb, Gm20275, and Esrrg, and the C1 group contained Bpifa2, Smr3a, Amy1, Scgb1b7, Obp1b, Scgb1b20, Mucl2, Scgb2b7, Scgb2b20, and Gm14743. The dot plot also showed high expression of LPO (lactoperoxidase) and Egf in the C0 group and elevated expression of Mucl2 (mucin-like2), Car6 (carbonic anhydrase VI), and Ltf (lactoferrin) in the C1 group (Fig. 2C). Most of these genes have been reported to be associated with acinar cell function.^12,32 The findings based on GO analysis and differential gene expression patterns revealed secretory functional heterogeneity of these two acinar cell subclusters. 

Figure 2.

 

The transcriptional heterogeneity of LG acinar cells from different conditions. (A) The t-SNE clustering of acinar cells. (B) Top 10 markers of acinar cells. (C) Dot plot showing the expression of representative protein for each cell type. (D) Differential genes between two groups for SCOP versus Ctrl conditions. (E) Representative shared GO terms for SCOP upregulation of two distinct subclusters in different acinar cell types. The color keys from yellow to red indicate the range of P values. (F–H) Violin plots showing expression levels of secretion, water transport osmotic stress, and two distinct subclusters in the Ctrl and SCOP groups. Two-tailed t-test P values are indicated.

Subsequently, we investigated cellular and transcriptional alterations in different acinar cell subpopulations under DED conditions. As shown in Supplementary Fig. S2C, the proportion of the C1 subcluster sharply decreased after treatment with SCOP, and the C0 group showed a slight reduction. More importantly, when compared with the Ctrl, 459 upregulated and 177 downregulated genes were identified in the C0 subcluster of the SCOP-treated group, as well as 1411 upregulated genes and 153 downregulated genes in the C1 subcluster (Fig. 2D). We observed decreased expression of LPO, Mucl2, and Car6 in both the C0 and C1 subclusters following SCOP induction (Supplementary Fig. S2D). Through GO analysis of the genes upregulated in the SCOP group, we observed significant enrichment in the positive regulation of the mitotic cell cycle in the C0 subcluster and multiple enriched GO terms in the C1 subcluster, including positive regulation of mitotic cell cycle, positive regulation of catabolic process, negative regulation of translation, negative regulation of protein phosphorylation, and negative regulation of cell development (Fig. 2E). In contrast, downregulated genes in the C0 and C1 subclusters were enriched both in cytoplasmic translation and ribosome biogenesis (Supplementary Fig. S2E). These results demonstrate reduced tear secretion function. Moreover, when compared with the Ctrl group, the score of negative regulation of secretion in both the C0 and C1 subclusters was dramatically increased after SCOP treatment (Fig. 2F), as well as a lowered score of water transport (Fig. 2G) and a higher score of response to osmotic stress (Fig. 2H), suggesting dysfunction of acinar cells in the DED scenarios. These results demonstrated that the LG contained two subpopulations of acinar cells whose pathological alterations contributed to the pathogenesis of DED. 

Transcriptomic Heterogeneity of LG Ductal Cells From Control and SCOP-Induced DED Mice

Similar to acinar cells, ductal cells were classified into two separate subclusters (Fig. 3A), with 61.57% C0 subclusters and 38.43% C1 subclusters. These two subclusters were uniformly distributed in different proportions in the Ctrl and SCOP samples (Supplementary Fig. S3A). GO analysis revealed that the significantly expressed genes in subcluster C0 were functionally associated with the regulation of ion transmembrane transport and response to osmotic stress, whereas subcluster C1 was mainly enriched in the circulatory system process, urogenital system development, and epithelial tube morphogenesis (Fig. 3B). These results demonstrate the heterogeneity of LG ductal cells, including two subpopulations with different functions. 

Figure 3.

 

The transcriptomic profiling of LG ductal cells under different conditions. (A) The t-SNE clustering of ductal cells. (B) Expression patterns of the top 50 differentially expressed genes for each subtype (left) and GO terms of all differentially expressed genes in each subtype (right). (C) Differential genes between two groups under SCOP versus Ctrl conditions. (D) Representative shared GO terms of upregulated genes of two distinct subclusters after SCOP treatment. The color keys from yellow to red indicate the range of P values. (E) Heatmap of abnormal lacrimal ductal morphology–related genes; the genes highlighted in red boxes were selected for validation. (F) Immunofluorescence revealed that Ndufb11 is expressed only in the acinar cells of Ctrl LGs, whereas Fgfr2 and Sox6 are expressed exclusively in the LGs of SCOP mice.

Following SCOP treatment, there was a noticeable increase in the proportion of both C0 and C1 subclusters of ductal cells (Supplementary Fig. S3B). Under SCOP conditions, 1409 upregulated genes and 151 downregulated genes were observed in C0 ductal cells, and 1125 upregulated genes and 159 downregulated genes were observed in C1 ductal cells (Fig. 3C). Furthermore, significantly enriched GO terms of the genes upregulated in C0 ductal cells were identified, including histone modification, cell–cell signaling by Wnt pathways, positive regulation of the cell cycle, regulation of chromosome organization, regulation of microtubule-based processes, and negative regulation of phosphorylation (Fig. 3D). The same terms were also observed for the upregulated genes in C1 ductal cells (Fig. 3D). Nevertheless, the downregulated genes in both the C0 and C1 subclusters were mainly enriched in cytoplasmic translation, ribosome biogenesis, ribosome assembly, and rRNA processing (Supplementary Fig. S3C), demonstrating a lower protein synthesis ability. Moreover, compared to the corresponding Ctrl subcluster, several abnormal lacrimal ductal morphology–related genes were observed in the SCOP-treated C0 and C1 subclusters, including Ndufb11, Sox6, Col3a1, Tfap2a, and Tet3 (Figs. 3E, 3F). These findings systematically revealed prominent alterations in the transcriptional profiles of different ductal cell subclusters under DED scenarios. 

Transcriptomic Characterization of LG MECs From Normal and SCOP-Induced DED Mice

As one of the major LG cell types, MECs surround the clusters of acinar cells and form acini together with acinar cells, playing an important role in tear secretion.⁴⁰ To investigate their heterogeneity under physiological conditions, MECs were divided into three distinct subclusters (Fig. 4A, Supplementary Fig. S4A). The top 10 DEGs in each subcluster are shown in Figure 4B. Several representative MEC-related genes were highly expressed in all three subclusters, including Acta2, Myh11, and Mylk (Fig. 4C).^41,42 

Figure 4.

 

The snRNA-seq analysis of LG MECs under different conditions. (A) The t-SNE clustering of MECs. (B) Top 30 markers of MECs. (C) Dot plot showing the expression of representative proteins for each cell type. (D) Differential genes between two groups under SCOP versus Ctrl conditions. (E) Functional enrichment map of DEGs shared by three subclusters of MECs. (F) Gene scoring analysis using muscle contraction related genes. ****P < 0.0001 (two-sided Wilcoxon rank-sum test).

Next, we investigated characteristic alterations in different MEC clusters under SCOP-induced DED. We found that the percentage of each MEC subcluster in the SCOP group was elevated compared to the corresponding subcluster in the Ctrl group, suggesting a significant change in the composition ratio of MEC under DED conditions (Supplementary Fig. S4B). A total of 1892 upregulated and 213 downregulated genes were identified in at least one MEC subcluster (Fig. 4D). Among these DEGs, 27.8% were shared in at least two subclusters, and the majority of DEGs were cell type specific (Fig. 4D). GO analysis revealed that the upregulated genes were closely associated with mRNA processing, muscle system processing, and muscle tissue development. However, the downregulated genes were mainly enriched in cytoplasmic translation, ribosome biogenesis, and the regulation of cellular amide metabolic processes. Moreover, the score of the DEGs responsible for muscle contraction in each MEC subcluster from the SCOP group was higher than that of the corresponding subcluster from the Ctrl group (Fig. 4F). These findings demonstrate a remarkable change in the transcriptional profile of MEC in SCOP-induced DED. 

Heterogeneity of LG Immunocytes Clusters From Normal and SCOP-Induced DED Mice

Subsequently, we explored the cellular diversity and heterogeneity of immune cells in Ctrl and SCOP-treated LGs. The immunocytes in the control LGs were categorized into four subclusters based on representative cell markers: mono/mac cells, natural killer (NK)/T cells, plasma cells, and dendritic cells (DCs) (Fig. 5A). The specific markers corresponding to the different clusters are shown in Figure 5B, with mono/mac cells identified by Csf1r, Itgam, and Adgre1 (F4/80)^43,44; NK/T cells by Klrb1, Cd3d, and Cd3e^45–47; plasma cells by Jchain, Igkc, and Mzb1^48,49; and DCs by Ccr7, Cacnb3, and Fscn1.^50–52 Moreover, the proportion of mono/mac cells in LC immunocytes was as high as 71.4%; NK/T cells accounted for 22.8%, and plasma cells and DCs accounted for only 4.58% and 1.34%, respectively (Fig. 5C). 

Figure 5.

 

Transcriptional characterization of different immunocyte populations under different conditions. (A) The t-SNE clustering of immune cells. (B) Characteristic markers of various immune cells. (C) Proportion of different cell types in normal LGs (fan chart) and proportion of different cell types in the SCOP and Ctrl groups (histogram). (D) The number of genes expressed by different cells under SCOP and Ctrl conditions. (E) Violin plots showing the expression levels of abnormal inflammatory response in different immunocyte subclusters.

However, under SCOP-induced DED conditions, the percentage of most immunocyte subclusters in the LG dramatically increased, including mono/mac, NK/T cells, and DCs (Fig. 5C). Through single-cell transcriptomic analysis, the DEGs of each immunocyte subpopulation in the SCOP LGs were identified, including 1879 upregulated and 170 downregulated genes in mono/mac cells, 1525 upregulated and 143 downregulated genes in NK/T cells, 409 upregulated and 81 downregulated genes in plasma cells, and 603 upregulated and 58 downregulated genes in DCs (Fig. 5D). Inflammatory signature analysis revealed that the score of each immunocyte subpopulation in the SCOP LGs was much higher than that in the Ctrl LGs, suggesting an increased abnormal inflammatory response in SCOP-induced DED (Fig. 5E). These findings demonstrate that the altered percentage and DEG profile of each immunocyte subcluster probably resulted in an inflammatory response caused by SCOP/DS. 

Cell–Cell Communication Networks Among Multiple Immunocyte-Types With Different Parenchymal Cells in Mouse LGs with Different Conditions

To systematically understand the pathogenesis of SCOP/DS-induced DED, we dissected the ligand–receptor interactions between different cell subclusters based on cell–cell interaction analysis, especially the associations between different immunocyte subpopulations and LC parenchymal cells. As shown in the chord diagram, a large number of interactions between immunocyte subpopulations and LG parenchymal cells were observed in the Ctrl and SCOP samples, including 520 interactions in the Ctrl and 1052 in the SCOP group (Fig. 6A). Compared to the Ctrl group, the number of connections between immunocyte subpopulations and LG parenchymal cells (mainly acinar cells, ductal cells, and MECs) was obviously increased in the SCOP group (Fig. 6B). 

Figure 6.

 

The interactions of different immunocyte-clusters with parenchymal cells in LGs based on ligand–receptor analysis. (A) Chord diagram of cellular interactions between immune cells and acinar cells, ductal cells, and MECs from the Ctrl and SCOP groups, separately. The number of interaction pairs and cell types are annotated. (B) The number of cell interactions between immune cells and acinar cells, ductal cell, and MECs from the Ctrl and SCOP groups. (C) Network plots showing the ligand–receptor interaction events between immune cells and acinar cells/ductal cells/MECs. Cell–cell communication is indicated by the connected line. The thickness of the lines and the size of the nodes are positively correlated with the number of ligand–receptor interaction events. (D) Dot plot of erased ligand–receptor interactions associated with acinar cells, ductal cells, and MECs in the SCOP group compared to those in the Ctrl group. Acinar cells, ductal cells, and MECs express receptors and receive ligand signals from immune cells. The row represents a ligand–receptor pair, and the column defines a cell–cell interaction. The samples of the Ctrl and SCOP groups are indicated by brackets. P values and means were calculated using the CellPhoneDB pipeline. The dot size reflects the P values for cell type specificity, and the dot color denotes the mean of the average ligand–receptor expression in the corresponding interacting cell types.

Moreover, network plots revealed remarkable alterations in the communication between immunocyte subpopulations and LG parenchymal cells under DED conditions, along with the gain of numerous interaction pairs and the loss of other interaction pairs (Fig. 6C). Based on the crucial connections of mono/mac cells with LG acinar cells, ductal cells, and MECs (Fig. 6C), we evaluated the complicated communication between mono/mac cells and different LG parenchymal cells under DED scenarios (Fig. 6D). Notably, the connection between insulin-like growth factor 1 and its receptor was strengthened in all cell types in the SCOP-treated group. This correlation has been verified to be linked with DNA damage repair, cellular aging, and antioxidative stress effects, which probably contribute to the development of dry eye.^53,54 The interaction of transforming growth factor β1(TGF-β1) with different receptors was also observed across a wide range of cell types, which was implicated in tissue fibrosis and epithelial–mesenchymal transition. Remarkably, these biological processes have been reported in the LGs of mice with dry eye induced by SCOP.^55,56 Altogether, disordered cell–cell communication between immunocytes and LG parenchymal cells provides novel insights into the pathogenesis of DED. 

Discussion

There was no substantial disparity in the ability of scRNA-seq and snRNA-seq to distinguish among different cell types. However, for substantial tissues such as the brain, heart, and kidneys, scRNA-seq often leads to the separation of easily detachable cell types, resulting in the loss of less easily detached cells.^57,58 Additionally, some fragile cells may fracture due to excessive separation; consequently, the dissociation process is unable to adequately capture all cell types in the tissue, resulting in a notable decline in the accuracy of identifying cell proportions.⁵⁹ Currently, there is no research that analyzes LG cells at the single-nucleus level, and our study addresses the deficiency in LG cell categorization. In this study, we initially constructed a very accurate single-cell atlas that encompassed both the prevalent and rare cell types in the LG. Simultaneously, we clearly demonstrated the relative proportions of different cells in the mouse LG, a feature that cannot be achieved using conventional methods of scRNA-seq. 

As an essential organ controlling tear production, the LG plays a critical role in maintaining homeostasis of the ocular surface.^6,60 Due to its functional relevance in numerous ocular diseases, its cellular and functional heterogeneity under physiological and pathological conditions has attracted growing attention in the field of ophthalmology. In this study, we offer a comprehensive single-nucleus atlas of the LG under normal and DED scenarios. Normal LG showed cellular diversity, including acinar cells, ductal cells, MECs, mesenchymal cells, immune cells, endothelial cells, Kit⁺ cells, and a cluster of mitotic cells. However, LGs under DED conditions were characterized by a reduced proportion of acinar cells with reduced tear secretion ability, a higher percentage of ductal cells and MECs, and more immune infiltration with an increased inflammatory response. Moreover, dysregulated communication between different immunocytes and LG parenchymal cells was observed in SCOP-induced DED. This study provides a cellular outline for understanding the biology of LG and the pathogenesis of DED. 

Although several studies concerning scRNA-seq of the LG have been conducted,^10,12,57 the proportion of different cell types in the LG remains uncertain. Through snRNA-seq, we showed that acinar cells accounted for 72.4% of the LG under physiological conditions, suggesting that acinar cells are important cells in the LG for tear secretion.⁶¹ The C0 cell subtype predominantly releases LPO and the cytokine epidermal growth factor (EGF), both of which are linked to the antibacterial and defensive roles of the ocular surface.^62,63 Meanwhile, the Car6 protein produced by the C1 subtype is involved in maintaining the pH homeostasis of tears.^12,32 These findings indicate that these two kinds of acinar cells probably work together to regulate the makeup of the tear fluid. However, under SCOP/DS-induced DED conditions, several enriched pathways of the upregulated DEGs in both C0 and C1 cells were observed, including the mitotic cell cycle, inhibition of cell development, and excessive catabolic processes. More interestingly, pharmacological blockade of aberrant cell proliferation rescued the tear reduction in SCOP-induced DED.⁶¹ These abnormal pathways therefore are likely involved in the loss of acinar cells and decreased secretion function.^62–64 

Previous studies revealed that ductal cells were mainly involved in the release of water and ions into the tear fluid.⁶⁵ We found that the G0 subcluster was mainly responsible for regulating the ion transmembrane and responding to osmotic stress and may be the major cells responsible for the release of water and ions into the tear fluid. G1 was more inclined toward duct morphology and development, as supported by urogenital system development and epithelial tube morphogenesis. These results indicate the functional heterogeneity of LG ductal cells. However, the transcriptional profile and enriched pathways were greatly altered in both G0 and G1 phases in the SCOP-induced DED model, characterized by the enrichment of histone modifications, positive regulation of the cell cycle, cytoplasmic translation, and ribosome biogenesis, as well as the abnormal expression of several genes governing lacrimal ductal morphology. These results provide strong evidence for further investigation into the role of ductal cells in the pathogenesis of DED. 

As an important component of the secretory apparatus of the LG, MECs are capable of contracting to expel lacrimal fluid.⁵⁵ Due to the expression of muscarinic and purinergic receptors, MECs can be regulated by neural innervation.⁶⁶ Via snRNA-seq analysis, MECs were classified into three subclusters, with the C0 subcluster in a high proportion and greater expression of classic MEC marker genes (such as Acta2, Myh11, and Myk^41,67,68). In contrast, in the SCOP-induced DED model, the percentage of the three subclusters of MECs in the LG increased, especially for the C0 subcluster. SCOP has been reported to be a competitive antagonist of acetylcholine that lowers the contractile capacity of smooth muscles.^14,15 More surprisingly, functional analysis of the upregulated genes revealed a higher potential for muscle contraction and muscle tissue development, which suggests the existence of a negative-feedback mechanism in SCOP-treated MECs. The exact roles of different MECs in the pathogenesis of SCOP-induced DED must be explored in detail. 

Several studies have systemically revealed the existence of various immunocyte types in LG under both normal and DED conditions.^12,61,69 Through snRNA-seq analysis, the cellular composition of LGs under normal conditions was identified, including mono/mac cells, NK/T cells, plasma cells, and DCs. Their interactions with the LG parenchymal cells were also determined. However, in SCOP-induced DED, a remarkable increase in the proportion of various immunocytes in the LG was observed, as well as increased interactions with LG parenchymal cells, especially mono/mac cells. Existing evidence has demonstrated the involvement of immune cells in the development of DED,^70–72 which supports the important role of immunocytes in DED development. However, the exact roles of different immunocyte populations in the maintenance of LG homeostasis and pathogenesis of DED must be examined. 

We also evaluated intercellular interactions in the LG based on receptor–ligand interactions in the Ctrl and DED groups. Multiple interaction pairs have been identified between immunocytes and LG parenchymal cells, suggesting that different immunocytes play critical roles in the regulation of tear secretion and DED development. The number of receptor–ligand pairs after SCOP treatment was significantly increased but was accompanied by the loss or intensification of some interaction pairs. After SCOP treatment, a decreased population of acinar cells and an increased proportion of ductal cells and MECs were observed. These phenomena probably resulted from strengthened TGF-β signaling, which has been reported to be closely associated with acinar-to-ductal metaplasia in pancreatic related studies.^70–73 Nevertheless, further studies are imperative to address these complex interactions in maintaining LG homeostasis and regulating DED development. 

This study has several limitations. Although relevant characteristics and states regarding LGs in both control and DED have been observed using snRNA-seq, the findings and hypotheses must be validated at the gene or protein level. Moreover, the function of different cell types in regulating tear secretion under both control and DED conditions should be investigated in future studies, as well as the underlying mechanisms. Finally, the transcriptional and functional heterogeneity of other cell types identified in this study, such as endothelial cells, Kit⁺ cells, and a cluster of mitotic cells, should be investigated and consolidated in future studies. 

In summary, our study offers a systemic single-nucleotide transcriptomic profile for normal LGs, as well as cellular and functional heterogeneity of LGs between the Ctrl and DED groups. These findings provide new insights into LG homeostasis and the pathogenesis of DED, and our study may provide some potential guidance for future treatments of dry eye, especially aqueous deficiency dry eye, but further mechanistic research is necessary. 

Acknowledgments

Supported by grants from the Shandong Provincial Key Research and Development Program (2021ZDSYS14) and Academic Promotion Programme of Shandong First Medical University (2019ZL001). 

Disclosure: Y. Tang, None; S. Dou, None; C. Wei, None; Z. Sun, None; D. Sun, None; Q. Zhou, None; L. Xie, None