Retinal transcriptome profiling identifies novel candidate genes associated with visual impairment in a mouse model of multiple sclerosis

ABSTRACT

Visual impairment is occasionally observed in multiple sclerosis (MS) and its animal model, experimental autoimmune encephalomyelitis (EAE). Although uveitis and optic neuritis have been reported in MS and EAE, the precise mechanisms underlying the pathogenesis of these visual impairments remain poorly understood. This study aims to identify differentially expressed genes (DEGs) in the retinas of mice with EAE to identify genes that may be implicated in EAE-induced visual impairment. Fourteen adult mice were injected with myelin oligodendrocyte glycoprotein_35–55 to induce the EAE model. Transcriptomes of retinas with EAE were analyzed by RNA-sequencing. Gene expression analysis revealed 347 DEGs in the retinas of mice with EAE: 345 were upregulated, and 2 were downregulated (adjusted p-value < 0.05 and absolute log₂ fold change > 1). Gene ontology (GO) analysis showed that the upregulated genes in the retinas of mice with EAE were primarily related to immune responses, responses to external biotic stimuli, defense responses, and leukocyte-mediated immunity in the GO biological process. The expression of six upregulated hub genes (c1qb, ctss, itgam, itgb2, syk, and tyrobp) from the STRING analysis and the two significantly downregulated DEGs (hapln1 and ndst4) were validated by reverse transcription-quantitative polymerase chain reaction. In addition, gene set enrichment analysis showed that the negatively enriched gene sets in EAE-affected retinas were associated with the neuronal system and phototransduction cascade. This study provides novel molecular evidence for visual impairments in EAE and indicates directions for further research to elucidate the mechanisms of these visual impairments in MS.

KEYWORDS:

• Differentially expressed gene
• experimental autoimmune encephalomyelitis
• immune response
• RNA-sequencing
• phototransduction

Introduction

Multiple sclerosis (MS) is a neurodegenerative disorder that commonly affects young adults and is characterized by the inflammatory demyelination in the central nervous system (CNS) (Reich et al. 2018). Patients with MS may also experience optic neuritis and uveitis during disease progression (Ransohoff et al. 2015), which can lead to visual impairments such as decreased visual acuity (Balcer and Frohman 2010), visual field defects (Nakajima et al. 2010), altered color vision (Villoslada et al. 2012), and blurred vision (Casselman et al. 2021). These symptoms are frequently associated with structural impairments in the visual system (Talman et al. 2010; Rosenkranz et al. 2021) and have been reported in approximately 60% of patients with MS (Salter et al. 2013) and impact their quality of life (Noble et al. 2006).

Experimental autoimmune encephalomyelitis (EAE) is a widely used animal model for studying MS that is characterized by T-cell-mediated inflammation and demyelination in the CNS (Shin et al. 1995). In EAE, lymph node-derived encephalitogenic T cells circulate in the blood vessels and infiltrate the CNS parenchyma, resulting in a disease that resembles certain aspects of MS pathology (Sorensen and Ransohoff 1998). During the progression of EAE, the composition of cerebrospinal fluid (CSF) is altered, and pro- and anti-inflammatory cytokines in the CSF are upregulated (Borjini et al. 2016). Thus, exogenous elements affect the brain and spinal cord and induce dysfunction of these organs (Baxter 2007).

In addition to dysfunction of the primary target organs, visual system dysfunction has also been observed in the course of EAE (Shin et al. 2021). In the chronic stage of EAE, affected mice experience axonal damage of the optic nerve (Horstmann et al. 2013), ganglion cell loss with acute optic neuritis (Shindler et al. 2006), and activation of astrocytes and microglia in the visual system (Larabee et al. 2016; Jin et al. 2019; Tassoni et al. 2019), as seen in patients with MS (Green et al. 2010). Nevertheless, the molecular mechanisms underlying visual impairment in CNS autoimmune diseases, such as MS and EAE, remain unknown. In this study, gene expression profiles were analyzed in the retinas of mice in the chronic phase of EAE to obtain insight into the underlying mechanisms of visual impairment in this animal model.

Materials and methods

Animals and EAE induction procedure

Nine-week-old male C57BL/6J mice aged (n = 14 mice/group) were obtained from Daihan Biolink Co. (Chungbuk, Republic of Korea) for use in this study. The mice were housed in a room with a temperature of 23 ± 2°C, a relative humidity of 50 ± 5%, and a 12-h light/dark cycle. They were provided with standard rodent food (Envigo Teklad, Madison, WI, USA) and water ad libitum. All animal care and procedures were conducted in accordance with the ethical guidelines of the Institutional Care and Use Committee of Chonnam National University (approval number: CNU IACUC-YB-2022-113) and followed internationally agreed standards for laboratory animal use and care, as mandated by the National Institutes of Health (NIH).

To induce EAE, mice (n = 14) were immunized subcutaneously on the hind flank with 1 mg/mL myelin oligodendrocyte glycoprotein peptide 35–55 (MOG_35–55; Koma Biotech, Seoul, Republic of Korea) emulsified in complete Freund’s adjuvant (CFA; Sigma-Aldrich, St. Louis, MO, USA) and supplemented with 5 mg/mL Mycobacterium tuberculosis H37Ra (Difco Laboratories, Franklin Lakes, NJ, USA). On days 0 and 2 post-immunization (DPI), 500 ng of pertussis toxin (List Biological Laboratories, Campbell, CA, USA) was injected intraperitoneally into the mice. Control mice (n = 14) remained non-immunized.

Mice were monitored daily for changes in their weights and clinical signs of EAE, which were graded as follows: grade 0 (G.0), no signs; G.1, floppy tail; G.2, mild paraparesis; G.3, severe paraparesis; G.4, tetraparesis; and G.5, moribund or death (Kim et al. 2019). At 28 DPI, retinal tissues were collected from both control and EAE mice to analyze histopathological changes and gene expression profiles. All possible measures were taken to reduce the number of animals used in this study and to alleviate their suffering.

Tissue preparation

After 28 DPI, mice were sacrificed using CO₂ gas inhalation and profound anesthesia. The globes were removed with care. For histopathological examination (n = 3 mice/group), the samples were fixed in Tokuda-Baron fixative, which is composed of 100 ml of glacial acetic acid, 150 ml of 99% ethyl alcohol, 100 ml of 40% neutral buffered formalin, and 650 of ml distilled water, and stored at 4°C overnight (Tokuda et al. 2018). They were then paraffin-embedded and sectioned at a thickness of 5 μm using a microtome (Leica Biosystems, Nussloch, Germany). The sections were stained with hematoxylin–eosin, and immunohistochemistry was performed. For RNA-sequencing (RNA-seq; n = 4 mice/group) and reverse transcription-quantitative polymerase chain reaction (RT-qPCR; n = 7 mice/group) analyses, the retinal tissues from both control and EAE-affected mice were stored at −80°C immediately after sampling.

Immunohistochemistry

Immunohistochemistry was performed using the same protocol used in a previous study (Kim et al. 2019). The primary antibody (#019-19741; Wako Pure Chemical Industries, Ltd., Osaka, Japan) was diluted 1:1,000 against ionized calcium-binding adapter molecule 1 (Iba1), as a marker for microglia and macrophages.

RNA isolation and RNA-seq

The total RNA was extracted using an RNeasy® Mini Kit (#74106, Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The RNA concentration was determined using Quant-IT RiboGreen (#R11490, Invitrogen, Carlsbad, CA, USA). TapeStation RNA ScreenTape (#5067-5576, Agilent, Santa Clara, CA, USA) was utilized to assess the integrity of the total RNA, and only high-quality RNA preparations with an RNA integrity number greater than 7.0 were used for RNA library construction.

Following the protocol outlined in a prior study (Ahn et al. 2022), each sample underwent independently library preparation using 1 µg of total RNA and the Illumina TruSeq Stranded mRNA Sample Prep Kit (#RS-122-2101, Illumina, Inc., San Diego, CA, USA). In the first step, poly-A-containing mRNA molecules were purified using poly-T-attached magnetic beads. The mRNA was then fragmented into small pieces using divalent cations at an elevated temperature. The cleaved RNA fragments were copied into first-strand cDNA using SuperScript II reverse transcriptase (#18064014, Invitrogen) and random primers. Second-strand cDNA synthesis was performed using DNA polymerase I, RNase H, and dUTP. The cDNA fragments underwent an end-repair process, followed by the addition of a single ‘A’ base and subsequent ligation of the adapters. The final cDNA library was then generated by purifying and enriching the products with PCR.

The libraries were quantified utilizing KAPA Library Quantification kits for Illumina sequencing platforms in accordance with the qPCR Quantification Protocol Guide (#KK4854, Kapa Biosystems, Wilmington, MA, USA) and qualified using TapeStation D1000 ScreenTape (#5067-5582, Agilent). Paired-end sequencing reads were generated on the Illumina sequencing NovaSeq platform. Before initiating the analysis, Trimmomatic v0.38 (ILLUMINACLIP: 2:30:10:5 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:15 MINLEN:36) was used to remove adapter sequences and trim bases with poor base quality. The cleaned reads (average 98% reads at mean length 100 bp) were aligned to the Mus musculus (mm10) using HISAT v2.1.0 (Kim et al. 2015) based on the HISAT and Bowtie2 (Langmead and Salzberg 2012) implementations. For this alignment, we obtained the reference genome sequence from the NCBI Genome assembly and gene annotation data from the NCBI RefSeq database. Post-alignment, the generated data in SAM file format were sorted and indexed using SAMtools v1.9 (Danecek et al. 2021). Following alignment, transcripts were assembled and quantified using StringTie v2.1.3b (Pertea et al. 2015; Pertea et al. 2016).

Differentially expressed genes (DEGs), enrichment analysis, and protein–protein interaction (PPI) analysis

Differential expression was analyzed using R version 4.2.2 (R Core Team 2013) and the R package DESeq2 (Love et al. 2014) for RNA-sequencing. Differentially expressed genes were determined using the DESeq2 (Love et al. 2014) Wald test, and p-values were adjusted using the Benjamini–Hochberg procedure for controlling the false discovery rate (FDR) (Benjamini and Hochberg 1995). The genes with an absolute log₂ fold change greater than 1 and an adjusted p-value less than 0.05 were considered to be differentially expressed. To examine the gene ontology (GO) enrichment of the DEGs, we used ShinyGO 0.77 (Ge et al. 2020) (accessed on February 26, 2023) and analyzed three GO datasets including biological processes (BPs), cellular components (CCs), and molecular functions (MFs) to analyze upregulated genes (Hong et al. 2014).

To construct PPI networks for DEGs and to identify essential nodes in the network with high confidence (more than 0.7), we used the Search Tool for the Retrieval of Interacting Genes/proteins (STRING) database (Szklarczyk et al. 2021). We then applied the CytoHubba plug-in (Chin et al. 2014) in Cytoscape software (Shannon et al. 2003) to identify the hub genes from the PPI network based on five ranking algorithms including betweenness, closeness, degree, maximal clique centrality, and maximum neighborhood component. We selected the top 50 nodes with the highest scores of the five ranking algorithms and took the intersection of the outcomes of the five algorithms. Overlapping hub genes were used to select the central hub genes.

We additionally utilized gene set enrichment analysis (GSEA) version 4.2.3 (Subramanian et al. 2005) to analyze the gene sets enriched in the RNA-Seq data. The GSEA user guide was followed to configure the analysis parameters for the preceding comparisons. In this study, GSEA was employed to identify dysregulated signaling pathways of all DEGs in the control versus EAE-affected groups using curated gene sets from the Molecular Signatures Database (MSigDB). An enriched dataset expression matrix was generated using the curated gene sets at the cutoff levels of a false discovery rate (FDR) q-value < 0.05 and nominal (NOM) p-value < 0.05.

RNA extraction, cDNA synthesis, and RT-qPCR

The protocols for RNA extraction, complementary DNA (cDNA) synthesis, and RT-qPCR were taken from a previous study (Weerasinghe-Mudiyanselage et al. 2022). Briefly, total RNA was extracted from the retinas of control and EAE mice with QIAzol Lysis Reagent (#79306; Qiagen) based on the manufacturer’s instructions. cDNA synthesis was performed using the SuperiorScript III cDNA synthesis kit (#EZ405S, Enzynomics, Daejeon, South Korea). The resultant cDNA was diluted to a final concentration of 8 ng/µL with RNase-free water and stored at −80°C. RT-qPCR was conducted using the mic instrument (Bio Molecular Systems, Potts Point, Australia) and 2× Quantity SYBR Green (PhileKorea Co., Ltd., Seoul, Korea). The sequences of RT-qPCR primers are listed in Table 1. The annealing temperature for the reaction was 58°C. Amplification curves were generated, and threshold cycle (Ct) values were calculated using the built-in software. The expression levels of the target genes were normalized to the Gapdh reference gene using the 2^−ΔΔCT method, and the data are reported as mean relative values compared with the control group.

Table 1. List of primer sequences for RT-qPCR.

Gene	Primer Pair	Product size (bp)
C1qb	Forward-GAACTATGAGCCACGCAACG Reverse-CGAGATTCACACACAGGTTGC	101
Ctss	Forward-TGTTCTTGTGGTTGGCTATGG Reverse-GGTTTAGATTTCTGGGTAAGAGCAATA	163
Itgam	Forward-ACTGAACATCCCATGACCTTCC Reverse-GCTGTAGTCACACTGGTAGAGG	147
Itgb2	Forward-GTGACACTTTACTTGCGACCAG Reverse-CAGCTTCTTGACGTTGTTGAGG	141
Syk	Forward-AAACTACTACAAGGCCCAGACC Reverse-ATAGGAGAACGCTTCCCACATC	136
Tyrobp	Forward-TTAAGTCCCGTACAGGCCCA Reverse-TTGTTTCCGGGTCCCTTCCG	186
Hapln1	Forward-CATCACGGATCTTACCCTGGAG Reverse-TAAGATTGTAGCGTCCCAGTCG	143
Ndst4	Forward-AGAGCCAGTATTCGCAACTTGG Reverse-TACTTCCCTTTGCCACTATCGG	130
Gapdh	Forward-CATCACTGCCACCCAGAAGACTG Reverse-ATGCCAGTGAGCTTCCCGTTCAG	153

Abbreviations: c1qb, complement component 1q beta chain; ctss, cathepsin S; gapdh, glyceraldehyde-3-phosphate dehydrogenase; hapln1, hyaluronan and proteoglycan link protein 1; itgam integrin subunit alpha M; itgb2, integrin subunit beta 2; ndst4, N-deacetylase and N-sulfotransferase 4; syk, spleen tyrosine kinase; tyrobp, TYRO protein tyrosine kinase-binding protein.

Statistical analysis

The RT-qPCR results were evaluated using independent two-tailed Student’s t-tests in GraphPad (version 9.5.1, GraphPad Software Inc., San Diego, CA, USA) to ascertain any differences between the control and EAE-affected groups. All data are presented as the mean (M) ± standard error of the mean (SEM). In all analyses, a p-value less than 0.05 was considered statistically significant.

Results

Histopathological findings

Hematoxylin–eosin staining of the control retinas showed well-preserved retinal layers (Figure 1(A)). However, in the EAE-affected group, there was infiltration and/or activation of inflammatory cells and loss of retinal ganglion cells (RGCs) (Figure 1(B)). Immunohistochemical staining revealed that a few Iba1-positive cells were present in the control retinas (Figure 1(C)). In contrast, more Iba1-positive cells were present in the EAE retinas (Figure 1(D)). These findings indicate that EAE progression is associated with retinal inflammation.

Figure 1. Histopathological examination of the retinas in the control and EAE-affected mice (EAE). (A and B) Hematoxylin-eosin staining revealed infiltration and/or activation of inflammatory cells (B, arrows) and loss of RGCs (B, arrowheads) in the EAE-affected retinas. (C and D) Immunohistochemistry for ionized calcium-binding adapter molecule (Iba1) was performed to detect retinal microglia. A few Iba1-positive resident microglia were observed in the control retinas (C, arrow). However, an increased number of Iba1-positive cells were detected in the EAE-affected retinas (D, arrows). The scale bars in the images represent 20 μm. The retinal layers are labeled as follows: inner limiting membrane (ILM), nerve fiber layer (NFL), ganglion cell layer (GCL), inner plexiform layer (IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer (ONL), outer limiting membrane (OLM), layer of rods and cones (LRC), retinal pigment epithelium (RPE).

Illumina sequencing

To construct each of the eight sample libraries, a total of 529.8 million paired-end (PE) reads with a length of 101 bp were generated. This cumulative raw data equated to 53.5 Gb. Upon the implementation of the read trimming process, which entailed discarding the low-quality segments at the 3’ end of each read as well as filtering out singleton reads and reads shorter than 36 bp, the resulting dataset encompassed 520.6 million paired-end reads, each with a length of 100 bp, totaling 52.2 Gb. The GC (guanine–cytosine) content for the refined dataset exhibited a range of 48–49% across the retina samples from both control and EAE-affected mice (Table 2).

Table 2. Summary statistics of RNA-seq data alignment.

Sample	Total raw read count	Total trimmed read count	Number of drop read after trimming	Number of mapped reads	GC (trimmed) (%)
CON1	67,671,398	66,533,606	1,137,792	65,787,926 (98.88%)	49.51
CON2	62,453,424	61,457,382	996,042	60,805,747 (98.94%)	48.72
CON3	68,917,012	67,936,144	980,868	67,153,686 (98.85%)	48.32
CON4	62,792,462	61,788,244	1,004,218	61,122,200 (98.92%)	48.29
EAE1	66,186,784	64,909,862	1,276,922	64,198,802 (98.9%)	48.45
EAE2	68,687,018	67,436,380	1,250,638	66,765,716 (99.01%)	48.30
EAE3	66,656,828	65,319,138	1,337,690	64,641,732 (98.96%)	48.87
EAE4	66,391,774	65,185,590	1,206,184	64,492,936 (98.94%)	48.49

Abbreviations: CON, control; GC, guanine-cytosine content; Q30, ratio of bases that have Phred quality score of over 30.

DEGs in the retinas of EAE-affected mice

Total RNA-seq was performed during the chronic phase (28 DPI) of EAE to evaluate the gene expression profiles in the retinas of the affected mice. When compared with the expression levels in the control group, 347 genes were differentially expressed in the retinas of EAE-affected mice, with a cutoff value of absolute log₂ fold change > 1 (adjusted p-value < 0.05); 345 were upregulated, whereas only 2 (hapln1 and ndst4) were downregulated (Supplementary Figure S1).

Functional analysis of DEGs

The DEG analysis revealed significant differences in gene expression between the retinas of control and EAE-affected mice. To determine the biological relevance of these DEGs, functional annotation analysis was performed using ShinyGO (Ge et al. 2020). The GO analysis results indicated that the upregulated genes were associated with BPs including immune responses, responses to biotic stimuli, responses to external biotic stimuli, defense responses, responses to other organisms, defense responses to other organisms, and innate immune responses (Figure 2). Figure 2(B) shows the network of enriched gene sets in GOBP terms, where nodes are connected if they share more than 20% genes. The intensity of the color of the nodes reflects the significance of the enriched gene sets, with larger nodes representing bigger gene sets and thicker edges indicating overlapping genes. On the other hand, CCs, such as I bands, Z discs, sarcomeres, myofibrils, contractile fibers, and collagen-containing extracellular matrix (Supplementary Figure S2), and MFs, such as peptide antigen binding, beta-2-microglobulin binding, structural constituents of muscle, antigen binding, extracellular matrix structural constituents, and extracellular matrix binding (Supplementary Figure S3), were enriched in the EAE-affected retinas. Table 3 and Supplementary Table S1 provide detailed information regarding the GO analysis for upregulated genes in the retinas of EAE-affected mice.

Figure 2. GO term enrichment analysis of BPs for upregulated genes in the control and EAE-affected retinas during the chronic phase. (A) The dot plot shows the top ten enriched genes of the GOBP term with a log₂ fold change > 1 and an adjusted p-value < 0.05 in the EAE-affected retinas. (B) The relationship between enriched pathways is represented through an interactive plot, where nodes are connected if they share more than 20% of genes. The intensity of the node color indicates the level of significance of the enriched gene sets.

Table 3. GO term enrichment analysis of upregulated genes in the retinas of EAE-affected mice under GOBP, GOCC, and GOMF.

Pathway	Fold Enrichment	Enrichment FDR	Number of genes	Pathway genes
GOBP
Leukocyte-mediated immunity	4.80	5.41 × 10^−12	34	342
Responses to bacteria	3.69	3.12 × 10^−12	46	745
Innate immune responses	3.69	1.88 × 10^−12	47	742
Defense responses to other organisms	3.39	1.84 × 10^−13	56	1012
Immune responses	2.91	1.58 × 10^−15	80	1619
Responses to external biotic stimuli	2.90	3.37 × 10^−14	72	1403
Responses to other organisms	2.90	3.37 × 10^−14	72	1400
Responses to biotic stimuli	2.89	3.37 × 10^−14	74	1439
Defense responses	2.80	3.37 × 10^−14	76	1546
Biological processes involved in interspecies interaction between organisms	2.66	1.10 × 10^−12	74	1524
GOCC
I band	7.44	1.58 × 10^−15	28	152
Z disc	7.21	1.15 × 10^−13	25	138
Sarcomere	6.75	1.70 × 10^−17	35	213
Myofibril	6.44	1.70 × 10^−17	36	234
Contractile fiber	6.15	5.79 × 10^−17	36	245
Collagen-containing extracellular matrix	4.90	1.58 × 10^−15	40	362
External encapsulating structure	4.17	7.39 × 10^−14	42	492
Extracellular matrix	4.17	7.39 × 10^−14	42	490
Extracellular space	2.92	9.19 × 10^−16	74	1658
Cell surface	2.75	1.01 × 10^−9	51	926
GOMF
Peptide antigen binding	14.44	3.20 × 10^−5	7	17
Beta-2-microglobulin binding	13.41	1.21 × 10^−3	5	12
Structural constituent of muscle	10.06	1.21 × 10^−3	6	21
Antigen binding	7.22	2.42 × 10^−3	7	42
Extracellular matrix structural constituent	6.29	5.81 × 10^−8	19	136
Extracellular matrix binding	6.24	4.20 × 10^−4	10	57
Actin filament binding	3.46	1.21 × 10^−3	16	206
Structural molecule activity	2.67	3.16 × 10^−5	35	612
Actin binding	2.56	1.21 × 10^−3	25	443
Calcium ion binding	2.44	4.20 × 10^−4	32	709

Abbreviations: FDR, false discovery rate; GOBP, gene ontology biological process; GOCC, gene ontology cellular component; GOMF, gene ontology molecular function.

RT-qPCR validation of DEGs in EAE-affected retinas

To determine if the gene expression detected by RT-qPCR was similar to that detected through the RNA-seq analyses, hub genes and downregulated genes were analyzed based on gene expression and biological significance. Twelve hub genes were identified in the EAE-affected retinas by analyzing the PPI networks of DEGs using the STRING database (Figure 3 and Supplementary Figure S4). To validate the identified gene expression, six hub genes and two significantly downregulated genes among the DEGs were selected for RT-qPCR validation (Figure 4). The validated genes were complement component 1q beta chain (c1qb), cathepsin S (ctss), integrin subunit alpha M (itgam), integrin subunit beta 2 (itgb2), spleen tyrosine kinase (syk), TYRO protein tyrosine kinase-binding protein (tyrobp), hyaluronan and proteoglycan link protein 1 (hapln1), and N-deacetylase and N-sulfotransferase 4 (ndst4).

Figure 3. PPI network for the DEGs in the EAE-affected retinas. The PPI network was generated using the STRING database, with a confidence value of 0.7. The nodes in the network represent the genes, with larger nodes indicating hub genes. The hub genes c1qa, c1qb, ctss, itgam, itgb2, lyn, ptprc, spi1, syk, tlr2, ttn, and tyrobp are shown in brown-colored nodes that are larger than others. The abbreviations for the genes are as follows: c1qa (complement component 1q alpha chain), c1qb (complement component 1q beta chain), ctss (cathepsin s), itgam (integrin subunit alpha M), itgb2 (integrin subunit beta 2), lyn (lyn proto-oncogene src family tyrosine kinase), ptprc (protein tyrosine phosphatase receptor type c), spi1 (spleen focus forming virus proviral integration oncogene), syk (spleen tyrosine kinase), tlr2 (toll-like receptor 2), ttn (titin), and tyrobp (TYRO protein tyrosine kinase-binding protein).

Figure 4. Expression changes in eight selected DEGs in the retinas of mice with EAE. The bar graphs depict the relative expression levels of six upregulated (A) and two downregulated (B) genes based on RNA-seq data. The data are presented as the mean ± SEM (n = 7/group). Statistical significance is denoted by * p < 0.05, ** p < 0.01, and *** p < 0.001 compared with the control group.

The validation of gene expression revealed that in the EAE-affected retinas, the expression levels of c1qb (M = 2.07, SEM = 0.32, t(12) = 3.298, p = 0.006), ctss (M = 3.10, SEM = 0.84, t(12) = 2.479, p = 0.029), itgam (M = 4.37, SEM = 1.10, t(12) = 3.040, p = 0.01), itgb2 (M = 11.45, SEM = 1.66, t(12) = 6.271, p < 0.001), syk (M = 1.92, SEM = 0.41, t(12) = 2.251, p = 0.044), and tyrobp (M = 5.45, SEM = 1.95, t(12) = 2.279, p = 0.041) were significantly increased (Figure 4(A)), while hapln1 (M = 0.21, SEM = 0.22, t(12) = 3.550, p = 0.004) and ndst4 (M = 0.31, SEM = 0.17, t(12) = 3.835, p = 0.002) were significantly decreased (Figure 4(B)). These findings indicate that the RT-qPCR validation was consistent with the RNA-seq analysis and that the selected genes are likely implicated in the pathogenesis of EAE.

Additional assessment of functional differences using GSEA

GSEA was performed with the MSigDB curated gene sets (FDR q-value < 0.05, NOM p-value < 0.05) and showed that the top 10 positively enriched gene sets in EAE-affected retinas were related to several processes including the electron transport chain, respiratory electron transport ATP synthesis by chemiosmotic coupling and heat production by uncoupling proteins, citric acid TCA cycles and respiratory electron transport, respiratory electron transport, graft versus host disease D7 up, glis2 target up, nonsense-mediated decay independent of the exon junction complex, cytoplasmic ribosomal proteins, SRP-dependent cotranslational protein targeting to the membrane, and oxidative phosphorylation (Figure 5 and Supplementary Figure S5).

Figure 5. Curated gene sets enriched in the retinas of mice with EAE. The bar graphs depict the top 10 GSEA curated gene sets (FDR q-value < 0.05, NOM p-value < 0.05) that were positively enriched (shown in red bars) and negatively enriched (shown in green bars) in the EAE-affected retinas. The height of each bar in panels A and B indicates the strength of the normalized enrichment score (NES) and the size of the gene set, respectively.

On the other hand, the top 10 negatively enriched gene sets included processes related to cortical plate postmitotic neurons, neuron markers, transmission across chemical synapses, the neuronal system, protein–protein interactions at synapses, the phototransduction cascade, splicing factor NOVA-regulated synaptic proteins, subventricular zone and up neuron fate-committed cells, the neurotransmitter release cycle, and targets of ptch1 and sufu DN (Figure 5 and Supplementary Figure S6). Table 4 and Supplementary Table S2 provide detailed reports on the GSEA with curated gene sets from the retinas of EAE-affected mice.

Table 4. GSEA results according to the MSigDB curated gene sets.

MSigDB gene set	Size	NES	FDR q-Value	NOM p-Value
Positively enriched gene sets
Electron transport chain	82	2.78	0	0
Respiratory electron transport ATP synthesis by chemiosmotic coupling and heat production by uncoupling proteins	100	2.68	0	0
Citric acid TCA cycle and respiratory electron transport	131	2.68	0	0
Respiratory electron transport	80	2.61	0	0
Graft versus host disease D7 up	84	2.59	0	0
Glis2 targets up	68	2.58	0	0
Nonsense-mediated decay independent of the exon junction complex	72	2.58	0	0
Cytoplasmic ribosomal proteins	68	2.57	0	0
SRP dependent cotranslational protein targeting to membrane	72	2.56	0	0
Oxidative phosphorylation	44	2.55	0	0
Negatively enriched gene sets
Cortical plate postmitotic neurons	86	−3.06	0	0
Neuron markers	53	−2.95	0	0
Transmission across chemical synapses	157	−2.88	0	0
Neuronal system	224	−2.86	0	0
PPIs at synapses	39	−2.67	0	0
Phototransduction cascade	23	−2.66	0	0
Splicing factor NOVA-regulated synaptic proteins	34	−2.56	0	0
Subventricular zone and up neuron fate-committed cells	35	−2.53	0	0
Neurotransmitter release cycle	36	−2.52	0	0
Targets of ptch1 and sufu DN	64	−2.48	0	0

Abbreviations: ATP, adenosine triphosphate; D7, day 7; DN, down; MSigDB, Molecular Signatures Database; NOVA, neuro-oncological ventral antigen; FDR, false discovery rate; NES, normalized enrichment score; NOM, nominal; PPI, protein-protein interaction; SRP, signal recognition particle; TCA, tricarboxylic acid cycle.

Discussion

This study aimed to identify potential molecular pathways associated with visual impairment in EAE, an animal model of MS, by analyzing changes in gene expression in the retinas of mice. Using RNA-seq, the retinal transcriptomes of control and EAE-affected mice were compared and revealed a significant upregulation of immune response- and inflammation-related genes among DEGs. Intriguingly, the downregulated DEGs were predominantly represented by only two genes, both of which are related to retinal structure and photoreceptor function.

An interesting observation during this study was the apparent asymmetry in the quantities of upregulated and downregulated genes. This phenomenon could potentially be attributed to either the selection of the statistical significance threshold or the inherent biological characteristics of the EAE model itself. Notably, the decision to set a relatively stringent fold-change threshold (with an absolute log₂ fold change > 1) to identify significantly involved DEGs in the retinas of EAE mice might underlie the fewer observed downregulated genes than upregulated genes within the DEG dataset. In addition, a previous study showed that there was an imbalance in the expression of immune- and neural-related genes associated with the progression of the disease in non-obese diabetic mice with EAE (Baranzini et al. 2005). This existing knowledge could serve to elucidate the seemingly asymmetric expression of upregulated and downregulated genes observed within the current DEG dataset.

Based on the results of our study, the c1qb, itgam, and itgb2 hub genes were upregulated, which might be associated with the common characteristics of EAE and MS pathogenesis, particularly T-cell mediated inflammation and demyelination. It has been suggested that the complement system is involved in the pathology of MS and EAE (Hammond et al. 2020; Park and Jung 2022). Various investigations have highlighted increased C1q expression in cerebral white matter lesions and spinal cord regions (Ingram et al. 2014) as well as in the hippocampi of subjects affected by MS and EAE (Michailidou et al. 2015; Hammond et al. 2020). It was additionally reported that macrophage-1 antigen, which is composed of the complement receptor subunits ITGAM (CD11b) and ITGB2 (CD18), is implicated not only in myelin phagocytosis mediated by microglia and macrophages but also in T-cell-driven demyelination in EAE (Bullard et al. 2005). In accordance with these findings, our results consistently align with the characteristic pathology of both EAE and MS, particularly in the context of complement system involvement.

Prior investigations have delved into the analysis of DEGs in various tissues such as retinal inflammatory cells (Cruz-Herranz et al. 2021), hippocampi (Weerasinghe-Mudiyanselage et al. 2022), and olfactory bulbs (Kim et al. 2019) in EAE animal models as well as peripheral blood mononuclear cells of patients with MS (Tuller et al. 2011). Nonetheless, the present study is the first to focus exclusively on DEGs within EAE-affected retinas. Notably, within this study, the upregulation of tyrobp expression in the retinas of EAE-affected mice parallels its upregulation in the hippocampi, suggesting its involvement across select nervous tissues. In a distinctive context, the downregulation of the expression of ndst4 and hapln1 appears to be exclusive to the retinas of EAE mice, highlighting potential unique roles for these genes in the retinal response.

Microglial infiltration/activation and loss of RGCs were also observed in the histopathological analysis, which is consistent with previous findings in EAE (Manogaran et al. 2019). Previous studies have also shown the involvement of innate immunity and molecular changes in innate immune cells in the retinas during EAE progression (Cruz-Herranz et al. 2021). However, the molecular mechanisms responsible for retinal tissue damage in EAE-affected mice remain poorly understood. Therefore, we attempted to identify the molecular alterations in EAE-affected retinas that are responsible for EAE-related visual impairments.

This study utilized GO, GSEA, and STRING analyses to investigate the functional changes in the retinas of mice with EAE. Previous studies have suggested that retinal innate immune responses are significantly involved in EAE pathogenesis and tissue damage (Cruz-Herranz et al. 2021), which is consistent with our findings. GO analysis of DEGs in EAE-affected retinas revealed significant enrichment in immune system processes, including immune responses, leukocyte-mediated immunity, and innate immune responses. GSEA showed that the most positively enriched gene sets were related to respiratory electron transport and oxidative phosphorylation, which are associated with mitochondrial pathology in EAE-affected retinas. Previous studies have linked axonal degeneration to mitochondrial changes and ROS production caused by oxidative stress in the chronic phases of MS and EAE (van Horssen et al. 2011), and macrophage-derived ROS can cause mitochondrial pathology in EAE-affected spinal cords, leading to focal axonal degeneration (Nikic et al. 2011). A significant increase in mitochondrial respiratory chain complex activity has been reported in MS remyelinated axons (Zambonin et al. 2011). Conversely, the neuronal system, PPIs at synapses, and the phototransduction cascade were negatively enriched in the GSEA. The PPIs in the synapse pathway, which contains synaptic adhesion-like molecule (SALM)-related genes, such as synaptic adhesion-like molecule 1-3 genes (salm1-3) and post-synaptic density 95 (psd95) genes, were negatively enriched in the GSEA. The SALM family is a leucine-rich repeat-containing synaptic adhesion molecule subfamily that promotes neurite outgrowths and has synaptogenic properties (Nam et al. 2011). The SALM family is involved in interactions with PSD95 (Morimura et al. 2006), whose expression was decreased in the inner retinas of EAE mice compared with control mice (Jin et al. 2019). The GSEA also revealed negative enrichment of genes involved in the phototransduction cascade pathway, such as ATP-binding cassette subfamily A member 4 (abca4) and phosphodiesterase 6B (pde6b) genes. ABCA4 protein, which is expressed in the photoreceptor outer segment disc membrane, is associated with retinal dystrophy (Kim et al. 2022), and phosphodiesterase 6 (PDE6) enzymes, which are located in both rod and cone photoreceptor cells, are crucial for the phototransduction cascade. Dysfunction of PDE6 activity causes retinal dystrophies in retinal degeneration 1, a mouse animal model of retinitis pigmentosa (Sahaboglu et al. 2013). These negatively enriched pathways in the retinas of EAE-affected mice may contribute to dysfunction of the neuronal system and phototransduction in the retinas.

In this study, the retinas of EAE-affected mice have been found to contain several novel genes including ndst4 and halpln1 that may contribute to visual dysfunction. Among the immune response- and inflammation-related genes identified in the retina, certain genes have been identified as hub genes that play significant roles in retinal inflammation.

One such gene is c1q, which is involved in the classical complement pathway and plays a role in target recognition (Schartz and Tenner 2020). In particular, C1q is a critical component of the complement cascade that is involved in synaptic remodeling (Zabel and Kirsch 2013). C1q is also involved in neuroprotection (Benoit et al. 2012) and axon regeneration (Peterson et al. 2015). In EAE-affected optic nerves, gene expression of c1q was found to be increased in activated microglia (Jin et al. 2019). In other conditions, such as retinal ischemia/reperfusion injury (Kuehn et al. 2008) and diabetes (Church et al. 2022), c1qb expression has also been reported to be upregulated in the retina. While C1q promotes and degrades damaged RGCs and removes cell debris from apoptotic RGCs (Kuehn et al. 2008), it has been reported that complement-mediated synapse elimination is abnormally reactivated in neurodegenerative diseases of the CNS and retina (Stevens et al. 2007). The increased expression of c1qb observed in this study might be related to downregulation of the synaptic-related pathways in the GSEA. Further investigation is needed to better elucidate the mechanisms underlying the abnormal reactivation of synaptic elimination in relation to c1qb.

Cathepsin S (ctss) plays an essential role in antigen presentation, being involved in both the processing of antigens and the processing of the Invariant chain (li; CD74), which can act as an MHC class II chaperone (Anderson and Miller 1992; Riese et al. 1996). Cathepsin S has been largely linked to inflammatory responses, oxidative phosphorylation, coagulation, and epithelial–mesenchymal transition, as shown in diabetic retinopathy (Someya et al. 2022). In models of autoimmunity, inhibition of cathepsin S causes both suppression of inflammation and disease chronicity index (Tato et al. 2017). While cathepsin S inhibitors have been investigated as treatments for autoimmune diseases (Baugh et al. 2011), it remains unclear whether they can reduce the severity of visual impairment associated with MS and EAE. The increased expression of cathepsin S observed in this study is likely related to inflammation, oxidative phosphorylation, and endothelial-mesenchymal transition. Identifying the mechanisms underlying visual impairments is necessary to determine whether cathepsin S affects them.

Integrins, such as itgam and itgb2, have been identified as inflammatory factors through cell migration and adhesion (Hynes 2002). Upregulation of these genes has been reported in diabetic retinopathy in Zucker diabetic fatty rats, as seen in the CNS inflammation, including spinal cord injury (Lund et al. 2022) and hemorrhage (Weng et al. 2022). Therefore, upregulation of itgam and itgb2 might promote leukocyte recruitment in the retina, which might be a therapeutic target.

Tyrobp, a hub gene in the hippocampus in EAE (Weerasinghe-Mudiyanselage et al. 2022), has also been found to be upregulated in EAE-affected retinas. Triggering receptor expressed on myeloid cells 2/DNAX activating protein of 12 kDa (TREM2/DAP12) signaling promotes the proliferation, migration, and phagocytosis of apoptotic neurons in microglia (Konishi and Kiyama 2018). In addition, upregulation of TREM2/DAP12 via SYK increases proangiogenic cytokines in Müller cells, which contributes to diabetic retinopathy progression (Yerlikaya et al. 2022). Both upregulation of tyrobp and syk might be involved in the glial activation and be associated with chronic retinopathy in EAE.

On the other hand, only two downregulated genes were identified in the present study, namely, ndst4 and hapln1; ndst4 encodes for heparan sulfate N-deacetylase/N-sulphotransferase (NDST) enzymes, which are responsible for the sulphation and modification of heparan sulfate (Aikawa et al. 2001). Heparan sulfate plays a crucial role in regulating several interleukins and chemokines, which are known mediators of inflammation (Farrugia et al. 2018). It is also essential for the projection of axons from RGCs to the optic nerve head (Coulson-Thomas 2016). In addition, the downregulation of ndst4, which has been reported in cases of retinal ischemia (Neo et al. 2020) and in vitro model of ALS (Lo et al. 2021), may contribute to the inflammatory process and dysregulation of RGC axon guidance.

Furthermore, hapln1 encodes for hyaluronan and proteoglycan link protein 1 (HAPLN1), which is abundant in the interphotoreceptor matrix of the human retina (Ishikawa et al. 2015). Studies have suggested that the interphotoreceptor matrix may play a role in retinal degeneration. Downregulation of hapln1 has been previously reported to be involved in dysfunction of RGCs (Rocha et al. 2020) and Müller cells (Wang et al. 2013) in the retina. However, the effect of hapln1 downregulation in the retinas of mice with EAE remains unclear. The disruption of the interphotoreceptor matrix, which contains proteoglycans such as hapln1 (Ishikawa et al. 2015), may be responsible for the downregulation of hapln1 and may contribute to visual impairments resulting from retinal degeneration.

In conclusion, this study utilized RNA-seq data analysis to identify transcriptomic profiling and hub genes in the retinas of mice with EAE. The upregulation of c1qb, ctss, tyrobp, syk, itgam, and itgb2, which are involved in immune responses, and downregulation of ndst4 and hapln1, which are involved in dysfunction of RGCs, may contribute to visual impairments in EAE-affected mice. These genes could serve as potential targets for future therapeutic studies on visual impairment in EAE. The findings of this study provide novel insights into the molecular mechanisms underlying visual disturbances in EAE and highlight potential therapeutic targets for future investigation.

CRediT authorship contribution statement

SH, CM, and TS participated in research design. SH, PDEWM and SK performed the study. SH, CM, and TS contributed in analysis and writing first draft of manuscript. All authors read and approved the final manuscript.

Data availability

Data will be made available on request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was supported by the National Research Foundation of Korea (NRF- 2019R1A2C1087753).