Integrated DNA Methylation and Transcriptomics Analyses of Lacrimal Glands Identify the Potential Genes Implicated in the Development of Sjögren’s Syndrome-Related Dry Eye

Introduction

Sjögren’s syndrome (SS) is a chronic autoimmune disease characterized by lymphocytic infiltration of salivary and lacrimal glands (LGs) leading to a clinical picture of dry mouth and dry eye.¹ SS-related dry eye is a kind of refractory dry eye disease accompanied by secretory dysfunction of LGs, which may result in severe ocular surface damage and visual impairment.² At present, the etiology of SS-related dry eye is incompletely understood, which limits the development of diagnosis, management, and effective therapies for SS-related dry eye. Although genetic and environmental factors are thought to contribute to the development of SS-related dry eye,³ increasing evidence has demonstrated that epigenetic modifications, such as DNA methylation, histone modification, non-coding RNA, and m6A methylation, play a key role in the development of various autoimmune eye diseases, including SS-related dry eye.^1,4–6

DNA methylation is an important epigenetic modification that can block transcription by recruiting proteins involved in gene repression or by interfering with the binding of sequence-specific transcription factors to their target sequences in candidate genes and lead to changes in gene expression.⁷ It has been reported that abnormal DNA methylation is critically implicated in the development and progression of autoimmune diseases, including SS.^8–10 For example, Yu et al observed that CD4⁺ T cells of primary SS patients had a significantly hypermethylated FOXP3 promoter compared to healthy controls, and the mRNA and protein levels of FOXP3 were negatively correlated with the mean methylation level in the promoter region.¹¹ Besides, Konsta et al found decreased global DNA methylation in minor salivary glands from SS patients compared with those in non-SS patients. They further demonstrated that the reduction of global DNA methylation was associated with lymphocyte infiltration and the levels of circulating anti-SSB/La Ab, but not with disease activity.¹² In a recent paper published in 2023, Chi et al identified two subgroups of SS patients with significantly different DNA methylation levels at the MHC region, based on the genome-wide DNA methylation data from labial salivary gland tissue of SS patients.¹³ Although the change of DNA methylation in peripheral blood cells and salivary glands of SS patients has been widely studied, little is known about the DNA methylation profile in LGs as well as the correlation between DNA methylation and disease progression of SS-related dry eye.

The non-obese diabetic (NOD) mice, which spontaneously develop ocular surface dryness and autoimmune dacryoadenitis that resemble human SS-related dry eye, have been widely used in SS dry eye-related studies.¹⁴ For example, Shin S et al demonstrated that treatment with mesenchymal stem cells could ameliorate the clinical manifestations of SS-related dry eye in NOD mice by increasing tear secretion, decreasing corneal stain scores, recovering goblet cell counts in the conjunctiva, and attenuating foci infiltrations in the LGs.¹⁵ In our previous study, we found that Fenofibrate could alleviate SS–like dacryoadenitis in NOD mice via modulating Th1/Th17 and Treg cell responses.¹⁶ Furthermore, it has been reported that the NOD mice have significant genetic homology with humans, which is helpful in studying the pathogenic mechanism of SS-related dry eye.¹⁷ Given that male NOD mice primarily developed autoimmune dacryoadenitis¹⁸ while female NOD mice developed autoimmune sialadenitis first and were easily affected by diabetes,^19,20 we here used male NOD mice for the exploration of SS-related dry eye etiology.

In this paper, we investigated the DNA methylation and transcriptomics alterations in the LGs of male NOD mice from 4 weeks old to 16 weeks old by using RRBS and RNA-seq, and analyzed the correlation between methylation and gene regulation to identify the differential methylation-regulated genes correlated to the pathogenesis of SS-related dry eye.

Materials and Methods

Animals and Management

Male non-obese diabetic mice (NOD/ShiLtJ) were obtained from Gem Pharmatech Co., Ltd (Nanjing, China). Age- and gender-matched normal BALB/c mice, also from Gem Pharmatech Co., Ltd, served as controls.^21,22 All mice were bred under specific pathogen-free conditions at the animal center of Tianjin Medical University Eye Hospital. All animal experiments were conducted in accordance with the guidelines of the ARVO Statement on Use of Animals in Ophthalmic and Vision Research and approved by the Laboratory Animal Care and Use Committee of Tianjin Medical University.

The NOD mice were randomly assigned into different groups (4-, 8-, 12- and 16-week-old groups) and those accompanied by any pre-existing eye defects were excluded before the experimentation. Mice in each group were maintained for 4, 8, 12 and 16 weeks with food and water available and ad-libitum, respectively. The number of mice used for RRBS and RNA-seq was 2 in each group. And we used 6–9 NOD mice for each group to investigate the development of SS-related dry eye.

Tear Secretion Test

Tear volume was quantified by phenol red cotton threads (Jingming Co., Ltd, Tianjin, China). Briefly, the threads were held with forceps and inserted into the lower conjunctival fornix close to the lateral canthus for 30 seconds. The levels of basal tear production were expressed as millimeters of the length of thread wetted and turned red by tears. Each eye of the mice was measured and tested three times. The average value was determined as the final length.

Corneal Fluorescein Staining

1 µL of 2% sodium fluorescein was instilled into the lateral conjunctival sac of each eye of mice followed by 3 artificial blinks. After 2 minutes, fluorescein staining of corneal epithelial defects was visualized and photographed by a Phoenix MICRON IV Anterior Segment Slit Lamp system (Phoenix MICRON, New York, USA). The staining of cornea epithelial was scored according to the grading system as described previously.²³

Histological Assessment of LGs

For histopathological evaluation of LGs, tissues were fixated in 10% formalin, embedded in paraffin, and sectioned as described previously.¹⁶ The LG sections were stained routinely with hematoxylin and eosin (H&E) and then photographed by light microscopy (BX51; Olympus Corporation, Tokyo, Japan). Inflammation was quantified by standard focus-scoring with a focus defined as an aggregate of at least 50 lymphocytes and the focus score defined as the number of foci per 4 mm² of tissue. Tissue areas were calculated by ImageJ software (National Institutes of Health, Bethesda, MD, USA).^24,25

Immunofluorescence Staining

5 μm LG tissue paraffin slides were incubated at 70°C for 1 h and deparaffinized in xylene. After rehydration in a descending series of ethanol, the slides were subjected to citric acid antigen retrieval. Following blocking with 5% goat serum, LG sections were incubated with anti-CD4-rabbit IgG (1:300; Abcam, Boston, MA, USA) at 4°C overnight. Subsequently, tissue sections were incubated with Alexa Fluor 488–conjugated goat antirabbit IgG (1:500; ImmunoWay, Plano, TX, USA) for 2 h at room temperature being protected from light. DAPI (Solarbio) was used to label nuclear DNA. And images were captured using the laser scanning confocal microscope (LSM800, Zeiss, Germany).

Reduced Representation Bisulfite Sequencing (RRBS)

LG tissues were carefully dissected from NOD mice of 4, 8, 12 and 16 weeks of age, respectively. Left LG tissues were subjected to histopathological assessment, and the right LG were immediately immersed into liquid nitrogen, and stored at −80°C for subsequent investigation. Genomic DNA from pooled LG samples was extracted by TIANamp Genomic DNA Kit (TIANGEN) following the manufacturer’s protocols. A total amount of 1 µg genomic DNA spiked with unmethylated lambda DNA was firstly digested by MspI enzyme followed by libraries construction as the Illumina Pair-End protocol and reduced representation bisulfite sequencing (RRBS) described previously.²⁶ Briefly, the purified and fragmented DNA was processed to repair, blunt and phosphorylate ends. After adenylation, cytosine-methylated adaptors were ligated to the fragmented DNA according to the manufacturer’s instructions. The DNA fragments were then purified using MinElute PCR Purification Kit (Qiagen) and treated with ZYMO EZ DNA Methylation-Gold Kit™ to convert unmethylated cytosine into uracil before the PCR amplification was carried out. The insert size of the selected library was analyzed on the Bioanalyzer analysis system (Agilent, Santa Clara, USA), and library concentration was quantified by the real-time quantitative PCR (qRT-PCR), subsequently sequenced on Illumina sequencing platform.

RNA Isolation and RNA-Sequencing

The RNA-sequencing (RNA‑seq) libraries were generated by E-GENE Co., Ltd (Shenzhen, China) following the manufacturer’s protocols. Briefly, the total RNA was extracted from LG tissues using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocols. After treatment with RNase-free DNase I, 1 µg of total RNA was firstly purified into mRNA containing poly(A) using Oligo(dT) beads. The captured mRNA was then fragmented, reverse transcribed, converted to double-stranded cDNA, and adapter-ligated with unique indices. Finally, PCR was carried out to amplify the adapter-ligated to cDNA and the libraries were validated using Agilent Bioanalyzer 2100 for quality control, followed by qRT-PCR quantification and sequencing on the Illumina sequencing platform.

Data Analysis and Functional Enrichment Analysis

Adapter sequences were trimmed by using cutadapt. And the cleaned reads were mapped back to the genome using BSMAP software version 2.90.²⁷ Only unique mapped reads were used to calculate methylation ratios. Methylation ratios were extracted from BSMAP output (SAM) using a Python script (methratio.py) which is distributed with the BSMAP package. We choose cytosines in a CpG context with sufficient sequencing depth (greater than or equal to 5x coverage) for further analysis. Differentially methylated cytosines (DMCs) and differentially methylated regions (DMRs) were separately detected using MethylKit²⁸ and metilene²⁹ in de-novo mode among CpG sites with at least 5x coverage. Then detected DMRs were filtered according to the standard: (1) P-value must be less than 0.05; (2) methylation level difference must be greater than 0.1; (3) CpG number contained in DMR must be greater than 5; (4) the length of the DMR must be greater than 50bp.

Quality control of raw reads was performed using FastQC (version 0.11.5) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and clean reads were obtained following the removal of low‑quality reads with Trimmomatic (version 0.38) software.³⁰ We then use Hisat2 (v2.2.1)³¹ to perform clean reads mapping to Mus musculus genome reference (mm10). Gene expression profiling was based on the number of mapped reads. Transcripts Per Kilobase of exon model per Million mapped reads values were calculated by StringTie (v2.1.7)³² and were used to estimate the expressed values of genes. The differentially expressed genes (DEGs) among NOD mice groups were obtained using edgeR (v1.34.0)³³ with a Padj-value cutoff < 0.05 and an absolute log₂(fold-change) of > 1.

Principle component analysis (PCA) was performed to visualize the statistical robustness of expression data by utilizing the R language models (2.18.1). Data were clustered based on hierarchical clustering analysis (HC) using R language factoextra (1.0.7) and each row was normalized using a z-score. The circle plot was based on R language (version 4.1.1). Functional annotations were subsequently conducted to elucidate significant biological implications of differentially methylated-related genes (DMGs) and DEGs in different NOD mice groups according to Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genome (KEGG) pathway analyses, provided by Cluster Profiler (v4.4.2)³⁴ or DAVID Bioinformatics Resources (version 2021).³⁵ Padj-value < 0.05 was defined as significant in all of the enrichment analysis.

Correlation Analysis

DNA methylation is known to have a regulatory effect on gene expression. The DMR annotation results of CGs in the promoter region and gene body region were combined with the DEGs results for VENN analysis and single-sample gene-wide association analysis. Additionally, enrichment analysis of genes related to DMRs was implemented by DAVID Bioinformatics Resources (version 2021) to assess the potential biological functions. GO terms and KEGG pathways with Padj-value <0.05 were considered significantly enriched.

Quantitative Real-Time PCR (qRT-PCR)

The total RNA of LGs was isolated using EZ-press RNA Purification Kit (EZBioscience, Roseville, MN, USA) and reverse-transcribed to synthesize cDNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocols. qRT-PCR was performed using gene-specific primers and SYBR Green PCR Master Mix (Roche, Mannheim, Germany) on Roche LightCycler 480 II Analyzer. Gene expression relative to housekeeping gene GAPDH was calculated using the comparative threshold cycle (2-ΔΔCt) method. The sequences of qRT-PCR primers are listed in Table 1.

Table 1 mRNA Primer Sequences

Methylation-Specific PCR (MSP)

The methylation levels in the promotor region of Itgal, Icosl, Vav1 and Irf4 were further examined through MSP. The genomic DNA was extracted from the LGs of NOD mice using the Animal Genomic DNA Quick Extraction Kit (D0065 M, Beyotime, China) and further modified with sodium bisulfite using the EZ DNA Methylation-Gold Kit (Zymo, Foster City, CA) following the manufacturer’s protocols. Primers for methylated and unmethylated promoter sequences were designed using MethPrimer software and listed in Table 2. Subsequently, 2% agarose gels were utilized to analyze the PCR product.

Table 2 Sequences of Primers Used for MSP

Statistical Analysis

All statistical analyses of the data were performed with GraphPad Prism software (GraphPad Software, San Diego, CA, USA). The normality of continuous variable was evaluated by Shapiro–Wilk test. All data are expressed as mean ± SD. For normally distributed data, differences among groups were analyzed by using one-way ANOVA. Kruskal–Wallis test was used for nonparametric analysis according to the number of groups. The association analysis was evaluated by the Spearman rank correlation test. P-value < 0.05 was considered statistically significant.

Results

The Severity of Dry Eye Showed Progressive Exacerbation with Increasing Age of NOD Mice

To determine the development of SS-related dry eye in animal models, NOD mice from 4 to 16 weeks of age were selected in our study for the evaluation of corneal fluorescein staining, tear secretion and LGs histology^36,37 (Figure 1A). As shown in Figure 1B and C, corneal epithelial defects were detected in NOD mice as early as 8 weeks of age and the corneal staining score showed a significant increase from 12 weeks old compared to those of 4-week-old NOD mice. Tear secretion volume measured by phenol red thread exhibited significant reduction from 12 weeks of age and decreased continually in NOD mice (Figure 1D). In contrast, the corneal fluorescein staining and tear secretion did not show significant alterations in BALB/c mice from 4 to 16 weeks old (Figure 1B-D). Consistent with the clinical results, histological analysis of LGs revealed that lymphocytic infiltration was easily detected as early as 8 weeks old in NOD mice and increased significantly until 16 weeks old, while there was no significant lymphocytic infiltration in the LGs of 4 to 16-week-old BALB/c mice (Figure 1E and F). Furthermore, the correlation analysis revealed that the extent of inflammatory cell infiltration in LGs was positively correlated with the corneal fluorescein score and negatively correlated with tear secretion volume (Figure 1G and H). Altogether, these data suggested that the clinical manifestations and histopathological alterations of SS-related dry eye in NOD mice were observed from 8 weeks and gradually worsened until 16 weeks of age.

Figure 1 The severity of dry eye showed progressive exacerbation with increasing age of NOD mice. (A) Time-point selection and experimental design flow for the course studies of SS-related dry eye in NOD mice. n=6-9 mice per group. (B and C) The representative photographs (B) and grading scores of corneal fluorescein staining (C) in different age groups of NOD mice and healthy BALB/c controls. n=6 mice per group. (D) Phenol red cotton threads measurement of tear secretion in NOD mice and healthy BALB/c mice at different time points. n=6-9 mice per group. (E and F) Histopathological analysis of the LGs in NOD mice at different ages. Representative images (E) and quantification (F) show lymphocytic infiltration in LGs. n=3. (G) Spearman rank correlation between the lymphocytic infiltration in the LGs and corneal fluorescein staining scores. n=3. (H) Correlation analysis between the LGs inflammation and tear secretion volume in NOD mice. n=3. Values are represented by means ± SD. Kruskal–Wallis test for multiple comparisons was used. *p<0.05, **p<0.01, or ***p<0.001, compared with 4-week-old group.

DNA Methylation Levels Were Increased in the LGsof NOD Mice During the Development of SS-Related Dry Eye

Recently, increasing evidence has demonstrated the contributing role of DNA methylation in the initiation and progression of autoimmune diseases.³⁸ To investigate whether DNA methylation is involved in the development of SS-related dry eye in animal models, we isolated LGs of 4 to 16-week-old NOD mice and profiled DNA methylation across the genome using the RRBS technique. The bisulfite conversion rates of each sample were higher than 98% (Table 3) and the coverage rate of C sites was uniform for each group (Figure 2A and B), which indicated integrity, accuracy, and high uniformity of the sequencing data. Principal component analysis (PCA) and hierarchical clustering (HC) showed distinct DNA methylation profiles in LGs from different ages of NOD mice (Figure 2C and D). By comparing the global methylation rates of all C sites among the four groups, we found that DNA methylation mainly occurred on CG cytosines (mean methylation level, m= 62.8%) and the methylated rates of CG cytosines in mice at 8, 12, and 16 weeks old were higher in comparison to 4-week-old group (Figure 2E-G). Further, we observed that the differential CG sites were mainly located in regions 2-kb upstream of the gene body (promoter region) (Figure 2H). To define the specific position of methylated sites in the promoter, the promoter was further divided into five gene functional regions: CGI-upshelf, CGI-upshore, CG islands (CGI), CGI-downshore, and CGI-downshelf. As depicted in Figure 2I and J, the 8-week-old group showed significantly higher CG methylation levels in the CGI region than the other three groups. These data collectively suggested that the methylation levels of LGs were increased during the development of SS-related dry eye.

Table 3 Basic sequencing quality control, bisulfite conversion rates, and map rate of RRBS libraries

Figure 2 DNA methylation levels were increased in LGs of NOD mice during the development of SS-related dry eye. (A) DNA methylation rate under the contexts of CG, CHG, and CHH in NOD mice at different stages of SS-related dry eye (4, 8, 12 and 16 weeks old) by RRBS. (B) The genome coverage rate of C sites in LGs of different NOD groups. (C) Principal component analysis (PCA) based on individual cytosines. (D) Hierarchical clustering (HC) analysis based on methylation levels at CG sites. (E) Global DNA methylation levels of all C sites across the genome in LGs of NOD at different stages of SS-related dry eye. (F) Distribution of methylation levels at CG sites in different groups. (G) Methylation rates of CG sites at different disease stages in NOD. (H) DNA methylation levels along the gene body, 2-kb upstream of the TSS and 2-kb downstream of the TES for all genes. (I) DNA methylation levels of promoter CG island and surrounding areas. (J) Heatmap of methylation levels at CG sites in the functional region of the gene. Red represents hypermethylated; blue represents hypomethylated.

Functional Enrichment of the Differentially Methylated Regions-Related Genes (DMGs) Identified in LGs of SS-Related Dry Eye Mice

Next, we attempted to explore the association between methylation differences and the pathogenesis of SS-related dry eye by calculating the differentially methylated regions (DMRs) among 4-, 8-, 12-, and 16-week-old NOD mice under CG contexts. Applying the number of differentially methylated CG sites in the region ≥ 5 and the length of the DMR must be greater than 50bp (P-value < 0.05) returned 4132 DMRs (including 2956 hyper- and 1176 hypomethylated), 23,396 DMRs (17,927 hyper- and 5469 hypomethylated) and 26,135 DMRs (20,422 hyper- and 5713 hypomethylated) in 8-, 12-, and 16-week-old groups compared with 4-week-old control group, respectively (Figure 3A). Consistent with this, the circle plot showed that there were more hyper-DMRs than hypo-DMRs along the chromosome. Most DMRs were located on autosomes, while fewer DMRs were located on sex chromosomes (Figure 3B-D). Also, more DMRs were located in exons than promotor regions among different aged mice (Figure 3E-G). To define the function of DMRs-related genes (DMGs), we performed GO enrichment and KEGG pathway analysis. The GO analysis revealed that DMGs were significantly enriched in certain molecular function, including GTPase activator activity, nucleoside-triphosphatase regulator activity and GTPase regulator activity (Figure 3H), which have been considered vital for regulating immune cells migration in SS.^39–41 Furthermore, the DMGs were enriched in biological process involved in synapse organization, axonogenesis and cell junction assembly (Figure 3H). PI3K-Akt signaling pathway, MAPK signaling pathway and Ras signaling pathway were also found to be enriched through KEGG analysis (Figure 3I). Each of these pathways has previously been linked to the development of dry eye diseases.⁴²

Figure 3 Functional enrichment of the differentially methylated regions-related genes (DMGs) identified in LGs of SS-related dry eye mice. (A) The numbers of DMRs in 8-, 12- and 16-week-old NODs versus 4-week-old controls under CG contexts. (B-D) Comparison of DMRs on different chromosomes between groups at different stages of SS-related dry eye. (E-G) The bar graphs showed the genomic location (exons and promoters) of identified DMRs among different aged NOD groups. (H) GO enrichment analysis for DMGs. The color of the dot indicates the P-value and the size of the dot represents the enrichment score. (I) Top 10 KEGG pathways analysis of DMGs. (J-L) The changes of Dnmt1 (J), Dnmt3a (K) and Dnmt3b (L) in the LGs were determined by qRT-PCR during SS-related dry eye progression. n=3 mice per group. Data are expressed as the mean ± SD. One-way ANOVA followed by the Tukey’s test for multiple comparisons was used. *p < 0.05, **p < 0.01, compared with 4-week-old group.

DNA methylation is predominantly catalyzed through a family of DNA methyltransferases (DMNTs), which consists of three members: DNMT1, DNMT3A and DNMT3B.^43–45 Thus, we investigated the mRNA level of three DNMT subtypes (Dnmt1, Dnmt3a, and Dnmt3b) in LGs of different age groups of NOD mice. The results showed that the mRNA level of Dnmt1 and Dnmt3b in LGs was significantly elevated in 8-week-old group compared with 4-week-old group, while there was no significant difference in Dnmt3a expression (Figure 3J-L), suggesting that Dnmt1 and Dnmt3b may play a major role in DNA methylation regulation of LGs.

Transcriptome Analysis of LGs During the Development of SS-Related Dry Eye

As the DNA methylome has been altered during SS-related dry eye progression, it is necessary to study whether the transcriptional profiles are affected. To delineate possible dysregulation of gene transcription during the development of SS-related dry eye, RNA sequencing (RNA-seq) was conducted on LGs from the 4-, 8-, 12-, and 16-week-old NOD mice. High-quality data, characterized by more than 70% mapping rates of clean reads for all samples, were used for further analysis (Figure 4A and Table 4). The PCA and HC results showed a clear separation among the four groups (Figure 4B-C). Transcriptome analysis revealed 1321, 2549 and 3712 differentially expressed genes (DEGs) in the 8-, 12- and 16-week-old groups compared with 4-week-old NOD mice, respectively (Figure 4D). By performing Venn diagram analysis, we further identified 858 overlapping upregulated DEGs in different age groups (Figure 4E). These DEGs were mainly involved in immunological synapse, positive regulation of T cell activation, and MHC class II protein complex binding (Figure 4F-H). Besides, KEGG analysis revealed that the 858 genes were significantly enriched in cell adhesion molecules, Th17 cell differentiation, antigen processing and presentation, and T/B cell receptor signaling pathway (Figure 4I). These results indicated a wide range of immune response-related genes are upregulated in NOD mice during the development of SS-related dry eye.

Table 4 Basic Sequencing Quality Control and Map Rate of RNA-Seq Libraries

Figure 4 Transcriptome analysis of LGs during the development of SS-related dry eye. (A) The mapping rates of clean reads for all samples. (B and C) PCA analysis and HC based on gene expression in different groups. (D) Heatmap of differentially expressed genes (DEGs) in each group. Threshold log2(fold-change)>1 or log2(fold-change) <-1, p-value<0.05. The top of the graph represents the sample name, the rows represent the genes, the expression levels are indicated by colors (red indicates high gene expression level, blue indicates low gene expression level), the left side indicates the gene clustering, and the top side indicates the clustering among the samples. (E) VENN plots of upregulated DEGs in the three groups. (F-H) GO enrichment analysis of overlapping upregulated DEGs. GO analysis of predicted genes was performed according to cellular component (CC) (F), biological process (BP) (G), and molecular function (MF) (H). The left side of the graph represents GO categories, the color corresponds to the -log10 of the p-value, and the size of dot represents the number of enriched genes. (I) KEGG pathway analysis of DEGs. The top 20 pathways were summarized.

Correlation Analysis between DMRs and genome-Wide Gene Expression

To identify the potential relationship between DNA methylation changes and gene expression in LGs of NOD mice at different disease stages, we performed an integrated analysis of DNA methylation and transcriptomic profiles (Figure 5A). Given that the critical role of DNA methylation within promoter regions in gene transcriptional regulation,⁴⁶ we focused on DMRs within gene promoter regions for further analysis. As shown in Figure 5B, a total of 140 overexpressed genes with hypomethylated DMRs in the promoters were obtained by the Venn diagram in 8-,12- and 16-week-old NOD mice (2, 46, and 92 genes, respectively) compared to the control group. The results of GO analysis showed that these upregulated genes with hypomethylated DMRs were significantly enriched in positive regulation of T cell proliferation, T cell activation, and T-helper 17 cell lineage commitment (Figure 5C). And the top 5 pathways enriched by KEGG were intestinal immune network for IgA production, T/B cell receptor signaling pathway, Th17 cell differentiation signaling pathway and cell adhesion molecules (Figure 5D).

Figure 5 Correlation analysis between DMRs and genome-wide gene expression. (A) Study design diagram for integrated analysis of methylome and transcriptome. (B) VENN analysis of upregulated DEGs with hypomethylated DMRs. (C and D) GO and KEGG analysis of genes upregulated by promoter hypomethylation in different NOD groups.

Identification of Key DNA Methylation-Regulated Genes Associated with SS-Related Dry Eye

Based on the literature, GO, KEGG and PPI network analysis, eight upregulated genes with hypomethylation in the promoter region, including Icosl, Batf, Irf4, Itgal, Vav1, Tgfβ1, Jak2 and CD86, were screened (Figures 6A and B, S1A and Tables S1-S3) and validated by qRT-PCR. The data showed that the expression levels of Icosl, Irf4, Vav1 and Jak2 were significantly increased in LGs of 12-week-old NOD mice, and Itgal was significantly upregulated in LGs of 16-week-old mice compared to the control group (Figures 6C and S1B). On the contrary, Batf, CD86 and Tgfβ1 did not have significant alteration among the four groups (Figure S1B). Further correlation analysis revealed that the expression levels of Itgal, Vav1, Irf4, and Icosl were positively correlated with the extent of inflammatory cell infiltration in LGs (Figure 6D). To clarify the relationship between these four genes (Itgal, Vav1, Irf4, and Icosl) and T cell response in this disease context, we performed immunofluorescence assay of LGs and observed that CD4⁺ cell infiltration dramatically increased in LGs of SS-related dry eye mice compared with the control mice (Figure 6E). Also, correlation analysis revealed that the expression levels of Itgal, Vav1, Irf4, and Icosl showed significantly positive correlation with the extent of CD4⁺ cell infiltration in LGs (Figure 6F). In addition, using methylation-specific PCR (MSP), we observed hypomethylation of the Irf4 and Itgal promoters in NOD mice with SS-related dry eye compared to those 4-week-old mice without dry eye (Figure 6G). Therefore, these data indicate that Itgal, Vav1, Irf4, and Icosl are closely associated with the development of SS-related dry eye by modulating T cell response, and DNA hypomethylation of the promotor region may be responsible for the upregulation of Itgal and Irf4 in SS-related dry eye.

Figure 6 Identification of key DNA methylation-regulated genes associated with SS-related dry eye. (A) Flow chart of key gene screening. (B) The heatmap of methylation levels in promoters and expression levels of candidate functional genes related to SS-related dry eye. (C) Validation of the mRNA expression levels of candidate genes in LGs by qRT-PCR. n = 3 mice per group. One-way ANOVA and Tukey’s test were used for multiple comparisons. Values are represented by means ± SD. *p<0.05, **p<0.01, compared with 4-week-old group. (D) Correlation analysis between the expression levels of candidate genes and the lymphocyte infiltration in LGs. Spearman rank correlation tests were used. (E) Immunofluorescence labeling of DAPI and CD4 in LGs of the 4-, 8-, 12-, and 16-week-old groups. CD4 (upper panel, green), DAPI (middle panel, blue), and the merge (lower panel) are shown. Arrows show CD4+ cell infiltration lesions. **p<0.01, or ***p<0.001, compared with 4-week-old group. Scale bars = 20 μm. (F) Correlation analysis between the expression levels of candidate genes and the extent of CD4+ cell infiltration in LGs. Spearman rank correlation tests were used. (G) Representation of MSP for Icosl (n = 3), Irf4 (n = 3), Vav1 (n = 3) and Itgal (n = 3) in NOD mice with SS-related dry eye compared to 4-week-old mice without dry eye. Bands in lanes labeled “M” and “U” are PCR products amplified with methylation- and unmethylation-specific primers.

Discussion

Accumulating evidence has demonstrated that epigenetic factors are critically implicated in the pathogenesis of SS-related dry eye.⁴⁷ DNA methylation, the main focus of epigenetic modification, has been proven to be altered in peripheral blood cells and salivary glands of SS patients and closely associated with the development of SS.⁴⁸ Nevertheless, the DNA methylation profiling in LGs and their roles in SS-related dry eye remains unclear. Here, by performing RRBS, we found that the methylation levels were increased in LGs of 8-, 12- and 16-week-old NOD mice when compared with the 4-week-old mice. Further combined analysis of DNA methylome and transcriptome revealed that several genes related to T cell immunity, including Itgal, Icosl, Vav1 and Irf4, were upregulated during the development of SS-related dry eye. Also, hypomethylation of the Irf4 and Itgal promoters were observed in SS-related dry eye mice. Our results thus suggest that altered DNA methylation is closely linked to the occurrence and development of SS-related dry eye, possibly in part because it affects the expression of T cell-related genes in LGs.

It has been reported that CG methylations play a vital role in regulating gene transcription and are critically involved in various autoimmune diseases.^49,50 In this study, we found the methylation levels of CG sites were significantly elevated in LGs of NOD mice during the development of SS-related dry eye, suggesting that hypermethylation of genomic DNA in LGs may be closely associated with the pathogenesis of SS-related dry eye. Interestingly, the clinical manifestations of dry eye and inflammatory cell infiltration of LGs in 8-week-old mice were less severe than those in 12- and 16-week-old mice, but the methylome of LGs has changed on a large scale. This might be due to epigenetic changes taking place earlier than pathophysiological changes.^51–53 Supporting this, DNA methyltransferases, including Dnmt1 and Dnmt3b, were observed to be upregulated in 8-week-old NOD mice. The causes underlying the alterations in DNA methylation especially DNA methyltransferases await further study. In addition, several previous studies have established the association between DNA methylation changes in salivary glands and peripheral blood cells and SS, revealing that genomic DNA methylation levels were decreased in salivary glands⁵⁴ and peripheral blood cells of SS patients.⁵⁵ Differences in species and tissue types may account for the discrepancy in these studies and our findings.

Rab27a, a member of the superfamily of Ras-related small GTPases, has been demonstrated to promote the secretion of cathepsin S in lacrimal gland acinar cells and be closely linked to SS-related dry eye.^56,57 Here, our analysis revealed that Ras signaling pathway and its related molecules were enriched in both GO and KEGG analysis, which indicated that these DMGs may contribute to the occurrence and development of SS-related dry eye by affecting Ras pathway and its downstream molecules. In-depth research is warranted to investigate the changes of these DMGs in SS-related dry eye and their effect on the development of SS-related dry eye. T cells are essential for the immune pathogenesis of SS-related dry eye.⁵⁸ It has been demonstrated that LGs can express major histocompatibility complex class II (MHC II) or human leukocyte antigen-DR and costimulatory molecules, and thereby act as nonprofessional antigen-presenting cells for T cell activation, promoting the development of SS-related dry eye.⁵⁹ Administration of small molecules inhibitor targeting MHC II could reduce the infiltration of Th1 and Th17 cells in the exocrine glands, and improve SS-like symptoms, such as SS-related dry eye, in NOD mice.⁶⁰ Consistent with these findings, our study proved that 858 upregulated DEGs in LGs of 8-, 12- and 16-week-old groups were closely related to T cell immune response, such as positive regulation of T cell activation and peptide antigen assembly with MHC class II protein complex,^58,61,62 indicating that these upregulated genes may contribute to SS-related dry eye by directly promoting T cell activation or acting on MHC class II molecules in LGs to indirectly drive T cell response.

140 overexpressed genes with hypomethylated DMRs in the promoter based on combined RRBS and RNA-seq data were further studied. These genes were mainly involved in T/B cell receptor signaling pathway and Th17 cell differentiation signaling pathway, which have been identified to play key roles in the pathogenesis of SS-related dry eye.^63–65 Among these candidate genes, it was further confirmed that Icosl, Irf4 and Vav1 were upregulated in LGs of 12-week-old NOD mice, and Itgal was upregulated in LGs of 16-week-old mice compared to the control group, suggesting that upregulated Itgal, Icosl, Vav1 and Irf4 might participate in the development of SS-related dry eye. Indeed, we observed that the expression levels of Itgal, Icosl, Vav1 and Irf4 were positively correlated with the severity of CD4⁺ cell infiltration in LGs. By performing MSP assay, we further found decreased methylation levels in the promotor region of Irf4 and Itgal. This indicated that hypomethylation of promotor region may lead to overexpression of Irf4 and Itgal in LGs of SS-related dry eye mice and contribute to the pathogenesis of SS-related dry eye. Consistent with our observations is a previous study showing that the promoter of ITGAL was hypomethylated in T cells from lupus patients, which led to elevated mRNA expression of ITGAL.^66,67 ITGAL, also known as CD11a or LFA-1, has been confirmed to be critical for the occurrence and development of dry eye,^68,69 and the LFA-1 antagonist has been approved by the FDA in 2016 as a new drug for the treatment of dry eye.^70,71 Targeting LFA-1 in dry eye mouse model could reduce ocular surface epithelial cell damage and increase goblet cell density.⁷² Therefore, elucidating the upstream molecules of ITGAL can provide new ideas for the treatment of dry eye diseases. Besides, IRF4, VAV1, and ICOSL have been identified as key regulators of T cell response.^73–75 Our results provide a clue that Icosl, Irf4, and Vav1 may be new biomarkers for the diagnosis of SS-related dry eye. However, the specific role and mechanism of Icosl, Irf4, and Vav1 in SS-related dry eye await further investigation.

The link between epigenetic factors and autoimmune eye disease has gotten more and more attention. To date, few studies have been performed on the role of DNA methylation in LGs of SS-related dry eye. This work could constitute a baseline for other future studies focusing on the potential influence of DNA methylation on the pathogenesis of SS-related dry eye. However, the present study has several potential limitations. First, considering that the samples were taken from NOD mice, the results in this study need to be explored in SS-related dry eye patients by further research. And second, the biological mechanisms of the DMR-related genes and signaling pathways in the progression of SS-related dry eye should be confirmed by further mechanistic studies. Finally, given that the LG tissue consists of a variety of cells, using the whole LG tissue rather than distinguishing different types of cells for sequencing may affect data mining to some extent. Future application of single-cell techniques might provide more accurate results and further help to reveal the underlying pathogenesis of SS-related dry eye.

Conclusion

We demonstrated that genomic DNA methylation levels were increased in LGs of NOD mice during the development of SS-related dry eye. Itgal, Icosl, Vav1 and Irf4, which were upregulated in LGs of NOD mice, were closely associated with T cell-mediated immune response and may play vital roles in the pathogenesis of SS-related dry eye. Additionally, DNA methylation levels of the Itgal and Irf4 promotor region were reduced, which may be responsible for the upregulation of Itgal and Irf4 in SS-related dry eye. These findings can help to further understanding of SS-related dry eye etiology and provide novel targets for the diagnosis and treatment of SS-related dry eye.