Sleep Loss Causes Dysfunction in Murine Extraorbital Lacrimal Glands.

Purpose: Sleep loss markedly affects the structure and function of the lacrimal gland and may cause ocular surface disease as a common public health problem. This study aims to investigate the circadian disturbance caused by sleep loss leading to dysfunction of extraorbital lacrimal glands (ELGs).
Methods: A mouse sleep deprivation (SD) model for sleep loss studies was built in C57BL/6J male mice. After four weeks, the ELGs were collected at three-hour intervals during a 24-hour period. The Jonckheere-Terpstra-Kendall algorithm was used to determine the composition, phase, and rhythmicity of transcriptomic profiles in ELGs. Furthermore, we compared the non-sleep-deprived and SD-treated mouse ELG (i) reactive oxygen species (ROS) by fluorescein staining, (ii) DNA damage by immunostaining for γ-H2Ax, and (iii) circadian migration of immune cells by immunostaining for CD4, CD8, γδ-TCR, CD64, and CX3CR1. Finally, we also evaluated (i) the locomotor activity and core body temperature rhythm of mice and (ii) the mass, cell size, and tear secretion of the ELGs.
Results: SD dramatically altered the composition and phase-associated functional enrichment of the circadian transcriptome, immune cell trafficking, metabolism, cell differentiation, and neural secretory activities of mouse ELGs. Additionally, SD caused the ROS accumulation and consequent DNA damage in the ELGs, and the ELG dysfunction caused by SD was irreversible. Conclusions: SD damages the structure, function, and diurnal oscillations of ELGs. These results highlight comprehensive characterization of insufficient sleep-affected ELG circadian transcriptome that may provide a new therapeutic approach to counteract the effects of SD on ELG function.

Sufficient sleep is an indispensable part of human life to maintain good health., In modern society, sleep deprivation (SD) and its complications caused by various factors have become a growing public health problem., Chronic SD of less than seven hours per day is associated with a variety of diseases such as obesity,, diabetes,, hypertension, heart disease,, and stroke. Sleep deficiency also increases the risk of ocular diseases, such as macular degeneration and glaucoma. The latest studies on the influence of SD on lacrimal gland secretion show that poor sleep in human and experimental animal models is accompanied by reduced lacrimal gland secretion of varying degrees, changes in the structure of the lacrimal glands, and the occurrence of dry eye disease.,, Although these expand our understanding of the association between sleep quality and lacrimal gland function, little known is known about the mechanisms by which disruptions in the sleep-wake cycle affect the physiological rhythms of the lacrimal gland.

Because of the daily cycles of light on earth, mammals have evolved a robust rhythm in the physiological functions of their organs and tissues., In a steady state, this rhythm is mainly regulated by the retino-hypothalamic tract. As a circadian pacemaker, the suprachiasmatic nuclei (SCN), a tiny region in the hypothalamus of the brain, sets the timing of rhythms of various peripheral organs by regulating neuronal activity, body temperature, and hormonal signals and then completes its daily physiological functions through the coordination mechanism of these organs. The lacrimal glands maintain the health of the ocular surface by secreting an aqueous layer in the tear film of the eye together with many important chemicals, peptides, and proteins to lubricate and protect the ocular surface. Similar to other glands,, its secretion process is closely controlled by the body's circadian rhythm, including tear volume,– physical properties,, and chemical composition., Any factor that interferes with the circadian rhythm of the extraorbital lacrimal glands (ELGs), such as jet lag, hyperglycemia, high fructose intake, and aging, can affect its normal physiological function. In severe cases, it may even cause dry eye disease. Normally, sleep occurs in synchrony with endogenous circadian rhythms., Disruptions of the circadian rhythm influence the sleep-wake cycle, whereas interrupting the sleep-wake cycle also affects endogenous circadian rhythms. Several studies have shown that mistimed sleep has a profound effect on the rhythm of gene expression in central and peripheral tissues. For example, sleep restriction reduces the transcription of circadian clock genes in the brain by about 80% and severely affects the transcriptome of the mouse liver. Insufficient sleep reduces circadian clock gene transcription in the human blood from 8.8% to 6.9%. As mentioned earlier, ELGs have a robust circadian rhythm as well.,, Accumulated evidence shows sleep disturbance can cause different degrees of damage to the structure and function of ELGs dependent on the time and degree of disturbance. For example, lack of sleep quickly leads to increase of tear film osmotic pressure, shortened tear film breakup time, and decreased tear secretion., This may lead to dry eye disease in healthy people.,, Despite these, the mechanisms behind the sleep-wake cycle interference that reshapes the circadian rhythm of ELGs have not been directly explored.

Circadian rhythms are relevant to almost all physiological functions in mammals. Circadian rhythm disturbance causes disruptions in body functions through complex mechanisms and increases the risk of certain diseases., Substantial evidence suggests that normal circadian rhythms are closely linked to the oxidative stress–antioxidant defense system., When normal circadian rhythms are disturbed (either genetically or environmentally altered–), this results in the accumulation of reactive species induced by oxidative stress, which severely interferes with normal immune, metabolic, and neuroendocrine functions. A recent study confirmed that reduced sleep duration can cause accumulation of reactive oxygen species (ROS) in the gut and increase the incidence of death. Importantly, a recent studies have shown that ROS accumulation is involved in various pathological changes in the lacrimal gland and its physiological secretion of tears. Given the damaging potential of ROS, it is valuable to explore the reversible or irreversible consequences of SD-induced lacrimal gland dysfunction by evaluating ROS levels and tear secretion.

To address this question, we sought to verify two complementary hypotheses by using an SD mouse model: (1) whether SD impairs the ELG function by accumulating ROS and disrupting the metabolic, immune, and neural functions; and (2) whether SD-induced ELG dysfunction is reversible or irreversible impairments. This study highlights the interplay between sleep, circadian rhythm, and lacrimal gland secretory function and provides a new perspective for understanding how SD impairs the physiological function of lacrimal gland that may trigger ocular surface disease.

Materials and Methods

Overall Experimental and Analysis Workflow

As depicted in Supplementary Figure S1, we present an overall experimental and analysis workflow of ELGs for sleep-deprived (SD) and non-sleep-deprived (non-SD) C57BL/6J mice. In this experimental model, we focus on the effect of SD on C57BL/6J mice by using phenotype profiling and transcriptional profiling. The phenotype profiling includes mouse habits (locomotor activity, core body temperature, body weight, water intake, and pellet intake) and basic physiological characteristics of ELGs (secretion, mass, cell size, and immune cell trafficking) for the SD-treated and non-SD-treated C57BL/6J mice (Supplementary Figs. S1A, S1B, S1D). At the end of the experiment, we collected the mouse ELGs every three hours according to the zeitgeber time (ZT) and used high-throughput RNA sequencing (RNA-Seq) to obtain the diurnal transcriptome data of the ELGs in an LD cycle (Supplementary Figs. S1A–C). To discover the pattern behind the massive transcriptome data, we used the Jonckheere-Terpstra-Kendall (JTK) cycling algorithm to analyze the changes in circadian rhythm genes of mouse ELGs between the SD-treated and non-SD-treated C57BL/6J mice. Further transcriptional profiling was performed, which includes Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, function interaction networks, time series clustering analysis, and functional annotation with phase set enhanced analysis (PSEA) (Supplementary Fig. S1E). The four-week treatment of mice under normal conditions after four weeks of SD (Supplementary Fig. S1F) was also performed by us. Locomotor activity, core body temperature, body weight, tear secretion, ELG mass, and ROS were collected (Supplementary Figs. S1G, S1H).

Animals and SD Protocol

All animal experiments were implemented in compliance with the guideline set by the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research, and were approved by the Henan Province People's Hospital Institutional Animal Care and Use Committee. Wild-type C57BL/6J male mice (six to eight weeks of age; Nanjing University Model Animal Institute, Nanjing, China) and CX3CR-1 mice (six to eight weeks of age; Jackson Laboratory, Bar Harbor, ME, USA) were housed in light-tight circadian chambers and had free access to food and water under an LD cycle (lights on at 07:00 a.m. and lights off at 07:00 p.m.). The mice were housed in SD devices under an LD cycle as previously described.,, After two weeks of acclimatization, the mice were divided into four groups including the non-SD group, the SD group, a forced-recovery-in-normal-sleep-condition (SD-FR-NSC) mouse model for four weeks after the four weeks of SD, and a normal-sleep-condition (SD-NSC) mouse model for four weeks after the four weeks of SD. For the SD group, the sweep bar worked every 1.5 minutes during the light cycle (ZT0–12) for four weeks. For SD-FR-NSC group, the sweep bar worked every 1.5 minutes during the light cycle (ZT0–12) for the first four weeks and sweep bar worked every 1.5 minutes during the dark cycle (ZT12–24) for the last four weeks. For the SD-NSC group, the sweep bar ran every 1.5 minutes during the light cycle (ZT0–12) for the first four weeks, and the sweep bar shut off automatically during the light cycle (ZT0–12) for the last four weeks.

Behavioral Activity Monitoring and Analysis

Locomotor activity and core body temperature were monitored by using telemetry system (model ER-4000 E-Mitter; Mini Mitter, Sunriver, OR, USA) as previously described.– After the mouse was anesthetized, the miniature transmitter (ER-4000 E-Mitter) was surgically implanted into the abdomen of the mouse. The radio-telemetry receiver (ER-4000 receivers) was used to collect telemetry signals to monitor the locomotor activity in five-minute bins and core body temperature in 20-minute bins.

Tear Secretion Measurement

The tear volume secreted by C57BL/6J male mice was assayed as previously reported.,, Briefly, mice were lightly anesthetized and injected intraperitoneally with pilocarpine hydrochloride (4.5 mg/kg in saline solution, HY-B0726; MedChemExpress, Monmouth Junction, NJ, USA). To obtain the tear volume of mice, phenol red thread (no. 30059010; Tianjin Jingming New Technological Development Co., Ltd., Tianjin, China) in the inner canthus of the eye for 20 seconds was used at ZT0, ZT6, ZT12, and ZT18. The red length of the phenol red thread infiltrated by tear fluid in the non-SD group, the SD group, SD-FR-NSC group, and SD-NSC group was recorded and analyzed.

ELGs Immunohistochemistry and Quantitative Analysis

The bilateral ELGs and intestinal tract, including small intestine and large intestine, of the non-SD group, SD group, SD-FR-NSC group, and SD-NSC group were immediately stored in 4% paraformaldehyde for immunohistochemical investigations. The immunohistochemistry of ELGs was performed as previously reported.,, Briefly, after dewaxing, antigen retrieval, blocking endogenous peroxidase, and serum blocking, the paraffin sections of the ELGs were incubated with anti-β-catenin antibody (Cat. no. GB11015; Servicebio Company, Wuhan, China), anti-CD4 antibody (Cat. no. 25229s; Cell Signaling Technology, Inc., Danvers, MA, USA), anti-CD8 antibody (Cat. no. 98941s; Cell Signaling Technology, Inc.), anti-Ki67 antibody (Cat. no. GB111141; Servicebio Company), anti-CD64 antibody (Cat no. ab140779; Abcam, Cambridge, MA, USA), anti-γδ T antibody (Cat no. 553178; RD Systems, Minneapolis, MN, USA), anti-γ-H2Ax antibody (Cat no. 7631T; Cell Signaling Technology, Inc.), anti-nitrotyrosine antibody (Cat no. MAB3248, RD Systems), and anti-beta III tubulin monoclonal antibody (Cat. no. GB12139; Servicebio Company). Then, ELG sections were placed in PBS (pH 7.4, Cat. no. G0002; Servicebio Company) and washed 3 times for 5 minutes. After ELGs sections were dried slightly, a secondary antibody corresponding to the primary antibody was dropped into the circle to cover the tissue and was incubated at room temperature for 50 minutes. The CX3CR1+ macrophages were analyzed by using CX3CR-1 mice. Light microscopy (Panoramic 250/MIDI; 3DHISTECH Ltd, Budapest, Hungary) was used to examine and image the ELGs, and quantification was performed by using ImageJ (version 1.42q; National Institutes of Health, Bethesda, MD, USA).

ROS Measurement

The mice in the non-SD group, SD group, SD-FR-NSC group, and SD-NSC group were sacrificed, and tissues (ELGs for non-SD group, SD group, SD-FR-NSC group, and SD-NSC group, small intestine, and large intestine for non-SD group and SD group) were collected and embedded in O.C.T (Tissue-Tek, Cat. no. G6059-110ML; Servicebio Company). Frozen sections were rewarmed at room temperature and the autofluorescence quencher (Cat. no. G1221; Servicebio Company) was added for five minutes. The sections were rinsed with running water for 10 minutes and were incubated at 37°C in a dark incubator for 30 minutes with ROS dye solution (Cat. no. D7008; Sigma-Aldrich, St. Louis, MO, USA). The sections were placed in PBS (pH 7.4, Cat. no. G0002; Servicebio Company) and washed three times for five minutes each time. Then DAPI staining solution (Cat. no. G1012; Servicebio Company) was added into sections and incubated for 10 minutes at room temperature in the dark. After the sections were slightly dried, they were mounted with anti-fade mounting medium (Cat. no. G1401; Servicebio Company). To examine and image the ELGs, ortho-fluorescent microscopy (NIKON ECLIPSE C1; Nikon Inc., Melville, NY, USA) with slide scanner (NIKON DS-U3; Nikon Inc.) was used, and quantification was performed using ImageJ (version 1.42q; National Institutes of Health).

Tissue Collection, RNA Extraction and RNA-Seq

Euthanasia of animals was implemented by using cervical dislocation to obtain the mouse ELGs. All ELG samples were collected in the same two-week period in July 2020 (n = 3 per sampling time point every three hours). The bilateral ELGs were immediately stored in liquid nitrogen for total RNA extraction by using Trizol reagent (Invitrogen Inc, Carlsbad, CA, USA). According to the Trizol RNA extraction protocol, the total RNA from two ELGs of the same mouse was extracted using the RNAeasy spin column kit (Qiagen, Hilden, Germany) and qualified and quantified using Nano Drop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, Waltham, MA, USA).

High-throughput mRNA sequencing was performed according to our previous reports.,, Briefly, the extracted mRNA was purified by using oligo(dT)-attached magnetic beads and fragmented into small pieces. The first-strand cDNA and second-strand cDNA were synthesized sequentially by using random hexamer-primed reverse transcription with installed A-Tailing Mix and RNA index adapters to end repair. To get the final library, the splint oligo sequence was adopted to denature and circularize the double-stranded PCR products of cDNA fragments by using PCR amplification. Then, the mRNA library construction was carried out on the BGIseq500 platform (BGI Group, Shenzhen, China). The sequencing data of mRNA library was filtered by using SOAPnuke (version 1.5.2, https://github.com/BGI-flexlab/SOAPnuke). The clean reads were mapped to the reference genome and reference coding gene set by using HISAT2 (version 2.0.4, http://www.ccb.jhu.edu/software/hisat/index.shtml) and Bowtie2 (version 2.2.5, http://bowtiebio.sourceforge.net/%20Bowtie2%20/ind-ex.shtml), respectively.

Analysis of Rhythmic Gene Expression

To seek the 24-hour cycling rhythmic genes from ELGs transcriptome of the non-SD group and SD group, we used the JTK_CYCLE algorithm as we previously described.,,, Briefly, the expression genes from ZT0, ZT3, ZT6, ZT9, ZT12, ZT15, ZT18, and ZT21 (n = 3) were imported into the JTK_CYCLE algorithm on R software package. The period, phase, amplitude, and adjusted P value of the ELGs transcriptome were calculated by using JTK_CYCLE. According to our previous research,, we defined the rhythm gene with expression ≥0.1 and adjusted P < 0.05, nonrhythmic genes with expression ≥0.1 and adjusted P ≥ 0.05, and low-expressed genes with expression <0.1 over the 24-hour period, respectively.

Functional Annotation With KEGG and PSEA

To observe changes of phenotype in non-SD group and SD group, functional enrichment analysis with KEGG (https://www.kegg.jp/) of annotated rhythmic genes were performed by Phyper (https://en.wikipedia.org/wiki/Hypergeometric_distribution). The terms and pathways with Q value < 0.05 were defined as significant levels. Phase distribution of circadian pathways over the 24-hour period were calculated by using the PSEA version 1.1 tool as previously described., KEGG gene sets (c2.cp.kegg.v6.2.symbols.gmt) was used and downloaded from the Molecular Signatures Database. Each significant pathway was set as Kuiper Q < 0.05.

Time-Series Clustering Analysis

To reveal data sets of structures in rhythmic gene expression, the soft clustering tools using the fuzzy c-means algorithm on R software package was adopted., The tool is an R package termed Mfuzz and can be found at http://itb1.biologie.hu-berlin.de/∼futschik/software/R/Mfuzz/index.html. In this study, the number of clusters in the non-SD and SD treated rhythmic genes was set as 4 and core threshold was set as 0.7 in the R package.

Protein-Protein Association Networks

To disclose the interaction between genes and genes involved in metabolic genes, immune genes, and neural genes, protein-protein association networks (PPANs) functional enrichment analysis was determined using the STRING database (version 11.0, https://string-db.org/). The interaction sources were set as experiments and databases, meaning of network edges were set as confidence, minimum required interaction score was set as highest confidence (0.900), clustering method was set as kmeans clustering, and the number of clusters was set as 4.

Transcription Factor Enrichment Analysis

To predict which transcription factors (TFs) regulated genes associated with ROS, immunity, metabolism, cell differentiation, and neural activity, we used the ChIP-X Enrichment Analysis 3 (ChEA3) tool. The integration method of MeanRank and Fisher's exact test P value is used to rank the TFs and display the top 10 TFs.

Statistics and Software

All data are shown as mean ± standard error of the mean (SEM). GraphPad Prism software (Version 8.1.0; San Diego, CA, USA) and IBM SPSS Statistics software (Version 23.0; IBM Corp, Armonk, NY, USA) were used to present the bar charts, graphs, line charts, and data analysis. The Rayleigh vectors, phase-period distribution, and phase distribution with KEGG pathway of the oscillating genes were analyzed using Oriana software (Version 4.01; Kovach Computing Services, Pentraeth, Wales, UK) and phase set enrichment analysis algorithm. To visualize complex data, the heatmaps were created in R package (64-bit, version 3.5.3). Data were analyzed by using the Student t-test or one-way ANOVA, followed by the Bonferroni correction for multiple testing. Statistical significance was set at P < 0.05.

Results

SD Altered Animal General and Circadian Behavior

To assess the effect of SD on animal general behaviors of C57BL/6J mice, we collected the change curves of pellet intake, water intake, and body weight in SD-treated mouse model for four weeks. As shown in Figures 1A through 1F, the pellet intake, water intake, and body weight of C57BL/6J mice continue to increase over time in the non-SD group and SD group. After four weeks of SD, the pellet consumption (Fig. 1A) and body weight (Fig. 1E) of SD-treated mice were significantly lower than those in the non-SD-treated mice. However, the water consumption (Fig. 1C) of SD-treated mice was significantly greater than that of non-SD-treated mice. In the normal rhythmic state, mice consume most of their food during the nighttime (ZT12–24) and only a small amount during the daytime., As Figure 1 shows, the SD treatment significantly reduced nighttime (ZT12–24) pellet intake (Fig. 1B) but increased daytime (ZT0–12) water intake (Fig. 1D).

Figure 1.

SD treatment altered the general biological and circadian rhythmic behavior of mice. (A) The change curve of pellet intake in SD-treated mice after four weeks. *P < 0.05, **P < 0.01. (B) The change in pellet intake in SD-treated mice after four weeks in an LD cycle was recorded. *P < 0.05, ##P < 0.01, &P < 0.05, NS, statistically nonsignificant. (C) The change curve of water intake in SD-treated mice after four weeks. *P < 0.05, **P < 0.01. (D) The change in water intake in SD-treated mice after four weeks in an LD cycle was recorded. *P < 0.05, &P < 0.05, (E) The change curve of body weight in SD-treated mice after four weeks. *P < 0.05, ***P < 0.001, NS implies not statistically significant. (F) The body weight of non-SD-treated mice and SD-treated mice. ***P < 0.001. (G, H) The locomotor activity of non-SD-treated and SD-treated mice in an LD cycle were recorded continuously for four weeks. (I) The locomotor activity of non-SD-treated and SD-treated mice in an LD cycle were recorded at day 28. The gray shading indicates dark cycles. (J) The locomotor activity of non-SD-treated and SD-treated mice in the light and dark cycle were further analyzed at day 28. ***P < 0.001, ### P < 0.001, &&&P < 0.001. (K) The core body temperature of non-SD-treated and SD-treated mice in an LD cycle were recorded at day 28 (non-SD-treated mouse: F = 6.330, P < 0.001; SD-treated mouse: F = 57.432, P < 0.001). The gray shading indicates dark cycles. (L) The core body temperature of non-SD-treated and SD-treated mice in the light and dark cycle were further analyzed at day 28. *P < 0.05, ***P < 0.001, ###P < 0.001, &&&P < 0.001.

To record animal circadian behavior, we collected the locomotor activity of the non-SD-treated and SD-treated mouse models in a light-dark (LD) cycle continuously for four weeks. As shown in Figure 1G, the locomotor activity of the non-SD-treated mice in the dark cycle is greater than in the light cycle as previously reported., However, the locomotor activity of SD-treated mice has an opposite rhythm: greater in the light cycle than in the dark cycle (Fig. 1H). To further visualize and compare the effects of SD on the locomotor activity and core body temperature rhythm of mice, we presented the data at day 28 of non-SD-treated and SD-treated mice in an LD cycle (Figs. 1I–L). As shown in Figures 1I to 1L, the locomotor activity and core body temperature of the non-SD group in the dark cycle is greater than that in the light cycle as previously reported, which show a normal circadian rhythm change, that is, the core body temperature of the mouse model increases with locomotor activity. However, locomotor activity of the SD group in the light cycle is significantly greater than that in the dark cycle (Figs. 1G–I), the core body temperature of the SD group in the dark cycle is significantly greater than that in the light cycle (Figs. 1K, 1L), which show an abnormal circadian rhythm change, that is, the core body temperature of the mouse not increases with locomotor activity (Figs. 1I–L). Collectively, these results demonstrate that SD dramatically changed the animal behavior.

SD Comprehensively Alters the Rhythmical Transcriptome in ELG

To analyze the effects of SD on circadian transcriptomic profiles in murine ELGs, we investigated the expression of cycling transcripts in ELGs collected at eight time points at three-hour intervals. To analyze the transcripts, the JTK_CYCLE algorithm was adopted according to previous reports.,,, This identified 2624 (16.03%) and 2989 (18.26%) circadian transcripts (Supplementary Figs. S2A, S2B) from a total of 16,373 ELG transcripts in SD-treated mice and non-SD treated mice, respectively (Fig. 2A). As shown in Figure 2A, the numbers of rhythmic genes shared in the non-SD group and SD group are 950 (20.40%); the numbers of rhythmic genes unique to the non-SD group and SD group were 2039 (43.70%) and 1674 (35.90%), respectively. The heatmap shows the ELGs rhythmic genes with significantly different expression patterns between the non-SD group (Fig. 2B) the and SD group (Fig. 2C).

Figure 2.

Changes of rhythmic transcriptome of ELGs between the non-SD group and SD group. (A) The Venn diagram showing the comparison of rhythmic transcriptome of ELGs between the non-SD group and the SD group. (B) Heatmap showing the expression changes of rhythmic genes unique to the non-SD group () and SD group () arranged in a specific order according to the non-SD group in a 24-hour cycle. The expression range of genes were normalized to ±2. (C) Heatmap showing the expression changes of rhythmic genes unique to the SD group () and non-SD group () arranged in a specific order according to the SD group in a 24-hour cycle. The expression range of genes were normalized to ±2. (D) Venn diagram revealing the comparison of transcriptomes of ELGs between the rhythmic genes in the non-SD group (cyan), non-rhythmic genes in the non-SD group (pink), and low expression genes in the non-SD group (light orange). (E) Phase analysis of oscillating rhythmic genes unique to the non-SD group () and the SD group (). Gray shades indicate dark cycles, and white indicate light cycles. (F) Phase analysis of oscillating rhythmic genes shared between the non-SD group () and the SD group (). Gray shades indicate dark cycles, and white indicate light cycles. (G) The histograms show phase distribution of oscillating rhythmic genes unique to the non-SD group () and SD group (). The normal curves were shown as black lines. (H) The histograms show phase distribution of oscillating rhythmic genes shared between the non-SD group () and SD group (). The normal curves were shown as black lines. (I) Pie charts displaying phase analysis of oscillating rhythmic genes shared between the non-SD group and SD group.

To investigate the loss of rhythm genes in the non-SD group, we calculated the shared genes between nonrhythmic genes and low-expressed genes in the SD group and the rhythm genes in the non-SD group. As shown in Figure 2D, the overlapping numbers of nonrhythmic genes and low-expressed genes in the SD group and the rhythmic genes in the non-SD group are 2028 (13.80%) and 11 (0.10%), respectively. To investigate the sources of rhythm genes in the SD group, we calculated the shared genes between nonrhythmic genes and low-expressed genes in the non-SD group and the rhythm genes in the SD group. As shown in Supplementary Figure S2E, the overlapping numbers of nonrhythmic genes and low-expressed genes in the non-SD group and the rhythmic genes in the SD group are 1671 (11.70%) and 3 (0.10%), respectively.

To further assess the effect of rhythmically transcribed genes in ELG, the Oriana software was be used to analyze the phase, period, and Rayleigh vector of the oscillating rhythmic genes unique to the non-SD group and SD group and those shared between the non-SD group and SD group. The rose diagram showing the phase of 2989 rhythmic genes of non-SD group was mainly distributed in the light cycle (ZT0–ZT10), with a mean vector (µ) of 06:12 and a length of mean vector (r) of 0.61 (Supplementary Fig. S2C). However, the rose diagram showing the phase of 2624 rhythmic genes of the SD group was mainly distributed in the juncture of the light-dark cycle (ZT7–ZT17), with a mean vector (µ) of 12:07 and a length of mean vector (r) of 0.53 (Supplementary Fig. S2D). The rose diagram showing the phase of 2039 rhythmic genes unique to the non-SD group was mainly distributed in the light cycle (ZT0–ZT9), with a mean vector (µ) of 05:34 and a length of mean vector (r) of 0.673 (Fig. 2E, left). However, the rose diagram showing the phase of 1674 rhythmic genes unique to the SD group was mainly distributed in the juncture of the light-dark cycle (ZT9–ZT16.5), with a mean vector (µ) of 12:59 and a length of mean vector (r) of 0.62 (Fig. 2E, right). Analogously, the rose diagram shows the different distribution of the phase (ZT6-ZT12 vs. ZT7.5-ZT13.5), mean vector (µ, 07:34 vs. 10:01), and mean vector (r, 0.52 vs. 0.47) for 950 rhythmic genes shared in the non-SD group (Fig. 2F, left) and the SD group (Fig. 2E, right).

To further explore the effect of SD on the rhythmic genes of ELGs, we analyzed the phase distribution of the non-SD group and SD group. As shown in Supplementary Figure S2F, the mean values of the phase distribution for rhythmic genes of the non-SD group and SD group are 7.46 and 11.49, respectively. Analogously, the histogram shows the different mean values of phase distribution (6.47 vs. 12.25) for rhythmic genes unique to the non-SD group (Fig. 2G, left) and unique to the SD group (Fig. 2G, right); and the histogram shows the different mean values of phase distribution (8.66 vs. 9.66) for rhythmic genes shared in the non-SD group (Fig. 2H, left) and unique to the SD group (Fig. 2H, right).

Meanwhile, we investigated whether the phase was shifted by SD. Among the 950 rhythmic genes shared in the non-SD group and SD group, the phase shifts were further compared. Compared with the non-SD group, only 260 (27.37%) rhythmic genes have no phase shift, and 690 (72.63%) rhythm genes of the SD group have phase shift (Fig. 2I, left). Further analysis revealed that there are 148 (21.45%) rhythm genes in the SD group with advanced phase and 542 (78.55%) rhythmic genes in the SD group with delayed phase (Fig. 2I, right). Collectively, these data suggest that SD globally alters the circadian transcriptional profiling in ELGs.

SD Reprograms the KEGG and Phase-Set Enriched Pathways in ELGs

To recognize the changes of transcriptomic profiling, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was adopted by using the rhythmic genes unique to the non-SD group and SD group, and shared between the non-SD group and SD group. The results are shown in Figures 3A, 3C, and 3D and Supplementary Tables S1 through S3, respectively. There is no KEGG pathway in the SD group as in the non-SD group (Figs. 3A, 3C). Two enriched KEGG pathways of oscillating rhythmic genes unique to the SD group and shared between the non-SD and SD group, which are metabolic pathways and pathways for biosynthesis of secondary metabolites (Figs. 3C, 3D), were found. Therefore these results suggest that the SD reprograms normal circadian pathways in mouse ELGs.

Figure 3.

KEGG pathways and phase distribution of circadian pathways identify dynamic patterns of transcriptomic activity over 24 hours and KEGG pathways for the non-SD group and SD group. (A, C, D) Top 10 enriched KEGG pathways of oscillating rhythmic genes unique to the non-SD group (A), SD group (C), and shared between the non-SD and SD group (D). The dot size as shown in the figure shows the boundary for Q < 0.05. (B, E) Phase distribution of enriched circadian pathways (Kuiper Q < 0.05) of oscillating rhythmic genes unique to the non-SD group ( and SD group (. The red dotted line on the inner circle and outside circle represent the oscillating rhythmic genes and phase distribution of KEGG pathways, respectively. The gray shading indicates dark cycles. (F–M) The distributions of temporal genes unique to the non-SD (F–I, and SD group of ELGs (J–M, . The lines indicate the dynamic patterns used to show the expression of gene over an LD cycle, with the expression range normalized to ±2. Gray shading indicates dark cycles. The top 10 significant KEGG pathways of rhythmic genes in cluster 1-4 unique to the non-SD (F–I, and SD group (J–M, with P < 0.05 were shown.

To investigate the biologically related gene sets association with temporally coordinated expression, the PSEA was used. As shown in Figure 3B, phase distribution of the enriched circadian pathways of oscillating rhythmic genes unique to the non-SD group are increased in the light phase (ZT1.5 to ZT9). However, the phase distribution of enriched circadian pathways of oscillating rhythmic genes unique to SD group are increased in the junction of the light-dark phase (ZT9 to ZT15) (Fig. 3E). For the rhythmic genes unique to the non-SD group, the enriched circadian pathways can be divided into six categories: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems (Fig. 3B). For the rhythmic genes unique to the SD group, the enriched circadian pathways can be divided into six categories as well: (1) cellular processes: lysosome (non-SD group vs. SD group, ZT5.95 vs. ZT10.58), peroxisome (ZT6.76 vs. ZT11.16), regulation of actin cytoskeleton (ZT6.30 vs. ZT12.41), adherens junction (ZT6.56 vs. ZT12.44), endocytosis (ZT5.20 vs. ZT12.64), cell cycle (ZT5.23 vs. ZT13.48), and focal adhesion (ZT5.50 vs. ZT13.48); (2) environmental information processing: cytokine-cytokine receptor interaction (ZT5.62 vs. ZT12.65), calcium signaling pathway (ZT6.04 vs. ZT12.71), JAK stat signaling pathway (ZT4.85 vs. ZT13.50), WNT signaling pathway (ZT6.16 vs. ZT13.99), and MAPK signaling pathway (ZT5.53 vs. ZT14.49); (3) genetic information processing: proteasome, spliceosome, and ribosome (ZT7.28 vs. ZT15.61); (4) human diseases: small cell lung cancer (ZT5.45 vs. ZT12.64), pathways in cancer (ZT5.05 vs. ZT13.51), Alzheimer's disease (ZT7.05 vs. ZT13.74), Huntington's disease (ZT6.94 vs. ZT13.90), and Parkinson's disease (ZT7.24 vs. ZT13.97); (5) metabolism: valine, leucine, and isoleucine degradation, lysine degradation, propanoate metabolism, pyruvate metabolism, purine metabolism (ZT6.46 vs. ZT11.98), pyrimidine metabolism, glycolysis/gluconeogenesis, and oxidative phosphorylation (ZT8.11 vs. ZT13.08); (6) organismal systems: insulin signaling pathway (ZT6.51 vs. ZT11.27), GnRH signaling pathway (ZT6.17 vs. ZT12.12), PPAR signaling pathway, neurotrophin signaling pathway (ZT6.71 vs. ZT12.68), cardiac muscle contraction, leukocyte transendothelial migration (ZT5.56 vs. ZT13.99), and axon guidance (ZT5.69 vs. ZT14.25) (Fig. 3E). These results demonstrated that SD alters the biological pathways of ELGs in both KEGG and PSEA level.

SD Alters the Cluster-Dependent Transcriptomic Map and Rhythmic Transcript-Associated KEGG in ELGs

To discriminate the difference of structural features for large gene expression datasets affected by SD, the cluster-dependent transcriptomic map was presented by using soft Mfuzz in R package. Four different oscillation clusters of rhythmic genes were discovered over a 24-hour cycle in the non-SD (Fig. 3F-I, left) and SD group (Fig. 3J–M, left). For cluster 1 (non-SD group vs. SD group, 584 vs. 357 enriched rhythmic genes), the peaks were located at different ZT times (ZT1-ZT4 vs. ZT5-ZT11), and the troughs were located at different ZT times (ZT12-ZT18 vs. ZT15-ZT21). For cluster 2 (non-SD group vs. SD group, 535 vs. 146 enriched rhythmic genes), the peaks were located at a similar ZT time (ZT6-ZT12 vs. none), and the troughs were located at different ZT times (ZT15-ZT21 vs. ZT9-ZT15). For cluster 3 (non-SD group vs. SD group, 308 vs. 466 enriched rhythmic genes), the peaks were located at different ZT times (none vs. ZT12-ZT21), and the troughs were located at different ZT times (ZT3-ZT15 vs. ZT0-ZT6). For cluster 4 (non-SD group vs. SD group, 612 vs. 705 enriched rhythmic genes), the peaks were located at different ZT times (ZT3-ZT12 vs. ZT9-ZT15), and the troughs were located at different ZT times (ZT15-ZT21 vs. none). To further explore the rhythmic genes clustered by Mfuzz of biological processes, the KEGG enrichment analysis was adopted. The significant top 10 KEGG pathways of rhythmic genes unique to the non-SD group and SD group in clusters 1 through 4 show that SD reprogrammed the KEGG pathways as shown in Figure 3F–M (right). Consequently, these results demonstrate that SD formulates a novel transcriptomic map and KEGG pathways of rhythmic genes.

SD Does Not Elicit Core Clock Desynchrony

To determine the influence of SD on the expression pattern of core clock genes in the ELGs, transcript levels of the 12 canonical core clock genes, including Nr1d1(also known as REV-ERBα), Nr1d2 (also known as REV-ERBβ), Clock, Per1, Per2, Per3, Arntl (also known as Bmal1), Cry1, Cry2, Npas2, Rora, and Rorc were determined using three biological RNA-Seq data sets for mouse ELGs at every three-hour intervals for 24 hours. On SD treatment for four weeks, the expression patterns of core clock genes in ELGs were basically kept similar to the non-SD-treated ELGs despite some variation in temporal expression profiles (Fig. 4). Thus we concluded that interruption of sleep-wake cycle does not elicit core clock desynchrony in ELGs.

Figure 4.

The expression profiles and phase distribution of core clock genes in an LD cycle for the non-SD and SD group. The temporal expression profiles and phase distribution of twelve core clock genes that includes Arntl (non-SD-treated ELGs: F = 43.268, P < 0.001; SD-treated ELGs: F = 112.159, P < 0.001), Clock (non-SD-treated ELGs: F = 3.624, P < 0.05; SD-treated ELGs: F = 8.204, P < 0.001), Cry1 (non-SD-treated ELGs: F = 21.625, P < 0.001; SD-treated ELGs: F = 21.444, P < 0.001), Cry2 (non-SD-treated ELGs: F = 21.278, P < 0.001; SD-treated ELGs: F = 12.598, P < 0.001), Npas2 (non-SD-treated ELGs: F = 25.537, P < 0.001; SD-treated ELGs: F = 10.182, P < 0.001), Per1 (non-SD-treated ELGs: F = 9.647, P < 0.001; SD-treated ELGs: F = 13.495, P < 0.001), Per2 (non-SD-treated ELGs: F = 19.999, P < 0.001; SD-treated ELGs: F = 96.995, P < 0.001), Per3 (non-SD-treated ELGs: F = 58.984, P < 0.001; SD-treated ELGs: F = 41.468, P < 0.001), Rora (non-SD-treated ELGs: F = 2.884, P < 0.05; SD-treated ELGs: F = 6.687, P < 0.01), Rorc (non-SD-treated ELGs: F = 38.391, P < 0.001; SD-treated ELGs: F = 34.644, P < 0.001), Nr1d1 (non-SD-treated ELGs: F = 70.815, P < 0.001; SD-treated ELGs: F = 135.501, P < 0.001), and Nr1d2 (non-SD-treated ELGs: F = 25.223, P < 0.001; SD-treated ELGs: F = 98.103, P < 0.001). Cyan line represents the non-SD group and red line represents the SD group. Gray shading indicates dark cycles. < 0.05, < 0.01, < 0.001.

SD Triggers DNA Damage Response by the Accumulation of ROS in ELGs

To study the effect of SD on the expression of ROS-related genes in mouse ELGs, we first visualized the changes in the expression of the non-SD group and the SD group during an LD cycle by using a heatmap. As shown in Figure 5A, the expression of ROS-related genes regulated by SD is obviously different from the non-SD group according to the same sequence of genes. Next, to investigate the levels of oxidative stress in ELGs caused by SD, the marker for ROS, dihydroethidium, and the marker for protein nitration, nitrotyrosine, were used for immunostaining. As shown in Figure 5B, the level of ROS shows circadian oscillation peaking at ZT18 in non-SD-treated mouse ELGs. However, the SD-treated mice shown disordered circadian oscillation and a higher amount of ROS compared with the non-SD group (Figs. 5C, 5D). Accumulation of ROS in the gut was also observed as previously reported (Supplementary Fig. S3). And then, to further study the levels of mouse ELGs affected by SD, protein nitration was visualized by using anti-nitrotyrosine antibody through immunohistochemistry technique. As shown in Figures 5E through 5G, increased immunoreactivity to nitrotyrosine was observed in the lacrimal glands of the SD group compared with the non-SD group. To verify the harmful effect of oxidation on DNA, histone H2A phosphorylation (γ-H2Av, γ-H2Ax in mammals) was adopted to detect DNA double-strand breaks. Our data shown widespread DNA damage in the SD-treated ELGs with daily oscillation peaking at ZT0 (Figs. 5H–J). Finally, to predict which TFs regulated ROS-related transcripts, the ChEA3 tool was used. As Supplementary Figure S4A shows, the top 10 TFs were enriched in the SD-treated ELGs, that is, MTF1, SPI1, IRF8, TFEC, IRF5, STAT5B, FLI1, TET2, NFE2, and ELF4, respectively. Therefore these results demonstrate that SD triggers DNA damage response by the accumulation of ROS in ELGs.

Figure 5.

Accumulation of reactive oxygen species triggered by SD between the non-SD group and SD group in mouse ELGs. (A) Heatmaps of diurnal expression for reactive oxygen-related genes between the non-SD group and SD group in mouse ELGs. The expression levels of reactive oxygen-related genes were obtained from RNA-Seq and expression range was normalized to ±3. (B) Representative images of ROS levels are shown in the ZT0, ZT6, ZT12, and ZT18 after SD between the non-SD group and SD group in mouse ELGs. Scale bar: 50 µm. (C) Quantification of ROS levels in ELGs from SD and non-SD-treated mice at six-hour intervals (non-SD: F = 10.079, P < 0.001; SD: F = 0.442, P > 0.05). Each time point shows the median and SEM of biological samples. **P < 0.01, ***P < 0.001. (D) Average quantification of ROS levels in ELGs from SD and non-SD-treated mice. ***P < 0.001. (E) Representative images of nitrotyrosine levels are shown in the ZT0 after SD between the non-SD group and SD group in mouse ELGs. Scale bar: 20 µm. (F) Quantification of nitrotyrosine levels in ELGs from SD and non-SD-treated mice at six-hour intervals (non-SD-treated mice: F = 6.140, P < 0.01; SD-treated mice: F = 40.561, P < 0.001). Each time point shows the median and SEM of biological samples. **P < 0.01, ***P < 0.001. (G) Average quantification of nitrotyrosine levels in ELGs from SD and non-SD-treated mice. ***P < 0.001. (H) Representative images of γ-H2Ax levels are shown in the ZT0, ZT6, ZT12, and ZT18 after SD between the non-SD group and SD group in mouse ELGs. Scale bar: 50 µm. (I) Quantification of γ-H2Ax levels in ELGs from SD and non-SD-treated mice at six-hour intervals (non-SD-treated mice: F = 77.252, P < 0.001; SD-treated mice: F = 49.237, P < 0.001). Each time point shows the median and SEM of biological samples. **P < 0.01, ***P < 0.001. (J) Average quantification of γ-H2Ax levels in ELGs from SD and non-SD-treated mice. ***P < 0.001.

SD Alters the Phenotypes and Immunity Transcriptome Composition in ELGs

Inflammatory cells flow daily between the bone marrow, blood, and tissues, and the number of leukocytes shows circadian oscillation., To evaluate whether immune cells in the lacrimal gland show circadian oscillations, we used immunohistochemical methods and immunofluorescence technique to label a variety of immune cells. As shown in Figure 6, T cells (CD4+, CD8+, and γδ+, Figs. 6A-C, 6D-F, 6M-O) and macrophages (CD64+ and CX3CR1+, Figs. 6G-I, 6J-L) show a diurnal oscillation pattern as previous reported. For non-SD-treated mouse ELGs, the T cells (CD4+, CD8+, and γδ+, Figs. 6A, 6D, 6M) and macrophages (CD64+ and CX3CR1+, Figs. 6G, 6J) peaked at ZT12. However, the T cells (CD4+ and CD8+) peaked at ZT24 (0) (Figs. 6A, 6D), T cells (γδ+) and macrophages (CD64+ and CX3CR1+) peaked at ZT18 (Fig. 6M, 6G, 6J). Compared with the non-SD group, average relative abundance of immune cells, including T cells and macrophages, for mouse ELGs in the SD group was higher. Collectively, SD changes the circadian oscillation of immune cells in ELGs (Fig. 6P).

Figure 6.

The effect of SD on the immune-related cell and genes in mouse ELGs. (A) Diurnal oscillations of the relative abundance of CD4+ cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 8.312, P < 0.001; SD-treated mice: F = 7.197, P < 0.001). Each time point shows the median and SEM. *P < 0.05, ***P < 0.001. (B) Average relative abundance of CD4+ cells of ELGs from non-SD group and SD group. ***P < 0.001. (C) Representative immunohistochemistry images (CD4+ cells) of mouse ELGs at ZT0 from non-SD group and SD group. The sectional view of ELG structure from the non-SD group (left) and SD group (right). Scale bar: 50 µm. (D) Diurnal oscillations of the relative abundance of CD8+ cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 18.638, P < 0.001; SD-treated mice: F = 7.446, P < 0.001). Each time point shows the median and SEM. ***P < 0.001. (E) Average relative abundance of CD8+ cells of ELGs from non-SD group and SD group. ***P < 0.001. (F) Representative immunohistochemistry images (CD8+ cells) of mouse ELGs at ZT0 from non-SD group and SD group. The sectional view of ELG structure from the non-SD group (left) and SD group (right). Scale bar: 50 µm. (G) Diurnal oscillations of the relative abundance of CD64+ cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 9.532, P < 0.001; SD-treated mice: F = 4.159, P < 0.01). Each time point shows the median and SEM. **P < 0.01, ***P < 0.001. (H) Average relative abundance of CD64+ cells of ELGs from the non-SD group and SD group. ***P < 0.001. (I) Representative immunohistochemistry images (CD64+ cells) of mouse ELGs at ZT0 from the non-SD group and SD group. The sectional view of ELG structure from the non-SD group (left) and SD group (right). Scale bar: 50 µm. (J) Diurnal oscillations of the relative abundance of CX3CR1+ cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 4.239, P < 0.01; SD-treated mice: F = 18.920, P < 0.001). Each time point shows the median and SEM. *P < 0.05, ***P < 0.001. (K) Average relative abundance of CX3CR1+ cells of ELGs from the non-SD group and SD group. ***P < 0.001. (L) Representative immunohistochemistry images (CX3CR1+ cells) of mouse ELGs at ZT0 from non-SD group and SD group. The sectional view of ELG structure from the non-SD group (left) and SD group (right). Scale bar: 50 µm. (M) Diurnal oscillations of the relative abundance of γδ+ T cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 6.625, P < 0.001; SD-treated mice: F = 4.620, P < 0.01). Each time point shows the median and SEM. **P < 0.01, ***P < 0.001. (N) Average relative abundance of γδ T+ cells of ELGs from non-SD group and SD group. ***P < 0.001. (O) Representative immunohistochemistry images (γδ T+) of mouse ELGs at ZT0 from the non-SD group and SD group. The sectional view of ELG structure from the non-SD group (left) and SD group (right). Scale bar: 50 µm. (P) Schematic diagram shows diurnal changes in the number of different immune cells in non-SD-treated mouse ELGs (left) and SD-treated mouse ELGs (right). Reproduced from Huang. (Q) Heatmaps of diurnal expression for immune-related genes between the non-SD group and SD group in mouse ELGs. The expression levels of immune-related genes were obtained from RNA-Seq and expression range of DEGs was normalized to ± 3. (R) The top 10 KEGG pathways enriched histogram of immune-related genes with Q < 0.05 were displayed. (S) The PPANs and functional clusters (clusters 1–3) with relevant KEGG pathway of immune-related genes between the non-SD group and SD group in mouse ELGs (Q < 0.05).

To further study the effect of SD on the expression of immune-related genes in mouse ELGs, we visualized gene expression changes in the non-SD group and the SD group during the LD cycle using a heatmap. As shown in Figure 6Q, there are different expression patterns in the LD cycle for the non-SD and SD group. To further observe changes of functional annotation in the non-SD group and SD group, we performed functional enrichment analysis using KEGG of annotated DEGs and obtained 21 significantly enriched KEGG pathways (Supplementary Table S4). The significantly enriched top 10 KEGG pathways are shown in Figure 6R, and the PPANs were visualized by using STRING database to disclose the relevance of immune-related genes in mouse ELGs. Three functional clusters can be found (Clusters 1–3) and the significantly enriched KEGG pathways associated with immune-related genes are presented (Fig. 6S). Finally, to predict which TFs regulated immunity-related transcripts, the ChEA3 tool was used. As Supplementary Figure S4B shows, the top 10 TFs were enriched in the SD-treated ELGs, that is, PLSCR1, BATF2, RELB, IRF7, STAT2, NFKB2, STAT3, CEBPB, BATF3, and TRAFD1, respectively. Collectively, these results demonstrate that the SD dramatically changed the immunologic processes.

SD Alters the Rhythmicity of Mass, Cell Size and Differentiation in ELGs

SD has been proven to accelerate neuronal autophagy and apoptosis, corneal epithelial apoptosis, and neuropeptide-associated apoptosis in the hippocampus. To verify whether SD will adversely affect cell differentiation in the lacrimal gland, we evaluated the mass, cell size, cell cycle, cell growth, cell proliferation, apoptosis, and cell senescence of mouse ELGs through the measurement of their transcriptome composition expression levels and immunohistochemical methods. As previously reported,, our data indicated that the mouse ELG mass adjusted per body weight exhibits diurnal changes using covariance analysis (AWT). However, as shown in Figures 7A and 7B, the AWT in SD-treated ELGs significantly decreased at five timepoints and the peak shifted six hours in advance. Consistent with this, the overall longitudinal histological view of the SD-treated ELGs were shorter in size at ZT18 compared with that of non-SD-treated ELGs (Fig. 7C).

Figure 7.

The variation of cell differentiation-related genes and metabolism-related genes between the non-SD group and SD group in mouse ELGs. (A) The diurnal changes of mouse AWT of non-SD-treated and SD-treated mouse models in the LD cycle were recorded at day 28. The gray shading indicates dark cycles (non-SD-treated ELGs: F = 3.520, P < 0.05; SD-treated ELGs: F = 1.501, P > 0.05). **P < 0.01, ***P < 0.001. (B) The changes of mouse AWT of non-SD-treated and SD-treated mice in an LD cycle were recorded at day 28. ***P < 0.001. (C) Overall longitudinal view of the murine lacrimal glands from the non-SD group and SD group at ZT18. Scale bar: 500 µm. (D) The circadian oscillating pattern of ELG cell size for the non-SD group and SD group in the LD cycle (non-SD-treated mice: F = 19.983, P < 0.001; SD-treated mice: F = 2.125, P > 0.05). **P < 0.01, ***P < 0.001. (E) The average cell size of ELG cell size for the non-SD group and SD group in the LD cycle. ***P < 0.001. (F) The representative immunohistochemical images of ELG cell size for the non-SD group and SD group at ZT18. Scale bar: 20 µm. (G, H) The enriched cell cycle-related KEGG pathways and heatmap were created by GSEA in the non-SD group and SD group. (I-M) Heatmaps of diurnal expression for cell cycle (I), cell growth (J), cell proliferation (K), cell apoptosis (L), cell senescence-related genes (M) between the non-SD group and SD group in mouse ELGs. The expression levels of cell cycle-related genes were obtained from RNA-Seq and expression range of DEGs was normalized to ± 3. (N) Diurnal oscillations of the relative abundance of Ki67+ cells by using immunohistochemistry in the ELGs for the non-SD group and SD group at six-hour intervals (non-SD-treated mice: F = 2.535, P > 0.05; SD-treated mice: F = 6.062, P < 0.01). Each time point shows the median and SEM. *P < 0.05, ***P < 0.001. (O) Average relative abundance of Ki67+ cells of ELGs from the non-SD group and SD group. ***P < 0.001. (P-Q) Representative immunohistochemistry images (Ki67+ cells) of mouse ELGs at ZT6 from the non-SD group and SD group. The sectional view of ELG structure from the non-SD group (P) and SD group (Q). Scale bar: 50 µm. (R) Heatmaps of diurnal expression for metabolism-related genes between the non-SD group and SD group in mouse ELGs. The expression levels of metabolism-related genes were obtained from RNA-Seq and expression range of DEGs was normalized to ±3. (S) The top 10 KEGG pathways enriched histogram of metabolism-related genes with Q < 0.05 were displayed. (T) The PPANs and functional clusters (Clusters 1–4) with relevant KEGG pathway of metabolism-related genes between the non-SD group and SD group in mouse ELGs (Q < 0.001).

To understand the cause of the reduction in lacrimal gland mass and volume caused by SD treatment, we calculated the alteration in lacrimal gland cell size over a 24-hour cycle using immunohistochemistry as we previously described., The results showed that the circadian rhythmicity of cell size in SD-treated ELGs basically disappeared compared to a significant diurnal change of the mean cell size in the non-SD-treated ELG (Figs. 7D, 7E).

Next, to explore the molecular basis of these changes, we investigated the transcripts related to cell cycle, cell growth, cell proliferation, apoptosis, and cell senescence of ELGs using GSEA and heatmap analysis. The GSEA result showed that the cell cycle pathway was significantly enriched in ELGs (Figs. 7G, 7H). Compared with the non-SD group, the transcriptional profiling of cell cycle (Fig. 7I), cell growth (Fig. 7J), cell proliferation (Fig. 7K), apoptosis (Fig. 7L), and cell senescence (Fig. 7M) within 24 has been completely changed in SD-treated ELGs. To predict which TFs regulated cell-differentiation-related transcripts, the ChEA3 tool was used, and the top 10 TFs were enriched in the SD-treated ELGs (Supplementary Figs. S4C–F). It is evident that nine shared TFs—that is, FOXM1, CENPA, DNMT1, E2F1, E2F7, MYBL2, PA2G4, ZNF367, and ZNF695—may regulate cell-cycle, cell-growth, cell-proliferation, and cell-apoptosis transcripts. Finally, to evaluate cell proliferation using levels of Ki67+, SD-treated ELGs changed their diurnal oscillation pattern of Ki67+, peaking at different ZT and were shown to have higher expression levels compared to the non-SD-treated ELGs (Figs. 7N-Q). Collectively, SD changes the circadian oscillation of cell differentiation in ELGs.

SD Reprograms the Composition of Metabolism-Related Transcriptome in ELGs

Circadian rhythms are intertwined with metabolic processes. To investigate whether circadian rhythm disturbance is accompanied by SD, we compared the differential expression of metabolism-related genes and their enrichment pathways between the non-SD and SD groups. As shown in Figure 7R, the expression pattern of metabolism-related genes was significantly altered in SD-treated mice compared with non-SD-treated mice. Gene enrichment analysis revealed 78 significantly enriched KEGG pathways (Supplementary Table S5). The top 10 KEGG pathways in SD-treated ELG were shown in Figure 7S. The interplay network of these metabolism-related genes was further visualized through protein-protein association networks based on STRING database. As shown in Figure 7T, four metabolism-related functional clusters were identified (clusters 1–4). Finally, to predict which TFs regulated metabolism-related transcripts, the ChEA3 tool was used. As Supplementary Figure S4G shows, the top 10 TFs were enriched in the SD-treated ELGs, that is, CHCHD3, SNAPC5, ZNF511, GTF3A, ESRRA, CENPX, PREB, MLX, HMGN3, and CEBPA, respectively. In summary, these results suggest that SD significantly alters the metabolic processes of ELGs.

SD Alters the Murine ELG Transcriptome Profile of Neural Activity-Related Genes

As previously reported,,, the circadian disruption has an adverse effect on nerve-controlled response in mouse ELGs. To investigate how SD-associated circadian disruption affects the neural activities in mouse ELGs, a secretion response experiment was performed using cholinergic agonist pilocarpine as previously reported.,, As Figures 8A to 8C show, tear secretion of the non-SD group exhibited circadian oscillation peaking at ZT18 as our previous report.,, However, SD-treated mice show disordered circadian oscillation and a lower volume of tear secretion compared with the non-SD group (Figs. 8B, 8C).

Figure 8.

The variation of nerve-related genes between the non-SD group and SD group in mouse ELGs. (A) Schematic diagram of tear secretion test in mice using phenol red threads. (B) Tear secretion of mice in an LD cycle for the non-SD group and SD group (non-SD-treated mice: F = 17.660, P < 0.001; SD-treated mice: F = 3.096, P = 0.05). ***P < 0.001. (C) Comparison of average tear secretion between the non-SD group and SD group in an LD cycle. ***P < 0.001. (D) Heatmaps of diurnal expression for nerve-related genes between the non-SD group and SD group in mouse ELGs. The expression levels of nerve-related genes were obtained from RNA-Seq and expression range of DEGs was normalized to ± 3. (E) The significantly enriched KEGG pathways enriched histogram of nerve-related genes with Q < 0.05 were displayed. (F) The PPANs and functional clusters (clusters 1–4) with relevant KEGG pathway of nerve-related genes between the non-SD group and SD group in mouse ELGs (Q < 0.05). (G-H) The enriched nerve-related KEGG pathways and heatmap were created by GSEA in the non-SD group and SD group. (I) Typical image of anti-mouse beta III tubulin (red) nerve fibers between the non-SD group and SD group in mouse ELGs. Scale bar: 50 µm. (J) Average relative abundance of beta III tubulin (red) nerve fibers of ELGs from non-SD group and SD group. *P < 0.05.

To further identify the molecular signature of SD-associated influence on tear secretion, we identified the alteration of neural activity-related transcriptome profile in SD-treated ELGs over the LD cycle. The heatmap shows that there were 45 down-regulated and 133 up-regulated genes in SD-treated ELGs compared with non-SD-treated ELGs (Fig. 8D). All the significantly enriched KEGG pathways based on these neural activity–related genes were presented in Figure 8E. To further reveal an interaction network of neural-related genes, the PPANs were visualized based on the abovementioned genes. Four distinct clusters were displayed by the Kmeans clustering method (Fig. 8F). Moreover, GSEA analysis showed that neurotrophin signaling pathway was significantly enriched in ELGs (Figs. 8G, 8H). Finally, we quantified nerve density in ELGs using fluorescein isothiocyanate–conjugated anti-mouse beta III tubulin monoclonal antibody. The results showed that the nerve density of SD-treated ELG was significantly reduced compared with non-SD-treated mice (Figs. 8I, 8J). Finally, to predict which TFs regulated metabolism-related transcripts, the ChEA3 tool was used. As Supplementary Figure S4H shows, the top 10 TFs were enriched in the SD-treated ELGs, that is, MYPOP, MYT1L, CAMTA1, SNAPC5, SCRT1, HMGN3, THYN1, KCNIP3, RORB, and ARNT2, respectively. Therefore these results demonstrate that SD dramatically alters the nerve innervation, interacting networks of neural-associated genes, and secretory response of mouse ELGs.

SD Causes Irreversible Impairment of ELG Function Accompanied by General and Circadian Behavior

To determine whether locomotor activity, core temperature, ELG mass, and tear secretion is reversible or irreversible impairment after SD, the SD-treated mice were kept in cages for four weeks under the normal sleep conditions (SD-NSC) to obtain sufficient sleep and forced recovery in normal sleep condition (SD-FR-NSC) to rapid reversal of sleep-wake cycle. As Figure 9A shows, the weight of SD-FR-NSC-treated mice did not exhibit a statistically significant change compared with that of non-SD mice; however, the weight of SD-FR-NSC-treated mice exhibited a significant increase compared with that of non-SD mice. The ELG mass was also collected. The AWT of SD-NSC and SD-FR-NSC significantly increased and exhibited different oscillating rhythms compared with the non-SD group (Figs. 9B, 9C).

Figure 9.

The irreversible impairments in general and circadian behavior and ELG function of SD-treated mice. (A) Changes in body weight of non-SD, SD, SD-NSC and SD-FR-NSC groups. non-SD vs. SD, *P < 0.05, **P < 0.01, ***P < 0.001; non-SD vs. SD-NSC, < 0.05, < 0.01, < 0.001; non-SD vs. SD-FR-NSC, $P < 0.05, $$P < 0.01, $$$P < 0.001; SD vs. SD-NSC, #P < 0.05, ##P < 0.01, ###P < 0.001;. NS indicated not statistically significant. (B) The diurnal changes of mouse AWT of non-SD, SD, SD-NSC and SD-FR-NSC in the LD cycle were recorded. The gray shading indicates dark cycles (non-SD: F = 3.520, P < 0.05; SD: F = 1.501, P > 0.05; SD-NSC: F = 1.682, P > 0.05; SD-FR-NSC: F = 1.134, P > 0.05). (C) The changes of mouse AWT of non-SD, SD, SD-NSC and SD-FR-NSC in an LD cycle were recorded. (D) The tear secretion of non-SD, SD, SD-NSC and SD-FR-NSC in an LD cycle were recorded. The gray shading indicates dark cycles (non-SD: F = 17.660, P < 0.001; SD: F = 3.096, P = 0.05; SD-NSC: F = 0.106, P > 0.05; SD-FR-NSC: F = 2.263, P > 0.05). (E) Comparison of average tear secretion of non-SD, SD, SD-NSC and SD-FR-NSC in an LD cycle were recorded. (F) The locomotor activity of non-SD, SD, SD-NSC and SD-FR-NSC in an 12h/12h LD cycle were recorded. The gray shading indicates dark cycles. (G) The core body temperature of non-SD, SD, SD-NSC and SD-FR-NSC in an LD cycle were recorded. The gray shading indicates dark cycles. (H) The locomotor activity of non-SD, SD, SD-NSC and SD-FR-NSC in the light and dark cycle were further analyzed. (I) The core body temperature of non-SD, SD, SD-NSC and SD-FR-NSC in the light and dark cycle were further analyzed. (J) Quantification of ROS levels in ELGs from non-SD, SD, SD-NSC and SD-FR-NSC at six-hour intervals (non-SD: F = 10.079, P < 0.001; SD: F = 0.442, P > 0.05; SD-NSC: F = 0.442, P > 0.05; SD-FR-NSC: F = 2.893, P < 0.05). Each time point shows the median and SEM of biological samples. (K) Average quantification of ROS levels in ELGs from non-SD, SD, SD-NSC and SD-FR-NSC. (L) Representative images of ROS levels of non-SD, SD, SD-NSC and SD-FR-NSC are shown in the ZT0, ZT6, ZT12, and ZT18 in mouse ELGs. Scale bar: 50 µm.

To further evaluate whether the damage of SD to lacrimal gland function is reversible or irreversible, we further tested the amount of tear secretion. As shown in Figures 9D and 9E, compared with the SD group, the amount of tear secretion in the SD-NSC and SD-FR-NSC groups were significantly increased, but it did not return to a normal circadian rhythm. The locomotor activity and core body temperature of SD-NSC and SD-FR-NSC were also collected and visualized, which are significantly improved compared with SD mice, but did not return to normal levels (Figs. 9F-I).

To further elucidate the mechanism of lasting adverse effect of lacrimal secretion function after SD following four weeks of spontaneous recovery or forced recovery in normal sleep conditions, we examined the level of ROS in mouse ELGs. The ROS of the SD-NSC and SD-FR-NSC groups present a different oscillating rhythm compared with the non-SD group and SD group, but the average ROS levels are lower than SD and it did not return to a normal circadian rhythm (Figs. 9J-L). Collectively, these results show that the SD causes irreversible impairment in the secretory function of the lacrimal gland, at least for the short term.

Discussion

Clinical studies have shown that quality of sleep is strongly associated with the prevalence of dry eye., However, its incidence exhibits some bias depending on the population, region, study size, and accompanying patient mental status (e.g., during the COVID19 epidemic).,– For example, studies from Europe and Japan have shown that up to 45% of participants with dry eye reported poor sleep quality., In addition, people suffering from sleep disorders (such as sleep apnea) and who engage in shift work have a higher prevalence of dry eye than other populations., However, a direct link between these two conditions has not been well established at the cellular and molecular levels. In this study, using a mouse model of chronic SD, we found that SD causes severe disturbances in the circadian rhythm of the lacrimal gland, both at the transcriptomic level and at the physiological level. In addition, we found that SD causes accumulation of oxidative stress products in the lacrimal gland and DNA damage. Importantly, the disturbance of lacrimal gland secretion caused by SD are irreversible.

SD Changes the Rhythm Oscillation of ELGs

The circadian clock system consists of a central circadian clock and a peripheral circadian clock. Although SCN are considered the central or master clock, peripheral tissues also have circadian rhythm genes., Our previous studies show that the ELGs, as a peripheral tissue, also has circadian rhythm gene oscillations under the SCN control and is reprogramed by overnutrition caused by high fructose intake, hyperglycemia caused by streptozotocin, as well as changes in light cycle phase. In this study, we reveal that (1) SD significantly reduced ELG mass, and significantly changed the circadian oscillation pattern of ELG mass, which is consistent with the circadian oscillation pattern of the change in liver mass by the LD cycle; (2) the rhythm genes of ELGs in mice after four weeks of SD were reduced, which is consistent with that of human blood cells after suffering from sleep interruption,; (3) the phase of unique rhythm genes in the SD group was delayed by about seven hours compared with the phase of unique rhythm genes in the non-SD group, and SD changed the phases and diversity of the KEGG pathways of rhythm genes. This result is the same to blood findings in humans with one week of insufficient sleep, with noted reduced number of circadian expression genes and circadian amplitude,; (4) SD alters a cluster-dependent transcriptomic map and rhythmic transcript-associated KEGG pathway in ELGs. All these changes explain the adverse effects of SD on the transcriptional level of lacrimal function.

SD Leads to a Chronic Inflammatory Response in ELGs

Effective and adequate sleep time is intimately linked to the maintenance of a normal immune system. Disruption of the normal sleep-wake cycle will initiate dysregulation of the immune system. It has been found that SD interferes with the normal migratory pattern of certain immune cells in a bone marrow-circulation-peripheral tissue pattern by affecting the expression of certain adhesion molecules, chemokines, and their corresponding receptors on the vascular endothelium., The cumulative results of this chronic process can trigger the development of certain organ-specific, systemic diseases associated with immune dysfunction. Our previous studies have shown that a jet-lag model induced by altering the light cycle significantly alters the pattern of recruitment of immune cells to the lacrimal gland and dramatically alters the immune-related circadian transcriptome. Ultimately, a chronic sub-inflammatory state of the lacrimal gland is invoked. In line with this, similar results were obtained in the present study disrupting the normal sleep-wake cycle. Thus our data further highlight the close association between normal circadian rhythms and the local immune niche of the lacrimal gland.

SD Weakens Secretion Activity of ELGs Under Nervous System Control

Although insufficient sleep or interference with normal sleep can reduce the secretion of tears in humans and experimental animals, change the composition of tears, and even lead to dry eye disease,, the underlying mechanism is not completely clear. The secretory activity of the lacrimal gland is closely controlled by the nervous system, especially the parasympathetic nerve branch of autonomic nerves. We first performed an enrichment analysis, including KEGG and GSEA, of neural-related differentiation genes expression and found that neurotrophin signaling-related pathways were activated. Neurotrophins control many aspects of survival, development, and function of neurons in both the peripheral and the central nervous systems. SD can cause changes and even damage in the structure of the central and peripheral nerves. Therefore we compared ELG nerve density of the non-SD group and SD group. Our findings show that the nerve density of the ELGs in the SD group decreased significantly compared with the non-SD group. Neural activity mainly depends on changes in the quality and quantity of synapses formed on the surface of nerve endings and target cells. The synaptic size, strength, potentiation, transcriptome, and proteome of the nervous system have significant daily rhythm and dynamically change with sleep with synapses larger after waking up and smaller after sleep.

We found that compared with non-SD treated ELGs, different types of synaptic signaling pathways in SD-treated ELGs have undergone significant changes, especially dopaminergic synapse and cholinergic synapses. This is consistent with findings of abnormal synapses in the central nervous system after SD. In summary, all these abnormal activities of structures, circadian rhythms, and neural-related transcriptomes may be one of the main reasons for the significant decrease in ELG secretory activity caused by SD in this study. However, further studies on the molecular mechanisms of SD-altered ELG circadian dysfunction is needed.

SD Strengthens the Metabolism of ELGs

Circadian rhythm and metabolism have a bidirectional effect. Any change in circadian activities may rewrite the metabolic process. For example, BMAL1 knockout, Per2 knockout, and Cry1/Cry2 double knockout mice are prone to obesity and metabolic syndrome., The circadian rhythm disorder caused by SD has been confirmed to be closely related to the occurrence and development of metabolic diseases. After SD in healthy men, the levels of various metabolites in the blood have been shown to have a significant increase. SD causes significant changes in protein-related genes and free amino acids of tear secretion. Similarly, in this study, we found disorders in the expression of metabolism-related genes in SD-treated ELGs. The main feature is that expression levels, KEGG enrichment pathways, and PPANs of all metabolism-related genes have undergone major changes. The significance and mechanism of SD leading to disorders of metabolism in the ELGs need to be further studied.

SD Damages the Morphology and Structure of ELGs

SD can cause changes in the morphology and structure of many organs or tissues through changes in hormone levels, which may cause atrophy in organs such as the ventral prostate, muscle, spleen, and brain. Consistent with these studies, mouse ELGs suffering from SD in this study experienced a significant reduction in mass, volume, and cell size. Further analysis of transcripts related to cell cycle, growth, proliferation, senescence, and apoptosis of ELGs showed that the expression of transcripts related to these pathways was significantly increased in SD-treated mice compared with the non-SD group. The results of immunohistochemistry showed that SD changed the circadian rhythm of Ki67 expression in ELGs and significantly increased its expression, which was consistent with the increase in transcripts related to cell proliferation in ELGs. The significance and mechanism of SD leading to disorders of morphology and structure through cell differentiation in ELGs need to be further studied.

SD Does Not Affect the Core Clock System

So far, existing data suggest that the core circadian clock network in cells is a relatively stable system. In general, if the retino-hypothalamic tract system is not disturbed, the core clock machine remains stable, whether reversing the timing of feeding, under metabolic stress, or in the context of aging., This study further confirmed the conclusion that the core clock system is a relatively stable system by using the model that interferes with sleep but does not change the light cycle. In contrast, rhythmic clock genes or clock-controlled genes are sensitive to different environmental alterations such as nutritional challenges. Therefore further exploration of the mechanisms by which SD causes the decoupling of the core clock and the expression of downstream CCGs would be of potential therapeutic value in addressing abnormalities in lacrimal gland structure and function caused by SD.

SD-Induced Dysfunction of ELGs Is Linked to ROS and Oxidative Stress

In mammals, there is balance between oxidative stress and the antioxidant defense system. When this system is out of balance, it causes damage to proteins, lipids, and DNA of cells. Excessive and prolonged oxidative stress can cause diseases as atherosclerosis, Alzheimer's disease, Parkinson's disease, cancer, and even death., Similarly, several recent studies have shown that oxidative stress is also involved in the pathophysiological processes of the lacrimal gland, such as aging and various inflammatory responses, which ultimately manifest as a decrease in lacrimal gland secretory function. The present study found that SD and its disturbance of the circadian rhythm of ELGs were accompanied by the accumulation of oxidative products in the lacrimal gland. Furthermore, four weeks of treatment under normal sleep conditions did not significantly reduce ROS accumulation in the lacrimal glands of SD-treated mice, and lacrimal gland secretion did not return to normal levels. Although a small number of studies have shown a close correlation between the circadian system and the oxidative stress-antioxidative stress system, it is worthwhile to further explore how SD-triggered circadian rhythm disturbances contribute to the oxidative stress in the lacrimal gland and its impairment. One possibility is that an internal time lag is created during the decoupling between the core clock and the downstream CCGs, leading to oxidative stress. This deserves further exploration.

Study Limitation

This study faces certain limitations. First, there are different degrees of confounding effects in the study that depended on the time, intensity, and method of SD implementation. This study was limited to observing adverse changes in the mouse ELG after four weeks of SD treatment. Second, this study used a nocturnal animal mouse model, whereas humans belong to the diurnal population. The opposite circadian oscillation pattern is noteworthy for the interpretation of certain data. Third, bioinformatics analysis in this study was limited to the transcriptome level. Fourth, in this study, we observed only the effect of SD on ELG function in mice and did not examine the effect of SD on corneal homeostasis, especially with respect to the corneal nerve and sensation. With prolonged corneal exposure or reduced blink reflexes of the eyelids, SD may also reduce corneal sensitivity through local inflammation-mediated effects or systemic effects in the cornea. This altered sensitivity may affect the secretory function of the lacrimal gland through the lacrimation reflex., Future studies could address this problem by implementing tarsorrhaphy and/or nonvisual means of sleep deprivation. Finally, a multi-omics study would likely provide more comprehensive information on the effects of SD in lacrimal gland function.

Summary

In conclusion, we demonstrated that SD can cause lacrimal gland dysfunction by accumulating ROS and disrupting the metabolic, immune, and neural functions, in terms of changes in the circadian transcriptome and lacrimal gland structure. These impairments hold clinical relevance, particularly in contexts in which chronic insufficient sleep is prevalent, especially in populations such as emergency service personnel, medical professionals, and new parents. In addition, we found that the above pathological changes caused by SD were irreversible (Fig. 10). These data present a comprehensive platform for future research to explore potential therapeutic strategies to alleviate abnormalities in the physiological state of the ocular surface caused by chronic insufficient sleep.

Figure 10.

Graphical summary of SD-triggered ocular surface diseases. (This figure was created using the Servier Medical ART: SMART (smart.servier.com) according to a Creative Commons Attribution 3.0.