Mechanisms of Extraorbital Lacrimal Gland Aging in Mice: An Integrative Analysis of the Temporal Transcriptome

Abstract

Purpose: This study used high-throughput RNA sequencing (RNA-Seq) and bioinformatics analysis to investigate the altered transcriptome profile of aging lacrimal glands in mice that occurs over the course of a 24-hour cycle.

Methods: Male C57BL/6J mice aged 12 weeks (young) and 20 months (aging) were housed in a pathogen-free setting with a 12-hour light/12-hour dark cycle. Throughout a 24-hour cycle, mouse extraorbital lacrimal glands (ELGs) were collected at eight time points at three-hour intervals. To prepare for the high-throughput RNA-Seq, whole mRNA was extracted. Differentially expressed genes (DEGs) in the young and aging groups were subjected to bioinformatic analysis based on diurnal patterns. Furthermore, the cell populations in which significant DEGs express and signaling pathways occur were validated at the single-cell RNA sequencing (scRNA-seq) level.

Results: The total transcriptome composition was significantly altered in aging ELGs compared with that in young mouse ELGs at eight time points during the 24-hour cycle, with 864 upregulated and 228 downregulated DEGs, which were primarily enriched in inflammatory pathways. Further comparative analysis of the point-to-point transcriptome revealed that aging ELGs underwent alterations in the temporal transcriptome profile in several pathways, including the inflammation-related, metabolism-related, mitochondrial bioenergetic function–associated, synaptome neural activity–associated, cell processes–associated, DNA processing–associated and fibrosis–associated pathways. Most of these pathways occurred separately in distinct cell populations.

Conclusions: Transcriptome profiles of aging lacrimal glands undergo considerable diurnal time-dependent changes; this finding offers a comprehensive source of information to better understand the pathophysiology of lacrimal gland aging and its underlying mechanisms.

Aging is a universal biological process that results in a gradual decline in the body's tissue structure and physiological functions.¹ Aging is a major risk factor for death from all age-related chronic diseases in adults. Globally, the number of people aged over 65 years is expected to increase from 617 million to more than 2 billion by 2050, representing 20% of the world's population.¹ This demographic shift poses various challenges from both a healthcare and a socioeconomic perspective. Current evidence suggests that manipulation of the biological processes that occur with aging can ameliorate or delay many age-related diseases.² Therefore it is important to explore the mechanisms underlying the onset of aging and its associated diseases. 

As the body ages, the intraocular tissues and their appendages, important components of the visual system, undergo age-related pathologic changes.^3,4 The elasticity and transparency of the lens decrease with age, resulting in presbyopia and varying degrees of opacity. In the retina, the incidence of age-related macular degeneration increases with age.^5–7 Similarly, the incidence of glaucoma tends to increase significantly with age.^8–11 The meibomian glands, which produce the lipid layer of the tear film, undergo significant age-related changes at the in vivo structure,¹² histology,¹³ and transcriptome levels.¹⁴ Additionally, the integrity of the ocular surface is significantly altered pathologically with aging, including decreased tear film stability, reduced goblet cell numbers,¹⁵ impaired corneal barrier function,¹⁶ and altered immune cell distribution and response.¹⁷ 

The lacrimal glands also undergo varying degrees of degeneration with age, and when the secretion of the aqueous layer from the lacrimal gland decreases with age, age-related dry eye syndrome is manifested.¹⁸ A meta-analysis reported that the prevalence of dry eye disease in the elderly population aged 60 years or older increases from 1.3% to 9.2% compared to the younger population.¹⁸ The high prevalence of age-related dry eye disease is a significant concern that profoundly affects the vision quality of the elderly population.¹⁹ Therefore comprehending the process of lacrimal gland aging and its underlying mechanisms is crucial for enhancing the quality of life for the elderly population. 

Clinical and experimental evidence suggests that the major pathological changes associated with lacrimal gland aging include chronic inflammation, diffuse atrophy, and periductal fibrosis.^20–22 The underlying mechanisms of these pathological changes may be primarily related to the following: (1) chronic inflammation triggered by an imbalance between anti-inflammatory and pro-inflammatory mechanisms during aging, as chronic inflammation is a significant contributor to the aging process^23,24; (2) decreased nerve density and neurotransmission dysfunction²⁵ leading to the impairment of lacrimal gland's secretory function; (3) damage to DNA transcription and repair due to oxidative stress^26–28; and (4) dysregulated communication between the mammalian genome and the gut microbiome that occurs with age²⁹ and leads to damage to peripheral organs or tissues, including the lacrimal gland.³⁰ The above data and information provide important clues for our understanding of the aging mechanisms of the lacrimal gland. However, these findings are based on general pathological and biological techniques using model animals and human samples at single time points, and our understanding of complex mechanisms remains remarkably limited. Therefore it is necessary to use new technical tools to revalidate these findings. 

Bulk RNA sequencing (RNA-seq) data from humans and rodents in recent years have confirmed age-dependent changes in the transcriptome across multiple organs and different ages, highlighting in particular the relevance of age-related chronic inflammation enhancement and loss of genomic stability.^31,32 The concomitant age-related decline in circadian physiological behaviors and alterations in circadian transcriptome profiles have been widely recognized.³³ Results from several aging organs or tissues show that the core clock gene machinery remains largely stable in aged tissues, but the expression of core clock–driven output genes is dramatically altered.^34–38 Similarly, our previous findings in the aging mouse lacrimal gland show stability in core clock genes and dramatic alterations in circadian output genes.³⁰ However, there is extremely limited understanding of the diurnal time-dependent changes in the overall output transcriptome of the aging lacrimal gland. 

To better understand this issue, we analyzed the changes in the basic composition of the transcriptome of aging mouse extraorbital lacrimal glands (ELGs) by sampling the ELGs at eight time points at three-hour intervals in a 24-hour cycle. Subsequently, we systematically analyzed and compared the transcriptome data from young and aging mouse ELGs in a point-by-point, multifaceted manner based on these circadian time patterns. Our data revealed significant time-dependent changes in the output transcriptome of aging mouse lacrimal glands compared to those of young mice. In addition, we specifically validated some aspects of the above transcriptome-related information based on single-cell RNA sequencing (scRNA-seq) technology. Based on these data, we discuss the possible mechanisms leading to lacrimal gland senescence, which may provide a basis for further exploration of the pathogenesis of lacrimal gland aging and the development of novel time-based strategies to delay the lacrimal gland aging process. 

Material and Methods

Experimental Design

ELGs were collected from young (12-week-old) and aged (20-month-old) mice at eight time points at three-hour intervals in a 24-hour cycle. Bilateral ELGs from each mouse were pooled into one sample, and three biological replicates were collected at each time point. After extraction of total mRNA, high-throughput sequencing was performed (Fig. 1A). DESeq2 was used to screen for differentially expressed genes (DEGs) between the two sets of ELGs. The following bioinformatic analyses were performed: (1) the differences in transcriptome composition between young and senescent ELGs were investigated by dimensionality reduction of the data using principal component analysis (PCA); (2) the genes with the largest differences between the selected samples were visualized using heatmaps, where the magnitude (x-axis) and statistical significance (y-axis) of the differences were shown using volcano plots; (3) Kyoto Gene and Genome Encyclopedia (KEGG) annotation was used to understand the major pathways involved in these DEGs; (4) gene ontology (GO) analysis was used to understand the biological processes (BPs) involved in these DEGs; and (5) gene set enrichment analysis (GSEA) was used to determine whether predefined gene sets showed a statistically significant increase or decrease in young and old ELGs (Fig. 1B). Finally, some of the data were validated by immunohistochemistry, hematoxylin and eosin (H&E) staining, tear secretion measurement, and scRNA-seq (Fig. 1C). 

Figure 1.

 

Flowchart of experimental design and data analysis. (A) Sample preparation and high-throughput RNA-seq for ELGs. (B) Informatics analysis protocol. (C) Verification and analysis.

Experimental Animals

Specific pathogen-free and ocular disease-free male C57BL/6 mice were purchased from GemPharmatech LLC (Nanjing, China). To exclude the influence of environmental and dietary factors, all animals were housed in the same environment and received the same diet. All animals were housed in circadian chambers (Longer-Biotech Co., Ltd, Guangzhou, China) to maintain a 24-hour circadian state (12-hour light/12-hour dark). Mice aged 12 weeks were used as young controls, whereas mice aged ≥20 months were used as aging experimental groups.³⁹ Time was expressed using the zeitgeber time (ZT) scale as an indicator of rhythm phase, where ZT0 and ZT12 refer to light-on time (7 AM) and light-off time (7 PM), respectively.^40,41 Mice had as desired access to standard chow and water. All animal experiments were performed in accordance with the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Vision and Ophthalmic Research and were approved by the Institutional Animal Care and Use Committee of Henan Provincial People's Hospital (No. HNEECA-2023-02). Animals were euthanized by inhalation of isoflurane and cervical dislocation at the end of the experiment. 

Measurement of Lacrimation

Tear secretion was measured by intraperitoneal injection of pilocarpine hydrochloride (4.5 mg/kg) into mice (n = 6 at each time point) at ZT0, 6, 12, and 18, as previously described.^42,43 Ten minutes after drug administration, phenol red threads (no. 30059010; Tianjin Jingming New Technology Development Co., Tianjin, China) were placed in the inner canthi of the mice for 20 seconds. The length of the red part (mm) was measured to determine the amount of tear secretion. 

Behavioral Activity Monitoring and Analysis

Behavioral activity was monitored, recorded, and analyzed according to previously established methods.⁴⁴ Specifically, mouse locomotor activity and core body temperature were monitored using a telemetry system (model ER-4000 E-Mitter; Mini Mitter, Sunriver, OR, USA). The E-Mitter was initially implanted in the peritoneal cavity of anesthetized mice (n = 3 for each group). After two weeks of recovery, the radiotelemetry receiver was used to receive signals. Locomotor activity and core body temperature were recorded on a configured computer system. The former was monitored every five minutes, and the latter was monitored every 20 minutes. 

ELG Collection and Total RNA Extraction

ELGs were collected and their total RNA was extracted per our previously described protocol.^30,42,43 ELGs from euthanized animals were rapidly collected at ZT0, 3, 6, 9, 12, 15, 18, and 21 and quickly placed in liquid nitrogen for cryopreservation. Total RNA in the ELGs was extracted using the RNAeasy Spin Column Kit (Qiagen, Hilden, Germany) according to the instructions. To ensure that samples collected at dark-cycle time points (ZT12, 15, 18, and 21) were not affected by light, tissue collection was performed under a red light source (three-watt safety lamp, DG-20A; Jining Hengshuo Testing Instruments Co., Ltd.). 

Immunostaining of Murine ELGs

Histological observation of the lacrimal glands was performed as described previously.^30,42,43 Four lacrimal glands from young and aging mice were removed after euthanasia and immediately fixed in formalin solution for 24 hours. Lacrimal glands were sectioned sagittally along the longitudinal axis of the lacrimal gland, deparaffinized, and rehydrated in xylene and 75% to 100% ethanol. Sections were stained by overnight incubation at 4°C with anti-mouse CD3 antibody (no. GB16669; Servicebio Company, Wuhan, China) for T lymphocytes, with anti-mouse CD19 antibody for B lymphocytes (no. GB11061-1; Servicebio Company), with anti-NKp46 antibody for natural killer (NK) cells (no. MAB22252, clone 29A1.4; R&D Systems, Minneapolis, MN, USA), and with fluorescein isothiocyanate–conjugated anti-mouse beta III tubulin monoclonal antibody (no. GB12139; Service-bio) for the nerve fibers in ELGs. ImageJ software (version 1.42) was used to calculate the percentage of positive cells and green staining (n = 6). In addition, some sections were stained with H&E to evaluate morphology. 

Bulk RNA-seq

The bulk RNA-seq was performed as described previously.^30,42,43 The poly-A-containing mRNA from the total RNA isolated in the ELGs was first purified using poly-T oligo-attached magnetic beads, and the purified mRNA was then decomposed into small fragments using divalent cations at high temperature. The first strand cDNA was then synthesized using reverse transcriptase and random primers, and the second strand cDNA was ligated to the adapter and synthesized by adding polymerase I, ribonuclease H, a buffer, and dNTPs. These cDNA fragments were then enriched by PCR amplification, and the PCR products were quantified using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Inc, Waltham, MA, USA) and pooled into single-stranded DNA loops to form the final library. Sequencing was performed on the BGISEQ-500 platform (BGI, Guangdong, China) after checking the fragment size and concentration of the final library on the Agilent 2100 Bioanalyzer (Agilent DNA 1000 reagents; Agilent Technologies, Santa Clara, CA, USA). The sequencing strategy was single-end 50 bp with a sequencing depth of 20 million. The sequencing quality scores Q20 and Q30 were both greater than 90%. The data were filtered using the SOAPnuke software (version 1.5.2, https://github.com/BGI-flexlab/SOAPnuke) to obtain clean reads, which were mapped to the reference gene set and genome using Bowtie2 software (version 2.2.5, https://sourceforge.net/projects/bowtiebio) and HISAT software (version 2.0.4, http://www.ccb.jhu.edu/software/hisat/index.shtml), respectively. Last, the gene expression levels were calculated using RSEM software (v1.2.12, https://github.com/deweylab/RSEM). It should be noted that the analysis of the core clock genes and core clock-controlled output circadian genes regarding circadian rhythms in the data set of this study was reported previously.³⁰ The present study focused solely on the analysis of the diurnal time-dependent changes in the total output transcriptome of aging lacrimal glands (all raw data are available from NCBI's BioProject database under accession PRJNA1000710). 

ScRNA-Seq

Preparation of Single Cell Suspensions

All ELGs were harvested at the same time (ZT3–4) for each animal to avoid the influence of the circadian rhythm.^45,46 Freshly obtained ELG samples from young (12 weeks old, n = 3) and aged mice (≥20 months old, n = 3) were first chopped into approximately 1 mm² pieces of tissue. The tissue was then enzymatically digested into a single cell suspension using collagenase D (0.5 mg/mL) and DNase I (0.2 mg/mL) and incubated in a shaker at 37°C for 30 to 45 minutes. The digested tissue was then passed through a 70 µm filter and treated with an erythrocyte lysis buffer. 

Preprocessing and Quality Control of ScRNA-Seq Data

Single-cell libraries were generated using the 10× Genomics Gel Bead Kit. Libraries were sequenced on the Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA). To generate normalized summary data across samples, Cell Ranger software (version 5.0.0, https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger) was used to process raw data, perform multiple decomposition of cell barcodes, and map and down-sample transcriptome reads. This resulted in a raw unique molecular identifier (UMI) count matrix, which was converted to Seurat objects using the R package Seurat software (version 3.1.1). 

To ensure data quality, low quality cells and likely multiplet captures were filtered out. Cells with UMI/gene counts exceeding a mean value of ±2 times of the standard deviation were excluded, assuming a Gaussian distribution of UMI counts per gene. Cells with more than 10% of counts belonging to mitochondrial genes were also removed. In addition, the DoubletFinder package (version 2.0.2) was used with default parameters to identify potential doublets.⁴⁷ Interference signals from populations present in only one sample were also eliminated. The count matrix of each sample was merged using the cbind function in R, followed by log-normalization using the NormalizeData function in Seurat. To generate an integrated data matrix, the ScaleData function in Seurat was employed for scaling. To remove batch effects, mutual nearest neighbors analysis was conducted for dimensionality reduction using the batchelor package. Marker genes for each cluster were identified using the FindAllMarkers function in Seurat. Finally, a total of 74,728 cells were obtained in this experiment, with an average of 31,821 reads per cell and an average of 1109 genes per cell. 

Visualization, Clustering, and Cell Annotation

The count matrices for each sample were merged using the cbind function in R and logarithmically normalized using the NormalizeData function in the Seurat package. Finally, scaling was performed using the ScaleData function in the Seurat package to create a combined data matrix. To eliminate batch effects in the scRNA-seq data, mutual nearest neighbor analysis was performed on top of the dimensionality reduction of the scaled comprehensive data matrix using the batchelor package. The FindAllMarkers function in Seurat was used to identify marker genes for each cluster. 

Analysis of Transcriptome

DEGs

The screening of the list of key genes on pathways related to innate immunity and inflammation, metabolism, bioenergetics, neural activity, and growth signaling was mainly obtained from the KEGG database (https://www.genome.jp/kegg/) and GeneCards database (https://www.genecards.org/). DEGs between the young and aging ELGs were screened using the DESeq2 (differential gene expression analysis based on the negative binomial distribution) package based on R (version 1.32.0; https://bioconductor.org/packages/release/bioc/html/DESeq2.html). The gene expression matrix was obtained, and the proportion of genes with a zero expression value greater than 50% was removed. The input matrix was constructed using the DESeqDataSetFromMatrix function, and the difference analysis was performed using the DESeq function. The differential expression matrix was extracted from the differential analysis results using the results function, and information such as the fold change (FC) and P value of each gene was obtained. P values were adjusted using the Benjamini-Hochberg method to obtain the false discovery rate (FDR). DEGs were defined as FDR < 0.05. All downstream analyses were performed using genes with significantly different expression with FDR < 0.25.⁴⁸ 

Protein-Protein Association Networks

To reveal the interaction between the genes, a functional enrichment analysis of protein-protein association networks (PPANs) was performed using the STRING database (version 11.5, https://string-db.org/). The significance of network edges was set as confidence, the minimum required interaction score was set as high confidence (0.700), the clustering method was set as K-means clustering, and the number of clusters was set as two, three, or five. 

KEGG, GO, and GSEA

To elucidate the relationship between the ELG transcripts and potential biological functions, KEGG,⁴⁹ GO,⁴⁹ and GSEA^50,51 analyses were performed as described in previous research.^42,43 To obtain the most recent gene annotations of the KEGG pathway, the obtained ELG transcripts were mapped based on the KEGG GENES Database (https://www.genome.jp/kegg/genes.html) and the GeneCards database (https://www.genecards.org/). ELGs transcripts were mapped to the background set, and KEGG enrichment analysis was performed using the R package clusterProfiler (version 3.14.3). To understand the GO-BP of ELG transcripts, a subset of the data was also subjected to GO analysis. Candidate genes were mapped to each term in the GO database and the number of genes in each term was calculated. P adjusted (P value adjusted by the Benjamini and Hochberg method) for significance was set at 0.05. GSEA software (version 4.2.2, Broad Institute of MIT and Harvard, Cambridge, MA, USA) was used to further characterize the signaling pathways of ELG transcripts. The annotated reference gene sets, c5.go.bp.v7.2.symbols.gmt and c2.cp.kegg.v7.2.symbols.gmt were downloaded from MsigDB (https://www.gsea-msigdb.org/gsea/msigdb/). P < 0.05 and the absolute value of the normalized enrichment score (|NES|) >1 were considered significantly enriched for signaling pathways. 

Statistical and Other Analyses

GraphPad Prism 9.0.0 was used for statistical analysis and graphing. The stats, pheatmap, gseaplot2, and ggplot packages based on R (64-bit, version 4.2.1) were used to plot the PCA, gene expression heatmaps, GSEA diagram, and volcano plots, respectively. The data were tested for normal distribution before all statistical analyses. A Student's t-test or one-way ANOVA with Bonferroni correction were used for statistical analysis. A nonparametric Mann-Whitney U test was used to compare the differences between groups of data with non-normal distributions. P < 0.05 indicated a statistically significant difference, whereas NS (not significant) indicated no statistical difference. 

Results

Aging Alters Basic Physiological Behaviors in Mice

To investigate possible changes in the general physiological characteristics of aging mice, we conducted a comprehensive comparison of locomotor activity, core body temperature, and food and water consumption between young and aging mice. Our results clearly show that the activity level of total daily and 12-h interval locomotor activity in aging mice was significantly lower than that of young mice (Figs. 2A–C). In addition, the core body temperature of aging mice was significantly lower than that of young mice (Figs. 2D, 2E). Interestingly, we also found that food and water consumption were significantly higher in aging mice compared to those in young mice (Figs. 2F, 2G). These findings highlight potential alterations in the physiological characteristics of aging mice. 

Figure 2.

 

Aging alters basic physiological behaviors in mice. (A–C) Comparison of activity between young and aged mice over a light-dark cycle (A) Line graph showing the measurement of locomotor activity every five minutes for both groups of mice. (B) Bar graph showing the area under the curve (AUC) for total daily and (C) every 12-hour locomotor activity with nonparametric analysis to compare differences between groups (n = 6, *P < 0.05, **P < 0.01). (D) Curve showing the average measurement of core body temperature in young and aged mice every 20 minutes for both groups of mice over a light/dark cycle. (E) Bar graph showing the average measurement of daily core body temperature, with nonparametric analysis used to compare differences between groups (n = 6, ****P < 0.001). (F, G) Comparison of food intake (F) and water consumption (G) between young and aged mice, respectively (n = 6, **P < 0.01, ***P < 0.001).

Aging Alters the Basic Structure and Cellular Composition of Mouse ELGs

To investigate the gross structural changes of the aging ELGs used in this study, we collected ELGs from young and aging mice and performed histological observations. Compared with the lacrimal glands of young mice (Fig. 3A, left), several significant dark areas appeared in the sections of aging lacrimal glands (Fig. 3B, left). Local magnification showed that the dark areas consisted mainly of numerous cells clustered around the ducts (Figs. 3A, right, 3B, right). The main characteristic of aged lacrimal glands is chronic inflammatory changes characterized by lymphocytic infiltration.^20,52,53 To understand these features, we first compared the changes in the number of CD3⁺ lymphocytes (Figs. 3C, 3D), CD19⁺ B cells (Figs. 3E, 3F), and NKp46⁺ NK cells (Figs. 3G, 3H) under the microscope. The results showed that the number of these cells in the lacrimal glands of the aged group was significantly higher than that of the young group. In the following sections, to further validate the detailed characteristics of age-related changes in the lacrimal gland cell population and potential cell populations of significantly activated signaling pathways, we collected ELGs from young (n = 3) and aged (n = 3) mice and performed a preliminary analysis of the cell population composition under scRNA-seq. Based on classical marker genes, the major cell types in mouse ELGs were identified as belonging to 11 categories (Fig. 3I), including glandular epithelial cells, fibroblasts, myoepithelial cells (MECs), vascular endothelial cells, pericytes, and various immune cells (classical T cells, γδ T cells, B cells, NK cells, macrophages, and dendritic cells) as shown by t-distributed stochastic neighbor embedding (t-SNE) visualization (Supplementary Figs. S1A, S1B; Supplementary Table S1). Compared with the young group (Fig. 3J), the aging group showed a significant increase in classical lymphocytes and γδ T cells but a significant decrease in the number of glandular epithelial cells, fibroblasts, and myoepithelial cells (Fig. 3K). Taken together, these results suggest that aging alters the basic structure and cellular composition of the mouse lacrimal gland. 

Figure 3.

 

Aging alters the basic histological structure and cellular composition of the mouse lacrimal gland. (A, B) Gross sections of the young and aging mouse ELGs under low magnification microscopy (scale bar: 500 μm; left) and a locally magnified image of the periglandular duct (scale bar: 100 μm; right). (C, D) Distribution (C) and quantitative analysis (D) of CD3⁺ cells in young and aging ELGs, visualized by immunohistochemical staining and microscopy (scale bar: 50 μm; n = 6, *P < 0.05). (E, F) Distribution (E) and quantitative analysis (F) of CD19⁺ cells in young and aging ELGs (scale bar: 50μm; n = 6, **P < 0.01). (G, H) Distribution (G) and quantitative analysis (H) of NKp46⁺ cells in young and aging ELGs (scale bar: 50 μm; n = 6, **P < 0.01). (I) The t-SNE plot shows the distribution of cell subpopulations in the pooled dataset of young (n = 3) and senescent (n = 3) ELG samples, based on classical marker genes. (J, K) The t-SNE plots show the composites of different cell types in young ELGs (J, left) and aging ELGs (K, left). The pie charts show the number of cells in each category as a percentage of the total number of cells in young ELGs (J, right) and aging ELGs (K, right). The colors in the graph correspond to the different cell types.

Aging Alters the Basic Composition of the Diurnal Transcriptome of Mouse ELGs

We investigated the diurnal effects of aging on the transcriptome profiles of ELGs by collecting transcriptomic data from eight time points at three-hour intervals in a 24-hour cycle as a whole unit to analyze the altered transcriptome composition of ELGs in the young and aging mice. The results of PCA were showed in Figure 4A presents the results of the PCA. The variance of the principal component (PC)1 and PC2 accounted for 29.53% and 6.78% of the total variance, respectively, indicating that the transcriptomes of the young and aging mouse ELGs were significantly altered in their overall composition over a 24-hour cycle. To further demonstrate the differences between the two groups, we used the DESeq2 algorithm, set the threshold to 1.5, and identified 1092 DEGs between the two groups (Supplementary Table S2), of which 864 genes were upregulated, and 228 genes were downregulated. The scatter plot in Figure 4B and the heatmap in Figure 4C further depict the changes in the expression of DEGs in the two groups. The STRING database was used to visualize PPANs and reveal the interaction among DEGs between young and aged mouse ELGs (Fig. 4D). Five functional clusters (clusters 1–5) were identified. In addition, the top five significantly enriched KEGG pathways from genes in clusters 1–5 were displayed. The most significantly enriched KEGG pathways in clusters 1, 3, and 4 were found to be mainly associated with metabolic pathways, whereas those associated with cluster 2 were primarily related to neural and muscle activities. Additionally, the pathways associated with cluster 5 were found to be primarily related to adaptive immune responses, particularly Th1, Th2, and Th17 (Fig. 4D). Compared with KEGG analysis, GSEA provides a more sensitive understanding of the functions involving these genes by aggregating information from a broad set of genes.⁵⁰ Therefore GSEA was performed on these DEGs to reveal that the top 10 signaling pathways were mainly associated with inflammatory responses (P < 0.05 and |NES| > 1; Fig. 4E). In conclusion, aging alters the basic composition of the ELG transcriptome to be primarily characterized by an inflammatory response. 

Figure 4.

 

Alterations in overall transcriptome in aging mouse ELGs. (A) The PCA analysis shows the transcriptome composition of the young and the aging ELGs in a 24-hour cycle. The orange square scatter represents the proportion of the aging group, and the blue circular scatter represents the proportion of the young group. (B) The scatter plot displays the transcriptomes of the young and aging ELGs for eight time points in a 24-hour cycle. The x-axis denotes the log10 (expression of the young ELGs) values, and the y-axis denotes the log10 (expression of the aging ELGs) values. The red and blue scatter points indicate changes greater than ±1.5-fold. (C) The heatmap displays the DEGs (|FC| > 1.5 and FDR < 0.05) in the young (left) and aging (right) ELG transcriptome for eight time points in a 24-hour cycle. Red denotes high expression, and blue denotes low expression. The colored bar shows the proportion used to indicate the expression of genes in each group at the eight time points, and the expression range is normalized to ±3. (D) PPANs and functional clusters (Clusters 1–5) with the top five relevant KEGG pathways of all DEGs in the young and aging groups of mouse ELGs (P-adjusted < 0.05). (E) The top 10 signaling pathway sets of DEGs (|FC| > 1.5 and FDR < 0.25) in young and aging ELGs according to GSEA for the eight time points at three-hour intervals in a diurnal cycle (P < 0.05 and |NES| > 1).

Aging Alters the Immunity-Associated Transcriptome Profile of Mouse ELGs

Aging is mostly accompanied by a characteristic chronic inflammatory response in the lacrimal glands in both humans and rodents.^21,52 We first screened immune-related DEGs (setting a difference threshold FC of 1.5 and FDR < 0.05) at eight time points at three-hour intervals over a 24-hour cycle in young and aging ELGs to understand the underlying molecular mechanism accompanying inflammatory responses in the aging lacrimal gland. There were 116 upregulated and nine downregulated genes (Supplementary Table S3). The heatmap in Fig. 5A shows the expression levels of these DEGs, and the volcano plot in Fig. 5B shows the FDR versus FC. For in-depth analysis, we focused on the differences between the mean expression of the chemokine family DEGs in two groups at eight time points in a 24-hour cycle, as shown by the violin plot (Fig. 5C). Table and Supplementary Figures S2A to S2J further reveal the cell subpopulations expressing the above differentially expressed chemokines and chemokine receptors under t-SNE dimensionality reduction analysis. To understand the correlation between immune-related genes expressed by aging ELGs, the STRING database was used to visualize the inter-gene network among DEGs. These genes automatically formed three gene clusters (Fig. 5D). We performed KEGG enrichment analysis to gain further insight into the functions involved in these gene clusters. The results indicated that the genes in Cluster 1 was mainly enriched in pattern-recognition receptor (PRR)–related pathways, which mainly include the signaling pathways of the toll-like receptors (TLRs), nucleotide oligomerization domain–like receptors (NLRs), C-type lectin receptors, and retinoic acid-inducible gene-I–like receptors (Fig. 5D). However, the genes in Clusters 2 and 3 were mainly enriched in innate immune-related pathways, primarily comprising signaling pathways related to leukocyte transendothelial trafficking, chemotaxis, and phagocytosis, and in adaptive immune-related pathways, the latter including antigen processing and presentation, T-cell differentiation, and activation (especially Th1, Th2, and Th17), and B-cell differentiation–related pathways (Fig. 5D). These PRR, innate, and adaptive immune-related pathways were mainly enriched by upregulated DEGs. We have demonstrated the temporal characteristics of these three pathway types using concentric circle plots (Figs. 5E–5G) and rose diagrams (Figs. 5H–5J) based on the KEGG data, which reveal that the PRR-related pathways were enriched only during the light cycle (Figs. 5E, 5H), while the innate immune pathways (Figs. 5F, 5I) and the adaptive immune-related pathways (Figs. 5G, 5J) were enriched during both the light and dark phases of the diurnal cycle. To further explore the detailed aspects of inflammation-related signaling pathways, GSEA was employed to gain insight into the significance of the up- or downregulation of these pathways (Figs. 5K–5P, top). Next, t-SNE dimensionality reduction analysis was used to infer the most likely cell types that would undergo the pathway based on the featured leading genes identified by GSEA (Figs. 5K–5P, bottom). The GSEA results demonstrate that the TLR signaling pathway (Fig. 5K, top), the NLR signaling pathway (Fig. 5L, top), the NK cell–mediated cytotoxicity pathway (Fig. 5M, top), Fc gamma R-mediated phagocytosis pathway (Fig. 5N, top), T-cell receptor–related signaling pathway (Fig. 5O, top), and B-cell receptor–related pathway (Fig. 5P, top) were significantly activated in aging ELGs. The results of the scRNA-seq dimensional reduction analysis indicate that the featured leading genes (Tlr1, 2, 4, 5, 6) of the TLR signaling pathway were mainly expressed in fibroblast cells, macrophages, endothelial cells, B cells, MECs, and dendritic cells. The featured leading genes (Antxr1, Nlrp1a, Nlrp3, Nod2, and Nlrp6) of the NLR signaling pathway were expressed fibroblast cells, pericytes, macrophages, MECs, and dendritic cells. The featured leading genes (Prf1, Klra4, and Klra8) of the NK cell–mediated cytotoxicity pathway were mainly expressed in NK cells and T cells. The featured leading genes (Fcgr1, Fcgr2b, and Fcgr3) of the Fc gamma R–mediated phagocytosis pathway were mainly expressed in fibroblast cells, NK cells, B cells, macrophages and dendritic cells. The featured classical genes (Cd3d, Cd3e, and Cd3g) of the T-cell receptor signaling pathway are mainly expressed in αβ T cells and γδ T cells. The featured classical genes (Cd19, Cd79a, and Cd79b) of the B-cell receptor signaling pathway mainly come from B cells (Figs. 5K-5P, bottom). In conclusion, aging ELGs undergo a multitude of chemokine-driven inflammatory responses, mainly triggered by PRRs and led by T-cell responses; temporally, the inflammatory responses accompanying aging function almost around the clock. 

Figure 5.

 

Aging alters the immunity-related transcriptome profile of mouse ELGs. (A) The heatmap shows immunity-related DEGs (|FC| > 1.5 and FDR < 0.05) in the young and aging ELGs at eight time points at three-hour intervals in a 24-hour cycle. The expression range of genes is normalized to ±2. (B) The volcano plot shows immunity-associated DEGs (|FC| > 1.5 and FDR < 0.05) for all time points in a 24-hour cycle in the young and aging mouse ELGs. (C) The violin plots show the chemokine and chemokine receptor genes with significantly different expression means at eight time points in a 24-hour cycle in the young and aging mouse ELGs (*P < 0.05, **P < 0.01, ***P < 0.001). (D) PPANs and functional clusters (Clusters 1–3) with the relevant KEGG pathways of immunity-related genes in the aging and young groups of mouse ELGs (P-adjusted < 0.05). (E-G) The concentric circle plots separately display the diurnally rhythmic distribution of three different types of upregulated immune-associated KEGG pathways, including PRR-associated, natural immunity, and adaptive immunity. Eight sampling points in a 24-hour cycle are shown in a clockwise order; color fill indicates P - adjusted < 0.05, and no fill indicates NS. ① TLR signaling pathway, ② NLR signaling pathway, ③ C-type lectin receptor CLR signaling pathway, ④ retinoic acid-inducible gene-I–like receptor RLR signaling pathway, ⑤ chemokine signaling pathway, ⑥ hematopoietic cell lineage, ⑦ Fc epsilon RI signaling pathway, ⑧ leukocyte transendothelial migration, ⑨ NK cell–mediated cytotoxicity, ⑩ Fc gamma R–mediated phagocytosis, ⑪ T-cell receptor signaling pathway, ⑫ B-cell receptor signaling pathway, ⑬ Th17-cell differentiation, ⑭ Th1- and Th2-cell differentiation, ⑮ antigen processing and presentation. (H–J) Rose diagrams based on the immunity-related KEGG data. Clockwise, the order is ZT0 to ZT21. The length of the y-axis represents -log10 (P-adjusted). (K–P) GSEA (top) shows the significance of the immunity-related pathway enrichment in the aging ELGs (P < 0.05 and |NES| > 1). The cell type (bottom) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of the GSEA were expressed under a t-SNE by scRNA-seq.

Table.

 

Cell Types With High Expression of Chemokines and Chemokine Receptors in Aging Mouse ELGs

 

Aging Alters the Metabolism-Associated Transcriptome Profile of Mouse ELGs

To understand the changes that occur in the metabolic processes of aging mouse ELGs, we screened for differences between the total metabolism-related DEGs of young and aging mouse ELGs at eight time points in three-hour intervals in a 24-hour cycle (setting a difference threshold FC of 1.2 and FDR < 0.05). We identified 145 DEGs between the two groups, 94 upregulated and 51 downregulated (Figs. 6A, 6B, Supplementary Table S4). The STRING database was used to predict the protein-protein network between metabolism -associated genes in young and aging ELGs (Fig. 6C). The results indicate that these genes formed three gene clusters. The KEGG enrichment analysis suggests that these gene clusters were enriched to different degrees in amino acid metabolism, lipid metabolism, and biodegradation of the xenobiotics and metabolism pathways (Fig. 6C). Concentric circles (Figs. 6D, 6E) and rose diagrams (Figs. 6F–6H) were plotted to show the distribution of the above-mentioned pathways at eight time points to understand the diurnal distribution of these functional pathways. The results show that (1) in the light phase, the upregulated genes were enriched in the lipid metabolism pathway (Fig. 6D), whereas the downregulated genes were enriched in the amino acid biosynthesis, nucleotide metabolism, xenobiotics biodegradation, and other metabolism pathways (Fig. 6E); (2) in the dark phase, the upregulated genes were enriched in the amino acid metabolism pathways (Fig. 6D). Considering the role of certain lipid metabolic pathways in age-related inflammatory responses,^54,55 we specifically focused on GSEA of DEGs between young and aging ELGs for certain lipid pathways. GSEA demonstrates that linoleic acid metabolism, arachidonic acid metabolism, and ether lipid metabolism pathways were significantly activated in aging ELGs (Figs. 6I-6K, top). To further explore the cell populations in which these pathways may occur, we examined the cell populations expressing the featured leading genes in each pathway by t-SNE plot analysis. The results indicate that the linoleic acid metabolism pathway (Pla2g2d, Pla2g5, and Alox15) was mainly expressed in fibroblast cells, macrophages, and MECs; the arachidonic acid metabolism pathway (Ephx2, Gpx3, and Cbr2) was mainly expressed in fibroblast cells, endothelial cells, pericytes, macrophages, dendritic cells, and MECs; and the ether lipid metabolism pathway (Enpp2, Lpcat4, and Lpcat1) was mainly expressed in fibroblast cells, endothelial cells, and pericytes (Figs. 6I–6K, bottom). 

Figure 6.

 

Aging alters the metabolism-related transcriptome profile of mouse ELGs. (A) The heatmap illustrates metabolism-related DEGs (|FC| > 1.2 and FDR < 0.05) in young and aging ELGs at eight time points at three-hour intervals in a 24-hour cycle. The expression range of the genes is normalized to ±3. (B) The volcano plot shows metabolism-associated DEGs (|FC| > 1.2 and FDR < 0.05) for all time points in a 24-hour cycle in the young and aging ELGs. (C) PPANs and functional clusters (Clusters 1–3) with the relevant KEGG pathways of metabolism-related genes in the young and aging ELGs (P-adjusted < 0.05). (D, E) The concentric circle plot shows the diurnal distribution of the KEGG pathway enriched by the metabolism-associated upregulated (D) and downregulated (E) genes. The eight sampling points in a 24-hour cycle are shown in a clockwise order; color fill indicates P-adjusted < 0.05, and no fill indicates no significance. ① Tryptophan metabolism, ② histidine metabolism, ③ arginine and proline metabolism, ④ glycerolipid metabolism, ⑤ linoleic acid metabolism, ⑥ arachidonic acid metabolism, ⑦ glycerophospholipid metabolism, ⑧ ether lipid metabolism, ⑨ sphingolipid metabolism, ⑩ biosynthesis of amino acids, ⑪ purine metabolism, ⑫ pyrimidine metabolism, ⑬ drug metabolism–cytochrome P450. (F–H) Rose diagrams based on the metabolism-related KEGG data. Clockwise, the order is ZT0 to ZT21. The length of the y-axis represents -log10 (P-adjusted). (I–K) GSEA (top) shows the pathway enrichment of the linoleic acid metabolism, arachidonic acid metabolism, and ether lipid metabolism in the aging ELGs (P <0.05 and |NES| > 1). The cell type (bottom) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of GSEA were expressed under a t-SNE by scRNA-seq.

Aging Alters the Mitochondrial Bioenergetic Function–Associated Transcriptome Profile of Mouse ELGs

The differences in total mitochondrial function–related DEGs (with a set difference threshold FC of 1.2 and FDR < 0.05) were compared between the young and aging ELGs at eight time points at three-hour intervals in a 24-hour cycle to gain insight into mitochondrial function–related changes in the ELGs of aging mice. Sixty-four DEGs were identified, of which 53 were upregulated, and 11 were downregulated (Fig. 7A, Supplementary Table S5). The volcano map demonstrates the FDR and FC of DEGs directly involved in mitochondrial production (Fig. 7B). The STRING database was used to visualize the inter-gene network between these genes to comprehend the correlations between mitochondrial function-related genes expressed by ELGs in aging mice. These DEGs formed two gene clusters (Fig. 7C). The KEGG enrichment analysis was performed to gain further insight into the functions involved in these gene clusters. Genes in aging ELGs were enriched in phosphoinositide 3-kinase/protein kinase B (PI3K-Akt), hypoxia-inducible factor-1 (HIF-1), forkhead box O (FoxO), adenosine 5ʹ-monophosphate–activated protein kinase (AMPK), and mammalian target of rapamycin (mTOR) signaling pathways. Concentric circle plots (Fig. 7D) and rose diagrams (Fig. 7E) were plotted to learn the diurnal distribution pattern of these pathways. KEGG enrichment analysis was performed on the upregulated genes. The results indicate that the PI3K-Akt, HIF-1, and FoxO signaling pathways were enriched in the light and dark cycles; however, the AMPK and mTOR pathways were enriched only in the light phase. Because of the central role of mTOR signaling pathway in the mitochondrial nutrient-sensing signaling network,⁵⁶ we specifically focused on the activity status of the mTOR signaling pathway in aging ELGs and showed that this pathway was significantly activated (Fig. 7F, left). In addition, the t-SNE plot of the mTOR signaling pathway leading genes (Pik3cg, Akt3, and Pik3) were expressed in almost all immune cell subpopulations in the ELG (Fig. 7F, right). In conclusion, the aging ELGs exhibit mitochondrial bio-energetical dysfunction. 

Figure 7.

 

Aging alters the transcriptome profile of genes related to mitochondrial function in mouse ELGs. (A) The heatmap shows mitochondrial function–related DEGs in the young and aging ELGs adjusted for |FC| > 1.2 and FDR < 0.05 at eight time points at three-hour intervals in a 24-hour cycle. Gene expression range is normalized to ±2. (B) The volcano plot shows all DEGs (|FC| > 1.2 and FDR < 0.05) related to mitochondrial function at all time points in a 24-hour cycle in the young and aging ELGs. (C) PPANs and functional clusters (Clusters 1–2) with the relevant KEGG pathways of mitochondria-related genes in the aging and young groups of mouse ELGs (P-adjusted < 0.05). (D) The concentric circle diagram shows the diurnal distribution of upregulated mitochondrial function-associated KEGG. Eight time points in a 24-hour cycle are shown in a clockwise order; color fill indicates P-adjusted < 0.05, no color fill indicates no significance. ① PI3K-Akt signaling pathway, ② HIF-1 signaling pathway, ③ FoxO signaling pathway, ④ AMPK signaling pathway, ⑤ mTOR signaling pathway. (E) Rose diagrams based on the mitochondrial function–related KEGG data. Clockwise, the order is ZT0 to ZT21. The length of the y-axis represents −log10 (P-adjusted). (F) GSEA (left) shows enrichment of the mTOR signaling pathway in aging ELG (P < 0.05 and |NES| > 1). The cell type (right) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of the GSEA were expressed under a t-SNE by scRNA-seq.

Aging Alters the Neurotransmission-Associated Transcriptome Profile of Mouse ELGs

The lacrimal gland, an exocrine gland, is tightly controlled by the nervous system in terms of its secretory activity and processes.⁵⁷ To determine the effects of aging on lacrimal nerve density and lacrimation, we compared the differences in nerve density and secretory activity between young and aging mouse ELGs. The results demonstrate that the density of fluorescein isothiocyanate–conjugated anti-mouse beta III tubulin monoclonal antibody–labeled lacrimal nerves was significantly reduced in ELGs from aging mice compared to that of the lacrimal nerves from young mice (Figs. 8A, 8B). In addition, measurements of tear secretion in mice receiving pilocarpine stimulation with phenol red cotton wool indicated that stimulated tear secretion in aging mice was significantly reduced at ZT12 and ZT18 during a 24-hour circadian cycle compared to that of the young mice (Fig. 8C). To further understand the molecular mechanism behind this phenomenon, we analyzed the neural-related DEGs, comparing the young and aging ELGs at eight time points in a 24-hour cycle. The results revealed 49 DEGs (setting a difference threshold FC of 1.2 and FDR < 0.05), of which 37 were upregulated, and 12 were downregulated (Fig. 8D, Supplementary Table S6). The volcano map illustrates the FDR and FC of DEGs that were closely associated with neuronal activity (Fig. 8E). The STRING database was used to visualize the network relationships among these DEGs to unravel the protein-protein interactions among the neuronal activity–related genes expressed by aging mouse ELGs. The neural activity–related DEGs expressed by aging ELGs were grouped into two gene clusters (Fig. 8F). The common feature of these two gene clusters is that the enriched pathways occur mainly in the synaptic, neurotrophin, and retrograde endocannabinoid pathways. Concentric circles (Fig. 8G) and rose diagrams (Fig. 8H) were plotted to show the distribution of the above pathways at different time points to investigate the diurnal distribution pattern of these pathways. Neurotrophin and long-term depression signaling pathways were enriched only in the light cycle, while other pathways were enriched in almost the entire light/dark cycle. The GSEA indicate that the neurotrophin signaling pathway was significantly upregulated (Fig. 8I, left). To explore the potential subpopulations of cells in which this pathway may occur, the t-SNE plot of featured leading genes (Ntrk2, Ntf3, and Ngfr) in the neurotrophin signaling pathway indicates that this pathway was only expressed in fibroblast cells, pericytes, and MECs (Fig. 8I, right). Collectively, these data suggest that aging alters the transcriptome profile associated with neural activity in the lacrimal gland. 

Figure 8.

 

Aging alters the transcriptome profile associated with neural activity in mouse ELGs. (A, B) Representative images for lacrimal gland nerve fibers in cross-sections of the young and aging ELGs (anti-mouse beta III tubulin, green; scale bar: 50 μm; A) and their quantification by fluorescent staining (B; *P < 0.05, n = 6). (C) The line graph shows the amount of lacrimal secretion (measured with phenol red threads) in mice during a diurnal cycle, with gray shading indicating the dark phase (*P < 0.05, n = 6). (D) The heatmap shows the expression levels of neuronal activity–associated DEGs (|FC| > 1.2 and FDR < 0.05) in the young and aging ELGs at eight time points in a 24-hour cycle. The expression range of genes is normalized to ±2. (E) The volcano plot displaying upregulated (red) and downregulated (blue) DEGs (|FC| > 1.2 and FDR < 0.05) closely engaging in neuronal activity at eight time points in a 24-hour cycle in the young and aging ELGs. (F) PPANs and functional clusters (Clusters 1–2) with the relevant KEGG pathways of neural activity–related genes in the aging and young groups of mouse ELGs (P-adjusted < 0.05). (G) The concentric circle diagram shows the diurnal distribution of the upregulated neural-associated KEGG pathway. Eight sampling points in a 24-hour cycle are shown in a clockwise order; color fill indicates P-adjusted < 0.05, no fill indicates no significance. ① Cholinergic synapse, ② retrograde endocannabinoid signaling, ③ GABAergic synapse, ④ glutamatergic synapse, ⑤ serotonergic synapse, ⑥ neurotrophin signaling pathway, ⑦ long-term depression. (H) Rose diagrams based on the neuronal activity-related KEGG data. Clockwise, the order is ZT0 to ZT21. The length of the y-axis represents -log10 (P-adjusted). (I) GSEA (left) shows enrichment of the neurotrophin signaling pathway in the aging ELG (P < 0.05 and |NES| > 1). The cell type (right) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of the GSEA were expressed under a t-SNE by scRNA-seq.

Aging Alters the Cellular Process-Associated Transcriptome Profile of Mouse ELGs

Characteristic cellular processes can exist in any organ or tissue. Our team has recently explored the circadian rhythmicity and transcriptomic profiling of physiological functions in mouse ELGs^43,45,46; however, little is known about the cellular processes that occur at the transcriptional level in aging mouse ELGs. Therefore we screened the DEGs (with a difference threshold FC of 1.5 and FDR < 0.05) associated with the respective cellular processes of young and aging mouse ELGs at eight time points at three-hour intervals over a 24-hour cycle and performed an overall analysis using the gene classification function in the BGI online analysis tool (https://biosys.bgi.com). Fifty-eight DEGs related to cellular processes were identified, of which 48 genes were upregulated, and 10 genes were downregulated (Fig. 9A, Supplementary Table S7). The volcano maps show that the FDR and FC of the DEGs that were directly involved in cellular processes (Fig. 9B). The STRING database was used to visualize the network relationships among these DEGs to unravel the protein-protein interactions among the cellular processes-related DEGs expressed by aging mouse ELGs (Fig. 9C). The cellular processes-related genes expressed by the aging ELGs were divided into two gene clusters. The common feature of the two gene clusters is that the enriched pathways occur mainly in autophagy pathways (endocytosis, phagosome, and lysosome; Fig. 9C). The upregulated genes were enriched in the autophagy pathways (endocytosis, phagosome, lysosome, and apoptosis) during the light cycle in a 24-hour cycle (Figs. 9D, 9E). The lysosome (Fig. 9F, left) and apoptosis (Fig. 9G, left) pathways were significantly activated in aging ELGs, indicating the dysregulation of the cellular processes in the aging lacrimal gland. To further explore the potential cell populations in which the lysosome and apoptosis pathways occur, we used the leading genes of the lysosome pathways (Napsa, Ctse, and Laptm5) and apoptosis pathways (Pik3cg, Bcl2, and Akt3) expressed in the t-SNE plot of the major cell populations. The results indicate that activation of the lysosome pathway occurs in almost all immune cell populations (Fig. 9F, right), whereas the apoptosis pathway is expressed in the fibroblast, endothelial, and pericyte populations, as well as in all immune cell populations (Fig. 9G, right). 

Figure 9.

 

Aging alters the transcriptome profiles related to cellular processes in mouse ELGs. (A) The heatmap shows DEGs (|FC| > 1.5 and FDR < 0.05) associated with cellular processes in the young and aging ELGs at all time points in a 24-hour light/dark cycle. The expression range of the genes is normalized to ±3. (B) The volcano plot shows cell process-related DEGs (|FC| > 1.5 and FDR < 0.05) for the young and aging ELGs at all time points in a 24-hour light/dark cycle. (C) PPANs and functional clusters (Clusters 1–2) with the relevant KEGG pathways of the cellular process–related genes in the aging and young groups of mouse ELGs (P-adjusted < 0.05). (D) The concentric circle plot shows the diurnal distribution of the cellular process–related KEGG pathways. Eight sampling points in a 24-hour cycle are shown in a clockwise order; color fill indicates P-adjusted < 0.05, no fill indicates no significance. ① Endocytosis, ② lysosome, ③ phagosome, ④ apoptosis, ⑤ necroptosis. (E) Rose diagrams based on the cellular processes–related KEGG data. Clockwise, the order is ZT0 to ZT21. The length of the y-axis represents −log10 (P-adjusted). (F, G) GSEA (left) shows the enrichment of the lysosome (F) and apoptosis (G) pathways in the aging ELGs (P < 0.05 and |NES|>1). The cell type (right) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of the GSEA were expressed under a t-SNE by scRNA-seq.

Aging Alters the DNA Processing-Associated Transcriptome Profile of Mouse ELGs

Aging is often accompanied by abnormal changes in genetics and DNA processing.^58,59 Therefore we screened the DEGs (with a difference threshold FC of 1.5 and FDR < 0.05) associated with DNA processing in young and aging mouse ELGs at eight time points at three-hour intervals in a 24-hour circadian cycle. We identified 230 DEGs, of which 177 were upregulated, and 53 were downregulated (Fig. 10A, Supplementary Table S8). The volcano plots further show the FDR and FC of the genes involved in the DNA processing and their significance (Fig. 10B). To visualize the interactions between these DNA processing–related DEGs by PPANs, the K-means clustering approach was used to reveal the existence of three distinct gene clusters for these DEGs (Fig. 10C). The GO enrichment analysis was performed on each of the three gene clusters to learn about the BPs in which they are primarily involved. Notably, the BPs associated with the genes in Cluster 1 were mainly related to muscle development, growth, and morphogenesis, whereas those associated with the genes in Clusters 2 and 3 were primarily involved in DNA recombination pathways (Fig. 10C). Concentric circles (Fig. 10D) and rose diagrams (Fig. 10E) were plotted to observe the diurnal distribution of the above pathways. We also performed GSEA on DNA processing–related genes. The results showed that the DNA recombination pathways were significantly activated in the aging ELGs (Fig. 10F, left). Notably, the t-SNE plot of the leading genes (Nono, Tcf3, and Fus) in the DNA recombination pathway indicates that this pathway was highly expressed in almost all cell populations (Fig. 10F, right). Therefore it appears that the aging process is accompanied by abnormalities in the DNA recombination of all the constitutive cells of the lacrimal gland. 

Figure 10.

 

Aging alters the DNA processing-related transcriptome profile of mouse ELGs. (A) The heatmap demonstrates significant DNA processing–related DEGs (|FC| > 1.5 and FDR < 0.05) of ELGs in the young and aging groups at eight time points at three-hour intervals in a 24-hour cycle. The expression range of genes is normalized to ±2. (B) The volcano plot presents the distribution of DNA processing–associated DEGs (|FC| > 1.5 and FDR < 0.05) in the young and senescent ELGs at eight time points at three-hour intervals in a 24-hour cycle. (C) PPANs and functional clusters (Clusters 1–3) with the corresponding GO pathways of DNA processing–related genes in the aging and young mouse ELGs (P-adjusted < 0.05). (D) The concentric circle plot shows the diurnal temporal distribution of the GO pathways enriched by regenerating genes. Eight sampling points in a 24-hour cycle are shown clockwise; colored padding indicates significant upregulation of the pathway (P-adjusted < 0.05), no padding indicates no significance. ① Regulation of DNA-binding transcription factor activity, ② positive regulation of DNA-binding transcription factor activity, ③ DNA recombination, ④ regulation of DNA recombination, ⑤ positive regulation of DNA recombination. (E) Rose diagrams based on the regeneration-related GO data. The clockwise order is ZT0 to ZT21. The length of the y-axis represents −log10 (P-adjusted). (F) GSEA (left) shows the enrichment of DNA recombination pathway (P < 0.05 and |NES| > 1). The cell type (right) in which the pathway is most likely to occur was inferred from the cells in which the leading genes of the GSEA were expressed under a t-SNE by scRNA-seq.

Aging Alters Fibrosis-Associated Transcriptome Profile of Mouse ELGs

One of the important characteristics of aging lacrimal glands is the fibrosis of the glandular parenchyma. Consistent with previous studies,^20,60 H&E staining revealed the presence of striated scar structures in the histological sections of aging ELGs in mice (Fig. 11A). To understand the molecular features underlying this fibrosis, we identified a total of 94 DEGs associated with fibrosis in the aging lacrimal gland compared with that in the young lacrimal gland (|FC| > 1.2 and FDR < 0.05), with 82 upregulated and 12 downregulated genes identified (Supplementary Table S9). The heatmap in Figure 11B displayed the expression levels of these DEGs, whereas the volcano plot in Figure 11C displayed the FC and statistical significance of the upregulated and downregulated genes. The PPANs analysis indicates that fibrosis-associated genes can be roughly clustered into two groups, and were mainly enriched in the TNF, Notch, Wnt, and VEGF signaling pathways (Fig. 11D). The concentric circle diagram and rose plot indicate that these pathways were mainly active during the light phase (Figs. 11E, 11F). Further GSEA revealed that the Notch, Wnt, and VEGF signaling pathways were significantly activated in the aging ELGs (Figs. 11G–11I). To identify the potential cell populations where these signaling pathways were mainly active, we searched for the cell clusters in the t-SNE dimensional reduction plot of scRNA-seq data using the featured leading genes of these pathways. The results indicate that (1) the Notch signaling pathway (Notch1, Lfng, and Dtx1) occurred in almost all cell populations except for the acinar epithelial cells (Fig. 11J); (2) the Wnt signaling pathway (Fzd1, Fzd9, and Fzd10) occurred in fibroblasts, endothelial cells, pericytes, and MECs (Fig. 11K); and (3) the VEGF signaling pathway (Kdr [Vegfr2]) mainly occurred in endothelial cells (Fig. 11L). Finally, we identified collagen genes that were highly expressed in the aging lacrimal gland and closely associated with fibrosis, including Col27a1, Col23a1, Col6a5, and Col9a1 (Fig. 11M). Through t-SNE dimensionality reduction analysis, we found that (1) Col27a1 was mainly expressed in MECs; (2) Col23a1 was mainly expressed in fibroblast cells; (3) Col6a5 was mainly expressed in fibroblast cells and dendritic cells; and (4) Col9a1 was mainly expressed in fibroblast cells and MECs (Figs. 11N-11Q). In summary, aging significantly alters the transcriptional profile associated with the lacrimal gland fibrosis process. 

Figure 11.

 

Aging alters the fibrosis-related transcriptome profiles in mouse ELGs. (A) Hematoxylin and eosin staining shows striated fibrillar-like structures (arrows) in the young (top) and aging (bottom) ELGs (scale bar: 50 μm). (B) The heatmap shows DEGs associated with fibrosis in the young and aging ELGs at all time points in a 24-hour light/dark cycle (|FC| > 1.2 and FDR < 0.05). The expression range of genes is normalized to ±2. (C) The volcano plot shows fibrosis-related DEGs for the young and aging ELGs at all time points in a 24-hour light/dark cycle (|FC| > 1.2 and FDR < 0.05). (D) PPANs and functional clusters (Clusters 1–2) with the relevant KEGG pathways of the fibrosis-related genes in the aging and young groups of mouse ELGs (P-adjusted < 0.05). (E) The concentric circle plot separately displays the diurnal distribution of upregulated fibrosis-associated KEGG pathways. Eight sampling points in a 24-hour cycle are shown in a clockwise order. Color fill indicates P-adjusted < 0.05, and no fill indicates no significance. ① TNF signaling pathway, ② notch signaling pathway, ③ Wnt signaling pathway, ④ VEGF signaling pathway. (F) Rose diagrams based on the fibrosis-related KEGG data. The clockwise order is ZT0 to ZT21. The length of the y-axis represents −log10 (P-adjusted). (G–I) GSEA shows the pathway enrichment of the Notch (G), Wnt (H), and VEGF (I) signaling pathways in the aging ELGs (P < 0.01 and |NES| > 1). (J–L) The most likely cell type in the Notch (J), Wnt (K), and VEGF (L) signaling pathways was inferred from the cells in which the featured leading genes of the GSEA were expressed under a t-SNE by scRNA-seq. (M) The violin plots show the collagen genes with significantly different expression means at eight time points in a 24-hour cycle in the ELGs of young and aging mice (**P < 0.01, ***P < 0.001) (N–Q) Cells expressed by highly expressed collagen genes associated with fibrosis in aging mouse ELGs under a t-SNE by scRNA-seq.

Discussion

This study performed a systematic analysis of the RNA-seq transcriptome profile of ELGs in aging mice using a diurnal model. The study is among the first to investigate the time-dependent characteristics of activities in various aging-related signaling pathways. The findings contribute valuable insights into the mechanisms underlying aging in ELGs and provide a reference dataset for future investigations into the aging process of ELGs. 

A Senescence-Associated Secretory Phenotype in Aging ELGs

The concept of senescence has undergone significant changes in recent years.^61,62 Senescent cells, although no longer actively proliferating, continue to maintain an active metabolism and express a senescence-associated secretory phenotype (SASP) in a robust manner.⁶³ The SASP primarily involves the overproduction of proinflammatory molecules, growth factors, and enzymes that aid in the remodeling of the extracellular matrix.^64–66 The observations made in this study are consistent with previous findings about SASP and reveal a marked activation of numerous relevant signaling pathways, particularly those involved in inflammation, metabolism, mitochondrial bioenergetics, cellular processes, neural activity, DNA processing, and fibrosis. The underlying causes of these SASP-related changes are likely multifaceted, involving factors such as cell proliferation arrest, oncogene activation, and persistent DNA damage.^64,67,68 However, it is worth noting that SASP can have both positive and negative effects.⁶² In the context of the aging lacrimal gland, SASP may promote chronic inflammation and fibrosis, while also maintaining a relatively stable secretory phenotype and regenerative capacity through adaptive and compensatory over-activation of corresponding signaling pathways. Although the role of SASP in the aging process has been extensively studied in recent years,^69–71 further investigation is needed to fully elucidate its role in the aging process of the lacrimal gland and the mechanisms underlying these effects. 

Chronic Inflammation in Aging ELGs

Increasing evidence supports the notion that aging is often associated with the development of chronic, sterile, low-grade inflammation, characterized by elevated levels of inflammatory markers in the blood and chronic infiltration of inflammatory cells in various aging organs. In the context of lacrimal glands, previous clinical and experimental studies have indicated predominant pathological changes with lymphocytic infiltration at the histological and transcriptomic levels.^21,30,52 Nevertheless, the underlying mechanisms responsible for this chronic inflammatory response remain unclear. It has been suggested that elevated basal stress in the immune system, caused by self-fragmentation and the release of autoantigens, such as mitochondrial damage-associated molecular patterns, may contribute to an increased number of damaged and dead cells during aging.⁷² In alignment with this concept, our study observed the activation of multiple signaling pathways associated with pattern-recognition systems, lysosomes, and apoptosis in the aging lacrimal gland. Remarkably, this activation was found to occur primarily during the light cycle, likely related to the increased metabolic activity of the lacrimal gland during this phase, leading to the accumulation of immunostimulatory metabolites.⁷³ These metabolites, regarded as "waste" molecules, could be sensed by antigen-presenting cells, such as dendritic cells, macrophages, and B cells, thereby enhancing the activation of antigen-presenting signaling pathways, as evidenced in our study. Previous studies in mice have also highlighted the role of age-related changes in antigen-presenting cell function in triggering inflammatory responses in the lacrimal gland and ocular surface.⁷⁴ 

The investigation of cellular mechanisms underlying age-related lacrimal inflammation has revealed aberrant activation of Th1, Th17, and regulatory T cells in the lacrimal gland.^30,52,74 Consistently, our study observed the activation of many T cell–related signaling pathways, such as Th1, Th2, and Th17, in a complete temporal pattern in the aging lacrimal gland. B-cell responses were also affected in senescent individuals, becoming more inflammatory and less responsive to infections and vaccines.⁶¹ Our transcriptomic profiling of aging lacrimal glands similarly indicated the activation of B-cell receptor signaling pathway, Fcγ-mediated phagocytosis, and NK cell-mediated cytotoxic effects, indicating the complexity of inflammatory response mechanisms in aging lacrimal glands.^62,64 

Notably, targeting inflammation to mitigate the aging process has garnered interest.^65,66 Immunosuppressive interventions, such as topical rapamycin drops, have shown promise in alleviating abnormal ocular surface changes triggered by the lacrimal gland by increasing the number of regulatory T cells in the lacrimal gland.⁶⁷ In conclusion, our findings reveal a time-dependent pattern of inflammatory response in the aging lacrimal gland. Striving to maintain a delicate balance between proinflammatory and anti-inflammatory processes during the aging process may potentially forestall or delay functional decline in the lacrimal gland. Nonetheless, further controlled clinical studies are warranted to validate these observations. 

Metabolic Disturbance in Aging ELGs

Aging is known to involve reprogramming of metabolic processes, which can impact various ocular tissues such as the lens epithelium, retina, and optic nerve, leading to reduced metabolic capacity.^68–71 Pharmacological approaches that improve metabolic processes have shown potential to decelerate the senescent rate of ocular tissues. For instance, supplementation with nicotinamide adenine dinucleotide precursors has been observed to modulate mitochondrial dysfunction and prevent the development of glaucoma and retinal degeneration in aging mice.^75,76 

However, our current understanding of altered metabolic processes in aging lacrimal glands remains extremely limited. In this study, significant metabolic alterations in the aging lacrimal gland were identified at the transcriptome level. The main findings include upregulation of certain amino acid (arginine and proline, histidine, phenylalanine, and tryptophan) metabolic pathways, as well as lipid (alpha-linolenic acid, arachidonic acid, ether lipid, and glycerolipid) metabolic pathways. Excessive production of linoleic acid, arachidonic acid, and ether lipid may potentially act as proinflammatory lipids, triggering the inflammatory response in the aging lacrimal gland through various links.^77–80 

Additionally, decreased metabolism of purines and pyrimidines, which are essential components for nucleic acid synthesis, may be associated with impaired DNA replication and cell repair in the aging lacrimal gland. Furthermore, decreased drug metabolism and synthesis of enzymes that degrade toxic substances (e.g., cytochrome P450) could potentially alter the composition of lacrimal gland secretions. Last, a decrease in retinol, amino acid biosynthesis, and fatty acid metabolic signaling pathways may also impact the composition of lacrimal gland secretory components. In conclusion, the findings from this study may provide insights into the altered biochemical composition of the lacrimal film secreted by the aging lacrimal gland.^81,82 However, it is essential to acknowledge that further research and investigation are required to fully understand the functional implications of these metabolic alterations and their precise role in the aging process of the lacrimal gland. 

Mitochondrial Dysregulation in Aging ELGs

Mitochondrial dysfunction is considered an important marker of cellular senescence and is associated with the aging of ocular tissues.^83–88 Changes in mitochondria during aging include a decrease in their number, alterations in membrane permeability, increased oxidation of structural and functional proteins, and instability of electron transport chain complexes. Imbalances in mitochondrial intermediate metabolism and electron transport may lead to altered energy production and increased production of reactive oxygen species, resulting in oxidative stress and damage to ocular tissues. Targeting mitochondrial function has shown promise in ameliorating the aging process in ocular tissue, highlighting the significance of mitochondria in delaying the aging process in the lacrimal gland.^89,90 In the current study, the DEGs in the lacrimal gland of aging mice were found to be enriched in various signaling pathways, including PI3K-Akt, HIF-1, AMPK, FoxO, and mTOR. Activation of the PI3K/Akt pathway may increase reactive oxygen species production through nicotinamide adenine dinucleotide phosphate expression. HIF and AMPK primarily respond to energy deficiency and hypoxic stresses, acting as evolutionarily preserved survival mechanisms. The FoxO pathway plays a role in aging by regulating longevity downstream of insulin and insulin-like growth factor signaling. The mTOR pathway integrates signals from growth factors, nutrient supply, energy status, and stressors, and significant upregulation of the mTOR signaling pathway was observed in senescent lacrimal cells. This finding aligns with previous studies using the mTOR inhibitor rapamycin to reduce the inflammatory response in the aging mouse lacrimal gland and improve clinical symptoms of dry eye.⁶⁷ It is important to note that there are direct and indirect interactions between these signaling pathways. For example, the PI3K/Akt axis can suppress the transcriptional activity of FoxO in metabolic organs. In conclusion, significant mitochondrial dysfunction is evident in the aging lacrimal gland, highlighting the potential value of measures targeting dysfunctional mitochondrial rejuvenation to improve the function of aging lacrimal glands. However, further research and investigations are needed to fully understand the functional implications of these signaling pathways in the aging process of the lacrimal gland and to explore potential interventions for mitigating age-related changes. 

Aberrant Transcriptomic Profile Associated With Neural Activity in Aging ELGs

The secretion process of the lacrimal gland is tightly regulated by the nervous system. Various stimulus signals from the ocular surface, including mechanical, chemical, and temperature signals, stimulate the lacrimal nucleus in the superior salivary nucleus within the central nervous system through afferent sensory signals distributed in the cornea and conjunctiva. The lacrimal nucleus, in turn, activates lacrimal secretion through sensory and autonomic nerves that innervate the lacrimal gland. This constitutes the lacrimal secretion reflex arc, which maintains lacrimal secretion under normal or stimulated conditions.⁹¹ The downregulation of lacrimal gland function during aging may have multiple factors involved. Based on the lacrimal secretion reflex arc described above, possible causes of altered lacrimation in aging animals include the following: (1) a significant decrease in the corneal sensory nerve density, reducing the intensity of afferent stimuli^92–94; (2) changes in the lacrimal nucleus in the brain, accompanied by significant structural and functional changes in aging brain tissue^95,96; (3) a decrease in the density of nerves in the efferent arc that regulates lacrimal secretion^25,30; and (4) changes in the secretory structure and function of the lacrimal gland itself.^30,97 

In the present study, the highly expressed DEGs associated with neural activity in the aging lacrimal gland mainly enriched the activation of multiple synaptic and neurotrophin signaling pathways. This may seem contrary to the findings of neuronal loss and decreased synapse number and function reported in previous research^30,97 and the present study. However, certain experimental studies have shown that lacrimal gland secretion is not reduced or even higher in aging mice compared to young mice.^16,98 The normal lacrimal secretion process is accomplished by the contraction of myoepithelial cells (MECs) around the secretory epithelium under the control of autonomic nerves.⁹⁹ Notably, the present study observed overactivation of signaling pathways associated with muscle contraction. This could potentially be a compensatory response to the reduced number of MECs in the aging lacrimal gland, although further validation experiments are required to confirm this. Therefore it is speculated that the mouse lacrimal gland may undergo adaptive compensation at some stage of aging to maintain a normal level of ocular surface function despite dysfunction.⁹⁹ 

Altered DNA Processing in Aging ELGs

During the process of chronological aging, both exogenous and endogenous factors contribute to continuous DNA damage. In active mammalian cells, DNA damage can occur up to 105 times per day.¹⁰⁰ As age increases, the rate of DNA damage accelerates and accumulates, leading to senescent genomic instability.¹⁰¹ Genomic instability manifests as a variety of DNA alterations, including point mutations, deletions, insertions, chromosomal rearrangements, and structural changes throughout the chromosome. These alterations irreversibly modify the informational content of the genome, eventually resulting in the loss of cellular function. Consequently, genomic instability is recognized as a hallmark sign of aging.¹⁰² The study's data supports the observation that many pathways associated with DNA processing, such as DNA replication, repair, and recombination, exhibit abnormal activation during aging. These processes can lead to the production of waste or toxic substances, potentially influencing cell fate.^103–105 

Age-Related Lacrimal Gland Fibrosis in Aging ELGs

Fibrosis of the lacrimal gland structure in humans and rodents typically initiates in middle age and is characterized by the deposition of excessive collagen and other extracellular matrix proteins in the lacrimal gland parenchyma, particularly around the ducts.^25,97 This fibrotic process can lead to dilation, tortuosity, and potential obstruction of the secretory ducts.^22,106 Fibrosis of the lacrimal gland is closely associated with the infiltration of inflammatory cells, especially lymphocytes.⁹⁷ The study not only confirms the aforementioned findings but also identifies significant signaling pathways associated with fibrosis, including the Notch signaling pathway in almost all cell groups except acinar epithelial cells, the Wnt signaling pathway in fibroblasts, endothelial cells, pericytes, and MECs, and the VEGF signaling pathway in endothelial cells. Recent studies have also implicated these pathways in fibrosis of other organs and tissues through various mechanisms.^107–109 Additionally, the study highlights that significantly upregulated collagen-related genes are mainly expressed in fibroblasts and MECs. Lastly, the data from the study suggest that the activation of fibrosis-related signaling pathways primarily occurs during the light cycle. These findings may hold potential for identifying targets that could potentially aid in preventing the fibrotic process associated with aging in the lacrimal gland. However, it is important to acknowledge that further research is required to fully comprehend the underlying mechanisms involved and to explore potential therapeutic approaches targeting these signaling pathways to prevent lacrimal gland fibrosis in aging. 

Time-Dependence of Transcriptome Activities Associated With Aging Lacrimal Glands

In this study, the time-dependent transcriptome activities associated with aging lacrimal glands were investigated. The findings indicate that the majority of the entire output genes of aging lacrimal glands, along with their enriched pathways, are more active during the light period, which is equivalent to nighttime in humans. These pathways include those related to pattern recognition, lipid metabolism, mTOR, neurotrophin signaling, autophagy, and DNA processing. This observation aligns with previous studies on circadian transcription genes in the mouse lacrimal gland,⁴⁵ which also demonstrated higher activity during the light phase, likely because of the nocturnal behavior of mice and their corresponding physiological processes. These findings may have potential implications for the treatment of age-related dry eye in humans. The data suggest that appropriate pharmacological interventions targeting the relevant pathways may have a greater therapeutic effect in aging dry eye patients before they enter the later hours of the day. For instance, the administration of cyclosporine A, an immunosuppressant approved by the Food and Drug Administration for dry eye in humans, could be a viable option. 

Limitations of This Study

The present study has several limitations that should be considered when interpreting the data. First, although significant changes in the transcriptomic profiles related to inflammatory response, mitochondrial function, neural activity, and DNA processing in the aging lacrimal gland were demonstrated, the causal relationship between these factors during lacrimal gland aging remains unclear. It is important to recognize that this study provides a snapshot of a light/dark cycle and does not provide information on the dynamic changes during different stages of life.¹¹⁰ Second, it is essential to acknowledge that only male aging animals were selected as study subjects because of sexual dimorphism in lacrimal glands.^111,112 As a result, the data obtained in this study do not represent the changes that may occur in the lacrimal glands of female aging mice. Third, this study primarily focused on a bioinformatic analysis of the transcriptome data. It is crucial to note that the translation of the transcriptome into the proteome is subject to complex and delicate regulation. Finally, in the scRNA-seq validation part of this study, the efficiency of dissociation may have been compromised because of the unique cellular morphological and functional properties of secretory epithelial cells. This could result in a loss of glandular epithelial cells and introduce potential bias in downstream analysis. Therefore optimizing the dissociation protocol and assessing the quality of single cell suspensions is crucial to ensure the validity and reliability of scRNA-seq data in future studies. 

Future Directions

In this study, a temporal characteristic chronic inflammatory response, abnormal metabolism, increased neurotransmission-related molecules, and abnormal DNA processing were observed in the aging lacrimal gland. However, the specific initiating relationships between these factors remain unclear. Identifying potential measures to delay age-related decline in lacrimal gland function by manipulating certain key factors to regulate downstream processes is an avenue worth exploring in future research. Understanding the key initiating mechanisms for the onset and development of aging processes in the lacrimal gland may hold promise for effectively decelerating or halting the aging process. 

Summary

The current study presents a diurnal RNA-seq transcriptome dataset of aging lacrimal glands and examines the transcriptional profile differences between young and aging lacrimal glands in a comprehensive manner. This dataset serves as a valuable resource for investigating the physiological processes and molecular mechanisms involved in lacrimal gland aging, particularly from a time perspective. Through the analysis of this dataset, researchers can potentially gain insights into new mechanisms of ontogenesis and disease pathways, as well as further understand existing mechanistic hypotheses. Moreover, the dataset provides potential new directions for exploring the mechanisms underlying lacrimal gland aging and the translational implications of these findings. 

Acknowledgments

Supported by the Basic Science Project of Henan Eye Institute/Henan Eye Hospital (grant number 21JCZD001; Z.L.), the National Natural Science Foundation of China (grant numbers 82171014, 81770962, Z.L.; 82101089, S.H.), Natural Science Foundation of Henan Province grant 212300410169 (D.Q.), Henan Provincial Medical Science and Technology Research Joint Co-construction Project grant LHGJ20190821 (D.Q.), and Basic Science Project for Youth of Henan Eye Institute/Henan Eye Hospital grants 20JCQN002 (D.Q.). 

Supplementary Data:Supplementary Table 1–10 to this article can be found online at https://iovs.arvojournals.org/. Additional data related to this article are available through NCBI's BioProject database under accession PRJNA1000710. 

Disclosure: J. Liu, None; H. Si, None; D. Huang, None; D. Lu, None; S. Zou, None; D. Qi, None; X. Pei, None; S. Huang, None; Z. Li, None