Short-term High Fructose Intake Reprograms the Transcriptional Clock Rhythm of the Murine Extraorbital Lacrimal Gland

Abstract

Purpose: To determine whether high fructose intake (HFI) influences the daily transcriptional clock rhythms of murine extraorbital lacrimal glands (ELGs).

Methods: Timed ELGs were collected from two groups of C57BL/6J mice subjected to a 12-hour light/12-hour dark (LD) cycle for 10 days; the first group received water-only feeding and the second received water with 15% fructose. Total RNA was extracted and subjected to RNA sequencing. A JTK_CYCLE algorithm and computational software were used to determine the periodicity, rhythmicity, and amplitude of the cycling transcripts. The rhythmic transcripts from different conditions were analyzed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Results: HFI feeding caused massive remodeling of the preexisting rhythmic genes in the normal control (NC)-fed ELGs. The induced transcripts in HFI-fed mice resulted in a profound reorganization of the coordinated transcriptional oscillations and KEGG pathways. Moreover, HFI feeding significantly altered the distribution of the KEGG pathways over an LD in the ELGs. Finally, we found that the ELGs have a robust core clock machinery and HFI feeding altered amplitude and the peak phase of clock gene transcriptional oscillation in ELGs.

Conclusions: Short-term HFI reprograms the daily transcriptomic oscillation of murine ELGs. This information may deepen our understanding of the outcomes of lacrimal glands altered by nutritional challenge.

The lacrimal glands are the primary producers of the aqueous component of tear fluid. The aqueous component forms a trilaminar tear film with a superficial lipid layer secreted by the Meibomian glands and an inner mucinous layer produced mainly by conjunctival goblet cells. The aqueous component contains not only an electrolyte content similar to that of plasma but also thousands of molecules, including proteins/peptides, lipids, and small molecule metabolites secreted from the lacrimal glands.^1,2 These molecules play a key role in protecting the ocular surface from microbial infections and maintaining ocular surface integrity.^3–5 Dysfunction of lacrimal gland lacrimation, caused by systemic diseases such as Sjögren's syndrome and ocular cicatricial pemphigoid, leads to aqueous-deficient dry eye, one of the most prevalent eye diseases, and results in corneal epithelial damage, significant vision loss, and a reduction in the quality of life.^6,7 Therefore, it is necessary to understand normal lacrimation processes and their underlying mechanisms in the lacrimal glands. 

Circadian rhythms are adaptations of organismal physiology that resonate with the day/night cycle and display an endogenous, entrainable oscillation of approximately 24 hours.⁸ In mammals, these 24-hour rhythms are driven by a circadian clock. The primary circadian clock in mammals is located in a pair of distinct groups of cells in the hypothalamus called the suprachiasmatic nucleus. Information on the lengths of the day and night signals from intrinsically photosensitive retinal ganglion cells in the retina reach the suprachiasmatic nucleus through a pathway called the retinohypothalamic tract. The suprachiasmatic nucleus orchestrates the peripheral clock in different tissues located throughout the body via hormonal and neural signaling.⁹ At the cellular level, the mammalian clock consists of a “core” transcriptional feedback loop with CLOCK and BMAL1 as positive factors and period (PER) and cryptochrome (CRY) as negative regulators. These components work in concert to generate a rhythm of roughly 24 hours' duration. 

All tissues and cells contain a specific intrinsic circadian clock. Circadian rhythms from this clock drive and align a wide variety of physiological events.¹⁰ Limited research has shown that the lacrimal glands also have a biorhythm clock similar to those in other peripheral tissues. The secretions from lacrimal glands and the renewal of tears have certain circadian rhythms and are controlled by a circadian clock (Vu CHV, et al. IOVS 2018;59:ARVO E-Abstract 4930).¹¹ Humans secrete more tears during the daytime than at night.¹² At night, due to lack of local and light stimuli, especially during sleep, the rate of tear secretion is significantly reduced.¹¹ Symptoms of ocular dryness vary widely during the day and are more severe at the end of the day.¹³ Some diurnal variations in other tear film parameters, including pH,¹⁴ osmolarity,^15,16 tear volume,^12,17,18 and cytokine levels have been studied.^19,20 However, little is known of the diurnal transcriptional changes of lacrimal glands. 

Nutritional challenges, such as time-restricted feeding and high-calorie diets, interfere with circadian metabolic rhythms by mechanisms that are not completely defined.^21,22 Fructose is a natural monosaccharide broadly used in modern society. Epidemiological studies have demonstrated that high fructose intake (HFI) is one of the main etiological factors of metabolic syndrome,^23,24 including high rates of unintentional exposure in children younger than 6 years. However, there is no definitive evidence as to whether HFI feeding can also affect the circadian rhythm of normal physiologic activities, especially in lacrimal glands. 

We hypothesized that the lacrimal glands, which have an active exocrine function, are driven by a circadian rhythmic transcriptional program and that changes in the clock rhythm of transcription levels occur due to metabolic stress. To validate this hypothesis, we examined the rhythm of lacrimal transcript levels in normal-feeding mice and HFI-feeding mice by analyzing the high-throughput transcriptome of the extraorbital lacrimal glands (ELGs) and determined the effect of HFI feeding on the rhythmic expression of the lacrimal gland clock gene and downstream-controlled genes. Our data offer the potential for understanding the pathophysiology of lacrimal gland disorders and developing new strategies to treat lacrimal dysfunction diseases. 

Materials and Methods

Animals

Forty-eight male wild-type C57BL/6J mice were housed under a 12-hour light/12-hour dark (LD) cycle. The mice had ad libitum access to food and water. All procedures followed the guidelines described in the ARVO Statement for the Use of Animals in Vision and Ophthalmic Research and were approved by the Henan Province People's Hospital Institutional Animal Care and Use Committee. Throughout this study, time is indicated using the zeitgeber time (ZT) scale as the indicator for the phase of the rhythm, whereby ZT0 refers to the time that the lights went on (6 AM), and ZT12 refers to the time that the lights went off (6 PM). 

HFI Feeding Protocol

After a 2-week acclimatization to an LD schedule, the mice were randomly divided into two groups (Fig. 1A). For the HFI experiments, mice were given a standard pellet diet with sterile tap water supplemented with 15% D (−)-fructose (>99.5% purity; cat. no. F3510; Sigma-Aldrich, St. Louis, MO, USA) at 8 weeks of age.^25,26 The age-matched mice in the control group were fed a standard pellet diet with sterile tap water as the normal control (NC). The mice's body weight, pellet intake weight, and fluid intake weight were measured every 2 days. 

Figure 1

 

Animal feeding protocol and changes in body weight, water consumption, and food consumption after feeding initiation. (A) Schematic diagram of the HFI-feeding regimen. (B) Consumption of food in male C57BL/6J mice during the 10 days after the start of the feeding experiment in NC-fed and HFI-fed mice. (C) Consumption of drinking water in male C57BL/6J mice during the 10 days after the start of the feeding experiment in NC-fed and HFI-fed mice. (D) Kinetics of weight gain in male C57BL/6J mice 10 days after starting the feeding experiment in NC-fed and HFI-fed mice. Data pooled from three independent experiments. Error bars represent the mean ± SEM.

Tissue Sample Collection and RNA Extraction

After 10 days of the above feeding regimen, the mice were killed by cervical dislocation, and ELGs were collected every 3 hours over the circadian cycle as previously described, with modifications.^27,28 The samples were rapidly frozen on dry ice. RNA was extracted from two pooled ELGs from one animal for each sample using a Trizol RNA extraction protocol followed by cleanup using the RNAeasy spin column kit (Qiagen, Hilden, Germany). To avoid the effects of seasonal changes, all the sample collections were completed within the same 2 weeks in January. 

RNA Sequencing (RNA-Seq)

The first step in the workflow involved purifying the poly-A-containing mRNA molecules using poly-T oligo-attached magnetic beads. Following purification, the mRNA was fragmented into small pieces using divalent cations at an elevated temperature. The cleaved RNA fragments were copied into first-strand cDNA using reverse transcriptase and random primers. This step was followed by second-strand cDNA synthesis using DNA Polymerase I and RNase H. These cDNA fragments then underwent the addition of a single A base and the subsequent ligation of the adapter. The products were then purified and enriched with PCR amplification. We then quantified the PCR yield with a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and pooled the samples together to make a single-strand DNA (ssDNA) circle, which gave the final library. DNA nanoballs were generated with the ssDNA circle by rolling circle replication to enlarge the fluorescent signals in the sequencing process. The DNA nanoballs (DNBs) were loaded into the patterned nanoarrays, and single-end reads of 50 base pairs were read through on the BGISEQ-500 platform for the following data analysis study. For this step, the BGISEQ-500 platform combined the DNA nanoball-based nanoarrays and stepwise sequencing using the Combinational Probe-Anchor Synthesis Sequencing Method. 

Analysis of Rhythmic Gene Expression

To determine the NC- and HFI-specific gene transcripts, the time-ordered fragments per kilobase of transcript per million mapped reads (FPKM) of actively transcribed genes were triplicated and inputted to the Jonckheere-Terpstra-Kendall (JTK)_CYCLE algorithm, an efficient nonparametric algorithm for detecting rhythmic components in genome-scale datasets and a function in the R package, as an eight-time-point dataset.²⁹ This allowed a thorough identification of the oscillation patterns starting at each different time point. Oscillating transcripts were defined as those with a JTK_CYCLE P ≤ 0.05 and an oscillation period within the range of 24 hours.^21,30 Based on the P value from the output of JTK_CYCLE, we determined the circadian transcription in both the NC-fed and HIF-fed mice. 

To intuitively visualize the oscillometric characteristics of the rhythmic genes, the three parameters (period, phase, and amplitude) acquired from the JTK_CYCLE calculation and online software (http://zuotu.91maths.com) were used to establish a mathematical model of the sinusoidal curve. The sine curve formula used in this study was as follows:  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}Y = AMP\cdot\cos \left( {{{2\pi } \over {PER}}\left( {x - LAG} \right)} \right) + n\end{equation}

where Y represents the y value corresponding to the function on the same rectangular coordinate system, and AMP, cos, PER, x, and LAG represent the amplitude, cosine, period, ZT time, and phase lag, respectively.  

Enrichment Analysis for Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Annotation

To assign the rhythmic genes to specific biological pathways, KEGG pathway annotation was performed using blastall software (National Center for Biotechnology Information [NCBI], Bethesda, MD, USA) against the KEGG database V81.0 (http://www.genome.jp/kegg/) and NCBI RefSeq GCF_000001635.25_GRCm38.p5 as the reference gene set. A value of Q < 0.05, which corrected for the false discovery rate, was used as a cutoff for significantly enriched terms.³¹ 

Statistical Analysis and Software

GraphPad Software (GraphPad Prism 6; La Jolla, CA, USA) was used for bar, scatter, line charts, and statistical analyses. Oriana software (Version 4.01; Kovach Computing Services, Pentraeth, Wales, UK) was used for analyzing the phase, period distribution, and Rayleigh vector of the oscillating gene. Power BI Desktop (Microsoft, Redmond, WA, USA) was used for visualizing proportional composition histograms. The Venn Diagram Plotter (https://omics.pnl.gov/software/venn-diagram-plotter) was used to compare the numbers of rhythmic genes in different groups. Heatmaps were generated by Pheatmap scripts in R (64-bit, version 3.43). All the data were tested for normal distribution. All the normally distributed data are displayed as the number of genes unless otherwise noted. Comparisons between the two groups were performed using an unpaired Student's t-test with Bonferroni correction for multiple testing. P < 0.05 was considered statistically significant. 

Results

Alterations in Body Weight, Water Consumption, and Food Consumption After HFI Feeding

At the initiation of feeding, NC-fed and HFI-fed animals had similar body weights. HFI-fed animals showed a little more, but not significantly more, pellet intake after the initiation of feeding up to 10 days of feeding, except day 4 after feeding, compared with NC-fed animals (Fig. 1B). The fluid intake of HFI-fed animals was significantly decreased at day 2 and day 4 after the initiation of feeding compared to that of NC-fed animals (Fig. 1C). After days 6 and 10 after feeding, HFI-fed animals had significantly greater body weights compared with NC-fed animals (Fig. 1D). 

HFI Alters the Composition of the ELG Transcriptome

To globally assess whether HFI feeding affects the transcriptome oscillation of ELGs, we first performed a transcriptional time course experiment cycle every 3 hours over an LD using the ELGs of NC-fed and HFI-fed mice (Fig. 2A). Totals of 24,341 and 24,538 genes were measured from the ELGs of the NC-fed and HFI-fed mice, respectively (Figs. 2B, 2C). To determine the differences in all the transcriptomes between the two groups' ELGs, we compared the transcripts expressed in both groups. Venn diagram plotting revealed that the overlap between the total transcripts in the NC-fed group and HFI-fed group accounted for 94.54% (23,758) of the total transcripts (24,938) from both groups together; HFI-fed ELGs lost 581 genes but induced the expression of 779 new genes that were not measured in NC-fed ELGs (Fig. 2D). To further understand the differences between the two groups, we divided all the transcripts from each group's ELGs into three categories: rhythmic, not-rhythmic, and low-expressed (FPKM < 0.001) (Figs. 2B, 2C). To identify rhythmic transcripts over an LD cycle, we used a JTK_CYCLE approach (a 24-hour oscillation period and P < 0.05 cutoff) to process all the transcripts from both groups. The data showed that of the 24,341 transcripts in the ELGs of NC-fed mice, approximately 3.11% (757) exhibited rhythmic expression, 67.02% (16,313) exhibited not-rhythmic expression, and 29.87% (7271) were low-expressed (Fig. 2B). Supplementary Table S1 lists all information about oscillating gene transcripts identified by RNA-Seq and JTK_CYCLE in NC-fed ELGs. However, of the 24,538 transcripts in the ELGs of HFI-fed mice, approximately 7.14% (1752) exhibited rhythmic expression, 71.31% (17,499) exhibited not-rhythmic expression, and 21.55% (5287) were low-expressed (Fig. 2C). Supplementary Table S2 presents all the information about oscillating gene transcripts identified using RNA-Seq and JTK_CYCLE in HFI-fed ELGs. These data suggest that short-term HFI feeding altered the composition of the ELG transcriptome. 

Figure 2

 

HFI alters the transcriptome composition of mouse ELGs. (A) Experimental design and sorting of rhythmic genes in NC-fed (24,341 genes) and HFI-fed (24,538 genes) ELGs. ELGs were collected every 3 hours. Extracted RNA was processed for RNA-Seq, and the JTK algorithm was used to determine daily rhythmic genes (P < 0.05). (B, C) Pie charts of transcriptome composition in NC-fed ELGs and HFI-fed ELGs. (D) The Venn diagram representing the relationships of the total transcripts of NC-fed and HFI-fed ELGs to each other.

HFI Reprograms Global Transcriptional Rhythmic Profiling in the ELGs

To determine the effects of HFI feeding on the nature of the aforementioned rhythmic genes, we performed heatmap and Venn plotting analyses on the rhythmic genes from both the NC-fed and HFI-fed ELGs. Heatmaps displayed rhythmic genes found exclusively in the NC-fed and HFI-fed groups, respectively (Fig. 3A). Figure 3A (left) shows the expression pattern of the rhythmic genes of NC-fed ELGs, whereas, based on the arrangement of the rhythmic genes from NC-fed ELGs, the rhythmic pattern in HFI-fed ELGs disappeared. Similarly, based on the rhythmic gene pattern from HFI-fed ELGs, the rhythmic pattern in NC-fed ELGs also disappeared (Fig. 3A, right). 

Figure 3

 

HFI feeding reprograms the rhythmic transcriptome of the murine ELGs. (A) Heatmaps displaying oscillating transcripts (Q < 0.05) exclusive to the NC-fed group (left) and exclusive to the HFI-fed group (right). The color bar indicates the scale used to show the expression of gene across eight time points, with the highest expression normalized to two. (B) KEGG pathway distribution in NC-fed (red) and HFI-fed (green) ELGs. The horizontal dashed line in the figure represents the boundary for Q < 0.05. (C) Gene annotation on KEGG pathways enriched in oscillating genes with Q < 0.05 in NC-fed mice only. (D) Gene annotation on KEGG pathways enriched in oscillating genes with Q < 0.05 in HFI-fed mice only.

To evaluate the effect of HFI feeding on the KEGG pathways of the circadian cycling gene in the ELGs, we first performed KEGG enrichment analyses on the NC-fed mice's rhythmic genes and the HFI-fed mice's rhythmic genes. As shown in Figure 3B, the general pattern of pathway distribution between NC and HFI feeding changed. Although both groups showed enrichment for the core clock machinery (circadian rhythm) and a variety of different key catabolic pathways, in HFI-fed ELGs many unique KEGG pathways were induced (Figs. 3C, 3D). Supplementary Table S3 shows the complete list of the most overrepresented KEGG pathways in both groups. Taken together, these data suggest that HFI feeding results in comprehensive alterations of biological pathways in a 24-hour oscillation of murine ELGs. 

To explore the effects of HFI feeding on rhythmic transcriptomes in the ELGs, Venn plotting was performed on the rhythmic genes between the two groups. The results showed that there was a total of 2399 transcripts from the two groups with rhythmic expression. Of these, 26.97% (647 genes) had rhythm only in the NC-fed mice (Fig. 4A, left). However, in HFI-fed mice, most of these rhythmic genes that had preexisted in NC-fed mice disappeared, and 1642 new rhythmic genes specifically expressed under HFI conditions appeared, accounting for 68.45% of the total rhythmic genes in the two groups (Fig. 4A, right). In addition, the NC-fed and HFI-fed groups share 4.59% (110) of their rhythmic genes (Fig., 4A middle). As 4.59% of all rhythmic genes oscillated in both NC and HF conditions, we analyzed their phase of expression. Of these 110 genes, 14.55% oscillated in phase, whereas 85.45% were phase-shifted by HFI (Fig. 4B); 36.17% showed a phase delay, whereas 63.83% of the shifted transcripts showed a phase advance (Fig. 4B). 

Figure 4

 

KEGG analysis of 24-hour period oscillating genes after HFI feeding. (A) Venn diagram displaying the total number (top) and ratio (bottom) of rhythmic genes in the ELGs from NC-fed and HFI-fed mice. (B) Proportion of the oscillatory transcriptome shared in both NC-fed and HFI-fed ELGs sets that is phase-shifted (left), and the direction of the phase shift (right). (C) The KEGG pathway distribution in shared oscillating genes from the NC- and HFI-fed ELGs (blue), oscillating genes in only NC-fed mice (gray), and oscillating genes in only HFI-fed mice (orange). The horizontal dashed line represents the boundary for Q < 0.05. (D) Gene annotation on oscillating genes with Q < 0.05 reveals pathways that are oscillatory in both the NC- and HFI-treated ELGs. (E) Gene annotation on enriched KEGG pathways in which oscillatory expression was lost by the HFI-fed group. (F) Gene annotation on KEGG pathways represented by genes that were oscillatory only in the HFI-fed ELGs. (G) Venn diagram displaying the overlapping number of rhythmic genes in HFI-fed ELGs with low-expressed genes in the NC-fed ELGs (red) and low-rhythmic genes in the NC-fed ELGs (yellow).

To further define the functional and biological pathways of resetting the rhythmic genes, we performed separate KEGG enrichment analyses on the three sets of genes (shared, NC-specific, and HFI-specific). We found that the oscillating genes in the three groups showed their own unique annotations (Fig. 4C). First, we performed gene function annotations for rhythmic genes shared by both the NC-fed and HFI-fed mice. The results revealed that shared oscillating genes had a unique annotation, terpenoid backbone biosynthesis (Fig. 4D, Supplementary Table S4). In the NC-group–specific oscillating transcripts, most of the pathways were linked to basic biological functions and processing as a unique oscillating route, mainly including ribosome, fatty acid metabolism (Fig. 4E). Supplementary Table S4 lists the pathways of specific rhythmic genes in NC-fed mice. However, in HFI-fed ELGs, pathways were associated with infection (human T-lymphotropic virus type I infection, toxoplasmosis, and salmonella infection), one carbon pool by folate, and adipocytokine signaling pathway (Fig. 4F). Supplementary Table S4 lists the pathways of specific rhythmic genes in HFI-fed mice. Altogether, these data suggest that HFI feeding reprogramed the rhythmically expressed transcriptional profiling and biological pathways. 

To understand the possible origins of HFI feeding–induced rhythmic genes, we used Venn plotting to compare the 1642 unique rhythmic genes induced by HFI feeding to the genes in the not-rhythmic and low-expressed categories in the NC-fed ELGs. The results showed that 1552 of the HFI-induced rhythmic transcripts overlapped with the not-rhythmic transcript category of the NC-fed ELGs and 55 overlapped with the low-expressed category; however, 35 transcripts did not overlap with either of above two categories (Fig. 4G). These data indicated that HFI feeding induced new rhythmic genes from that of NC-fed ELGs. 

HFI Increases Transcriptome Profiling of Circadian Period Oscillation

Most circadian rhythm activities in mammals are 24-hour cycle activities, including various physiological activities, behaviors, transcriptomes, and proteomes. To determine the effects of HFI feeding on rhythmic profiling of the 24-hour period in the murine ELGs, we selected those rhythmically expressed genes expressed only in the 24-hour period with the JTK_CYCLE, as before, and obtained the parameters for their oscillations. The data revealed that, compared with the NC-fed group (Fig. 5A, left), the HFI-fed group (Fig. 5A, right) had more 24-hour-period rhythmic genes with different amplitudes. The number of rhythmic genes increased from 757 in the NC-fed group to 1752 in the HFI-fed group. 

Figure 5

 

HFI alters characteristics of 24-hour period oscillation genes. (A) A scatter plot showing the amplitude distribution of the oscillating genes with a 24-hour period (P < 0.05). (B, C) Simulated mathematical model diagrams of the oscillating features of rhythmic genes with approximately 5 to 6 AMP and at different phases with different timings in the NC-fed (B) and HFI-fed (C) ELGs. (D) A rose diagram and Rayleigh plots of the peak phase distributions of rhythmically expressed genes at the 24-hour period in NC-fed ELGs. The blue line indicates the mean vector (μ) representing the mean of the phases. The length of the mean vector (r) represents the degree of synchrony of the group. (E) A rose diagram and Rayleigh plots of the peak phase distributions of rhythmically expressed genes at the 24-hour period in HFI-fed ELGs. The blue line indicates the mean vector (μ) representing the mean of the phases. The length of the mean vector (r) represents the degree of synchrony of the group.

To illustrate the oscillating characteristics of all the rhythmic genes in the two groups, we selected two parameters, amplitude (Amp 5–6) and phase, according to the JTK_CYCLE, and subjected them to the mathematical model formula for plotting. Figures 5B and 5C show the oscillation characteristics of the NC group and the HFI-fed group, respectively. As expected, the oscillation patterns between NC and HFI-fed ELGs were drastically altered. After HFI feeding, oscillations increased with different phases and high amplitudes. This suggested that short-term HFI feeding rapidly changes the pattern of a superimposed oscillation. 

To understand the phase characteristics of the 24-hour periodic oscillation, the peak phase distributions were plotted on a Rayleigh graph form for NC-fed (Fig. 5D) and HFI-fed mice (Fig. 5E). The results show that the oscillating phase of the 24-hour period occurred in the light cycle (from ZT0 [24] to ZT12) in the NC-fed group; the phase range in the HFI-fed group was wider compared with the NC-fed group; a small section of the phases occurred in the dark cycle (about ZT18); the average vector magnitude (r) of the oscillating phase between the two groups decreased from 0.521 in the NC-fed group to 0.342 in the HFI-fed group; and the direction mean vector (μ) decreased from ZT7:34 in the NC-fed group (Fig. 5D) to ZT3:24 in the HFI-fed group (Fig. 5E). 

HFI Disturbs the Distribution of Rhythmic Peak Genes and KEGG Pathways During Light and Dark Cycles

We hypothesized that the rhythmic gene functions clustered under an LD cycle might be informative for understanding the circadian aspects of ELGs governed by the peripheral clock. We performed a KEGG analysis with the highest enrichment and Q < 0.05 during the light period (ZT3–ZT12) and the dark period (ZT15–ZT24) in both the NC-fed group and HFI-fed groups. As expected, the active degree (Q values) of the peak gene under the KEGG pathway analysis of the oscillating genes during the dark cycle was higher than that during the light cycle in NC-fed mice; however, only some pathways during the day cycle were higher than those during the dark cycle in NC-fed mice (Fig. 6A). This is consistent with the peak energy demand and food intake during the active phase for nocturnal rodents. 

Figure 6

 

HFI feeding influences the distribution of rhythmic peak genes and KEGG pathways during the light cycle and the dark cycle. (A) Comparisons of KEGG pathway distributions with all cycling transcripts during the light cycle (red) and dark cycle (blue) in the NC-fed group (left side). The x-axis represents the −log2 Q value, and the y-axis represents the different biological pathways. (B) KEGG pathway distributions with all cycling transcripts during light cycle only (red) and dark cycle only (blue) in HFI-fed mice. (C) Comparisons of KEGG analysis with all cycling transcripts between light cycles in only NC-fed mice (red) and only HFI-fed mice (blue) (D) Comparisons of KEGG pathways with all cycling transcripts between dark cycles in NC-fed (red) and HFI-fed (blue) mice. The vertical dashed line in the figure represents the boundary for P < 0.05.

However, compared with NC-fed ELGs, the active degree (−log2 Q values), kind, and number of pathways of peak genes under the KEGG pathway analysis during the dark cycle in HFI-fed ELGs were obviously less than those during the light cycle (Fig. 6B). Strikingly, the active degree (−log2 Q values), kind, and number of pathways of peak genes under the KEGG pathway analysis during the light cycle in HFI-fed ELGs were obviously more than those during the dark cycle (Fig. 6B). 

Next, we examined whether HFI feeding respectively influences the biological pathways of both the light cycle and dark cycle. As shown in Figures 6C and 6D and Supplementary Tables S5 and S6, the enrichment pattern of the KEGG pathways of both the light cycle and the dark cycle after HFI feeding dramatically changed. These results suggest that HFI feeding dramatically remodeled the biological pathway distribution of the ELGs over an LD cycle. 

HFI Alters the Expression of Core Clock Genes in Mouse ELGs

The above-defined effect of HFI feeding on the cycling gene expression of ELGs and the limited data on clock gene expression in rodent lacrimal glands (Lorber M. IOVS 2003;44:ARVO E-Abstract 2504) (Vu CHV, et al. IOVS 2015;56:ARVO E-Abstract 3078) prompted us to formulate the hypothesis that the ELGs house a circadian clock system and that HFI feeding affects the oscillation of the core circadian gene in the ELGs. To test this, we used RNA-Seq and the JTK_CYCLE to analyze the periods, phases, and amplitudes of the diurnal expression of the transcripts encoding 10 core components (Clock, Arntl [Bmal1], Per2, Per3, Nr1d1, Nr1d2, Rorc, Npas2, Cry1, and Cry2) of the circadian clock genes in the ELGs from the NC- and HFI-fed mice. The data revealed that, consistent with recent findings in other peripheral tissues,³² more than 10 core components of circadian clock genes are rhythmically expressed in NC-fed ELGs (Fig. 7A). Because high-calorie diets influence the behavioral and molecular circadian rhythms of the mouse liver,^21,33 we hypothesized that HFI feeding might alter the transcription of core clock components in the ELGs. However, we found that 10 core circadian gene expressions we assayed were quite similar in the NC and HFI groups (Fig. 7A). However, the Rorc gene expression shifted backward by 3 hours; also, the amplitudes of the oscillation of the genes Per2 and Per3 were dampened at the peak time point, whereas the amplitude of oscillation of the genes Arntl and Npas2 were enhanced at the nadir points (Fig. 7A). 

Figure 7

 

Core clock gene transcription in the ELGs after HFI feeding. Transcripts of the core circadian clock genes (Clock, Arntl, Per2, Per3, Nr1d1, Nr1d2, Rorc, Npas2, Cry1, and Cry2) in the murine ELGs were analyzed using RNA-Seq data and JTK_CYCLE. ELGs were harvested every 3 hours from either NC- (red lines) or HFI-fed mice (blue lines) for 10 days. (A) The temporal abundance of the transcripts of nine core clock genes in an LD cycle. The x-axis represents the sampled time points; the y-axis represents the percent of the expression of a particular gene at a certain ZT time point, accounting for the total expression at all ZT time points. Gray shading indicates dark cycles. (Mean ± SEM, n = 3, each time point; Student's t-test with Bonferroni correction, *P < 0.05 between groups at each time point.) (B) The peak phase distribution in ZT for core clock genes.

Recent studies have found that phase peaks of clock gene expression in different organs and tissues occur at different ZT times.³² To characterize the phase change of the core clock genes, we plotted the data obtained based on the JTK_CYCLE analysis using Oriana software (Kovach Computing Services). As Figure 7B (left) shows, the peak phase distribution of the 10 clock gene transcriptions occurred at different ZT times. In the NC-fed group, seven clock genes distributed during light cycle (ZT0–12) had three peak phases: ZT1:30 (Arntl, Clock, and Npas2), 10:30 (Nr1d1), and 12:00 (Nr1d2, Cry2, and Per3). Three clock genes distributed during the dark cycle (ZT12–24) had the two following peak phases: ZT15:00 (Per2) and 21:00 (Rorc and Cry1). However, in response to 10 days of HFI feeding, six clock genes distributed during the light cycle (ZT0–12) had the four following peak phases: ZT1:30 (Arntl), ZT3:00 (Npas2 and Clock), 9:00 (Nr1d1), and 12:00 (Nr1d2 and Per3) (Fig. 7B, right). Four clock genes distributed during the dark cycle (ZT12–24) had the four following peak phases: ZT13:30 (Cry2), 15:00 (Per2), 19:30 (Rorc), and 21:00 (Cry1) (Fig. 7B, right). These data thus suggest that although the core clock machinery remains robustly oscillatory in both groups of mice, HFI feeding still induces some alterations in the amplitude and peak phase of the clock gene transcriptional oscillation in the ELGs. 

Discussion

We found that short-term HFI feeding induced the reprogramming of rhythmic genes in the transcriptomes of the lacrimal glands. This conclusion is mainly based on the following observations: HFI abrogated preexisting rhythmic genes but induced a large number of new rhythmic transcripts; HFI significantly increased the transcriptional oscillation deriving from not-rhythmic genes of the normal-feeding mice; and HFI altered the distribution of biological function pathways of the rhythm genes over a diurnal cycle. These findings indicate that the circadian system of ELGs is sensitive to nutritive challenge. 

As in other peripheral tissues, we found that approximately 3% of the transcripts in the transcriptome of the mouse ELGs are rhythmic. These rhythmic genes mostly encode for basic cellular functions, such as cellular metabolism, circadian rhythm, metabolic degradation, protein processing, and some signaling pathways.³² Similar to other nocturnal animals, the functions and pathways in the murine ELGs depend on the time of day. The dark cycle is the active phase for most physiological functions and activities. 

Based on the transcriptional-level analysis, we found that short-term HFI feeding remarkably remodels the ELG clock. Strikingly, HFI feeding abrogated preexisting rhythmic genes in the NC-fed group but can also generate new rhythmic genes where they do not normally exist in the NC-fed ELGs. These transcripts are mainly derived from preexisting not-rhythmic genes under normal dietary conditions. To understand the biological functions and pathways of these induced cycling transcripts, we analyzed the rhythmic genes specifically expressed by the NC-fed group only, the HFI-fed group only, and rhythmic genes shared by the two groups. The results showed that HFI feeding induces different KEGG pathways, participating in different metabolic pathways and exerting different physiological functions. All of these results suggest that the nutritive perturbation induces reprogramming of the transcriptional system in the lacrimal glands. These findings are consistent with the reprogramming of the mouse liver^34,35 and the suprachiasmatic nucleus induced by high-fat-diet feeding²² and pathologic changes of the rat salivary glands induced by high-sucrose diet.³⁶ 

We found that murine ELGs also rhythmically express most of the core components of the molecular clock system. This suggests that ELGs harbor a circadian clock system and its physiological activities are subject to the influence of an LD cycle similar to other peripheral organs.³² In addition to light cycle misalignment, nutrient stresses such as high-calorie Western diets also interfere with the rhythm of core clock gene transcriptions to a lesser extent.^37,38 Similarly, our data revealed that HFI feeding causes fewer alterations in the clock gene transcription of the lacrimal glands—attenuations, enhancements, or phase shifts—although the transcription is still robust. These data support the concept that circadian oscillation within the core clock genes is strongly resistant to perturbation, whereas clock output genes are more sensitive to food as a zeitgeber^33,39; however, it is unclear whether minor changes in these core clock genes are responsible for reprogramming other circadian output genes. Once elucidated, it is possible to reverse the negative effects of HFI feeding on the lacrimal glands by targeting the expression and function of core clock genes. 

Although our investigation provides a comprehensive understanding of the daily transcriptional rhythm of the murine ELGs, our data have certain limitations in explaining some of the functions of the human lacrimal glands. First, humans and rodents have different lacrimal gland composition and anatomical locations.³⁰ Humans have only one main lacrimal gland on the dorsal and lateral sides of the eyeball and small accessory glands on the conjunctiva. Although rodents also have a lacrimal gland system, it consists of two parts: the ELGs and the intraorbital Harderian gland.⁴⁰ Second, mammals have different circadian timing systems. Given that mice are nocturnal and humans are diurnal,⁴¹ transcriptional outputs of the human lacrimal glands over an LD cycle are hypothesized to be the opposite of that of the murine lacrimal glands. Finally, we observed only the ELG rhythmic gene transcriptome brought about by short-term HFI-feeding metabolic stress. Thus, the results do not reflect certain effects from long-term HFI feeding or from other sugars or carbohydrates on the ELGs. 

In summary, our results demonstrate that short-term HFI feeding, even for just for a few days, profoundly reprograms the circadian transcription of lacrimal glands. This finding suggests that the circadian transcriptome of the lacrimal glands has an intrinsic plasticity that is influenced by the stimulation of exogenous nutrients. This also may mean that metabolic stress caused by HFI can alter the normal physiological function of the lacrimal glands. Dysfunctional lacrimal glands may increase susceptibility to dry eye. Limited clinical evidence shows that patients with metabolic syndrome and type 2 diabetes mellitus present with lower tear volumes and a higher incidence of lacrimal gland hypofunction than age-matched controls.^42,43 Further analysis of the reprograming mechanism may provide insight into the onset and progression of disorders associated with the secretory function of the lacrimal glands. 

Acknowledgments

Supported by grants from the National Natural Science Foundation (nos. 81470603, 81770962, and 81700808) and NIH grant (5R01EY018239-08). 

Disclosure: D. Lu, None; C. Lin, None; X. Jiao, None; Z. Song, None; L. Wang, None; J. Gu, None; Z. Li, None