All animal experiments conducted in this study were approved in advance by Northeast Agricultural University. FGF21 knockout mice (FGF21 KO mice, KO group) were purchased from Cyagen Biosciences Inc. (Suzhou, Contract No. KOAIP180401WZ1+KOAIP180401WZ2). C57BL/6 mice of equivalent age (WT group) were obtained from Liaoning Chang Sheng Biotechnology Co., Ltd. (Liaoning, China). All animals were housed in a controlled environment with a 12-hour light/12-hour dark cycle, provided with ad libitum access to food and water, and shielded from external stimuli. The ambient temperature for housing was maintained at 22 $\pm$ 2 ℃. A total of 36 female mice were used in this study. Prior to the OCT imaging, mice were anesthetized with Avertin (1.25% w/v solution, 30 µL/g, MA0478, Dalian Meilun Biotechnology Co., Ltd, Dalian, Liaoning, China). At the end of the experiment, all animals were humanely euthanized. Euthanasia was performed by cervical dislocation under deep anesthesia (induced by the same Avertin solution (1.25% w/v, 30 µL/g, i.p.) used for OCT imaging). The absence of a pedal reflex was confirmed to ensure death before any tissue collection. Then the retinal/choroidal tissues were harvested for RNA-seq and other analyses.

To determine whether endogenous FGF21 deficiency affects retinal resilience under pathological stress, mice were allocated into two independent disease-modeling experiments (n = 6 per genotype per model). For diabetic retinopathy (DR) Model, because of the hypersensitivity of FGF21 KO mice to streptozotocin (STZ) (HY-13753, MedChemExpress, Shanghai, China), diabetes was induced using genotype-specific protocols. WT mice received intraperitoneal injections of STZ (30 mg/kg/day for 5 consecutive days) while maintained on a high-fat diet (HFD, 1135DM). FGF21 KO mice were fed the same HFD without STZ injection until their fasting blood glucose levels exceeded 11.1 mmol/L. Diabetic mice (fasting glucose $>$11.1 mmol/L from both groups) were maintained for an additional 8 weeks before analysis. For dry Age-related Macular Degeneration (dAMD) Model, both FGF21 KO and WT mice were allowed to freely consume drinking water containing 0.8% hydroquinone for three months to induce oxidative retinal damage. Optical coherence tomography (OCT) was performed at the endpoint of each model to assess retinal structural integrity (n = 6).

To investigate the direct effect of endogenous FGF21 deletion on retinal homeostasis in the absence of external stressors, a separate cohort of mice (n = 6 per genotype, WT and FGF21 KO) was examined at 6 weeks of age under normal physiological conditions. After baseline OCT imaging, mice were euthanized, and neural retina/choroid complexes were promptly harvested. Tissues from each group were pooled for bulk RNA sequencing (n = 3 mice), and other mice per group were used for subsequent qPCR validation (n = 3 retina/choroid complexes) and immunohistochemistry (n = 3 retina/choroid complexes).

Eyeballs were fixed with 4% paraformaldehyde, dehydrated, cleared, and embedded in paraffin. Sections (4–5 µm) were dewaxed and subjected to heat-induced epitope retrieval (HIER). After blocking endogenous peroxidase activity with 3% H₂O₂ and nonspecific sites with 10% horse serum, sections were incubated with primary antibodies against CD163 (bs2527R, 1:100, Bioss, Beijing, China). Following PBS washes, sections were incubated with a Goat Anti-Rabbit IgG H&L antibody (bs-0295G-HRP, 1:200, Bioss, Beijing, China), developed with 3,3-diaminobenzidine (DAB) (HY-W014212, MedChemExpress, Shanghai, China), and counterstained with hematoxylin (C0107, Beyotime, Shanghai, China).

For semi-quantification of CD163 staining, whole-section images were analyzed using ImageJ software (NIH, Bethesda, USA) with the “IHC Toolbox” plugin. The integrated optical density (IOD) of the DAB (brown) signal, which reflects the amount of CD163 protein within the tissue, was measured in a defined region of the entire retinal section.

Following a 12-h fast, mice were weighed and baseline blood glucose was measured prior to oral gavage of glucose (2 g/kg). Blood glucose levels were then monitored at 30-minute intervals for 120 minutes. Glucose tolerance was determined by calculating the area under the curve (AUC) for blood glucose over time using Origin 2021 (OriginLab Corporation, Northampton, MA, USA) software.

Mice were anesthetized via intraperitoneal Avertin (1.25% w/v solution, 30 µL/g, MA0478, Dalian Meilun Biotechnology Co., Ltd, Dalian, Liaoning, China). Pupillary dilation was achieved using tropicamide-phenylephrine eye drops (1–2 drops per eye, Santen Pharmaceutical Co., Ltd, Suzhou, Jiangsu, China). Throughout the procedure, ophthalmic hydroxypropyl methylcellulose (UV-Gel) (Shanghai Haohai Biological Technology Co., Ltd., Shanghai, China) was intermittently applied to the eyes of the mice to safeguard the cornea. OCT images were captured using the VG200D imaging camera system (SVision imaging, Ltd, San Jose, CA, USA).

RNA-Seq testing was commissioned to the Biotechnology Company (Shanghai, China).

Total RNA was extracted from tissues using the TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA). Then, the first-strand cDNA was synthesized using a HiScript II 1st Strand cDNA Synthesis Kit (R212-01, Vazyme, Nanjing, Jiangsu, China). The gene expression was detected with a Taq Pro Universal SYBR qPCR Master Mix Kit (Q712-02, Vazyme, Nanjing, Jiangsu, China), according to the manufacturer’s protocol. The expression for each gene was calculated using the expression 2^-⁣△△Ct method. The Real-time PCR primers used in this study were synthesized by Sangon Biotech (Shanghai) Co., Ltd. The primer sequences are shown in Table 1.

DESeq was employed for the analysis, with the screening criteria set as p value 0.05 and a multiplicity of difference $|$log2FoldChange$|$ $>$1 to identify significantly different genes. The results of the expression difference analysis were visualized using BioLadder (https://www.bioladder.cn). Disease-related gene sets were retrieved from three public databases: GeneCards (https://www.genecards.org), DisGeNET (https://www.disgenet.org), and the Therapeutic Target Database (TTD, https://db.idrblab.net/ttd).

The relevant genes obtained from the screening were imported into the STRING website (http://string-db.org) to construct the protein-protein interaction (PPI) network [12]. The results of the PPI network were then visualized and enhanced using Cytoscape (version 3.10.0, https://cytoscape.org, National Resource for Network Biology, USA). To comprehensively analyze the PPI network and identify key components, the cytoHubba [13] plugin was utilized. Specifically, the Maximal Clique Centrality (MCC) algorithm was employed to screen for key genes, facilitating the construction of sub-networks and identification of pivotal components within the PPI network.

The analysis was carried out by using R software (v.4.2.2) package clusterProfiler (v.4.5.0) [14] through Hiplot Pro (https://hiplot.com.cn/), a comprehensive web service for biomedical data analysis and visualization.

The correlation coefficients between the variables were computed through Pearson correlation analysis. A grid plot was generated to display the genes in the group with an absolute correlation coefficient greater than 0.8 and a p-value less than 0.05. This facilitated the visualization of interactions between significantly correlated feature nodes, aiding in the identification of relationships between the variables. The process of calculation and plotting was conducted using CNSknowall (https://cnsknowall.com/).

Immune Cell Abundance Identifier (ImmuCellAI) is a tool to accurately estimate immune cell abundance from gene expression datasets, including RNA-Seq [15]. In the present study, the Immune Cell Abundance Identifier for mouse (ImmuCellAI-mouse) [16] is used to estimate the abundance of 36 immune cell (sub) types in mouse retina/choroid RNA-Seq data.

Data are presented as the mean $\pm$ standard error of the mean (SEM). The sample size (n) for each experiment represents biological replicates. Statistical comparisons between two groups were performed using an unpaired, two-tailed Student’s t-test. For comparisons among multiple groups, one-way analysis of variance (ANOVA) followed by appropriate post-hoc tests was used. A p 0.05 was considered statistically significant.

To assess the intrinsic role of endogenous FGF21, we first examined retinal structure and systemic metabolism in unchallenged, 6-week-old mice. In vivo spectral-domain optical coherence tomography (OCT) revealed comparable retinal thickness, reflectivity, and overall architecture between wild-type (WT) and FGF21 knockout (KO) mice (Fig. 1A,B). Furthermore, FGF21 KO mice showed no significant differences in body weight, fasting blood glucose, or glucose tolerance compared to WT littermates (all p $>$ 0.05; Fig. 1G–I).

Fig. 1.

Fibroblast Growth Factor 21 (FGF21) knockout exacerbates retinal structural damage induced by hyperglycemia and oxidative stress in mouse models. (A–F) OCT detection results of the fundus of mice. (A) 6-week-old WT mice. (B) 6-week-old FGF21 KO mice. (C) WT diabetic mice. (D) FGF21 KO diabetic mice. (E) WT AMD mice. (F) FGF21 KO AMD mice. Green lines indicate the position of the OCT scan line. White arrows indicate drusen, the typical lesions of dry age-related macular degeneration. (G) Oral glucose tolerance test (OGTT) performed in 6-week-old FGF21 KO and WT mice. (H) Area under the curve (AUC) of blood glucose during the OGTT. (I) Body weight measurements of 6-week-old FGF21 KO and WT mice. KO: FGF21 knockout; WT: normal control; ns: p $>$ 0.05.

We next investigated whether the absence of FGF21 affects retinal resilience to pathological insults. Quantitative analysis of retinal layer thickness by OCT revealed that FGF21 KO mice exhibited significantly exacerbated structural damage under both pathological conditions.

When subjected to diabetic conditions or hydroquinone-induced oxidative stress (a model for dry age-related macular degeneration), FGF21 KO mice exhibited significantly more severe retinal structural damage than their WT counterparts, as quantified by OCT (Fig. 1C–F). Retinal thickness measurements were performed on the averaged image from three consecutive scans. Using the segmented line tool, three perpendicular measurements were taken from the internal limiting membrane (ILM) to the retinal pigment epithelium (RPE) to determine the total retinal thickness. The values were then averaged to yield a single thickness value per eye.

Under diabetic conditions, the total retinal thickness in FGF21 KO mice was thinner than in WT controls (p 0.05). Similarly, in the hydroquinone-induced AMD model, FGF21 KO mice exhibited more severe thinning (Supplementary Fig. 1). In addition to thickness, OCT analysis revealed discrete hyperreflective foci and irregularities in the Bruch membrane area of FGF21 KO mice, suggestive of drusen-like deposits (white arrow). These results indicate that endogenous FGF21 is crucial for maintaining retinal homeostasis, and its deficiency predisposes the tissue to injury under stress.

Transcriptomic analysis of retinal/choroidal complexes from 6-week-old FGF21 KO and WT mice revealed 449 differentially expressed genes (DEGs) (p 0.05, $|$log2FoldChange$|$ $>$1), with 293 up-regulated and 156 down-regulated in the KO group (Fig. 2A). Principal component analysis (PCA) demonstrated a clear separation between the genotypes along the first principal component (PC1, 95.7% of variance), indicating a fundamental shift in the global transcriptional landscape upon FGF21 deficiency (Fig. 2B). Unsupervised hierarchical clustering heatmaps of the DEGs also distinctly segregated the WT and KO samples (Fig. 2C).

Fig. 2.

Transcriptomic profiling reveals differential gene expression patterns in the retina/choroid of FGF21 KO mice compared to WT controls. (A) Volcano plot illustrating differentially expressed genes (DEGs) between KO and WT groups. Significantly up-regulated genes are highlighted in red, while significantly down-regulated genes are shown in green. Blue dots represent non-significantly expressed genes. (B) Principal component analysis (PCA) plot demonstrating the clustering pattern of samples based on the global transcriptome profiles, with PC1 and PC2 indicating the first and second principal components, respectively. (C) Hierarchical clustering heatmap of differentially expressed mRNA. Rows represent individual genes, and columns represent individual samples. Color scale indicates normalized expression levels, with red denoting high expression and blue denoting low expression. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs. KO: FGF21 knockout; WT: normal control.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of these DEGs highlighted significant involvement in immune and inflammatory pathways. The most enriched terms included Cytokine-cytokine receptor interaction signaling pathway, Motor proteins, Alcoholic liver disease, C-type lectin receptor signaling pathway, and TNF signaling pathway (Fig. 2D).

To explore the disease relevance of the transcriptomic changes, we intersected the 449 DEGs with known fundus disease-associated genes from the GeneCards database (relevance score $>$1; 3401 genes). This yielded 51 overlapping genes (Fig. 3A). A protein-protein interaction (PPI) network constructed from these 51 genes using the STRING database (species: Mus musculus) revealed significant interconnectivity (Fig. 3B). Further analysis of this network in Cytoscape using the Maximal Clique Centrality (MCC) algorithm identified six hub genes: tumor necrosis factor (Tnf), interleukin 1$\beta$ (Il1b), prostaglandin-endoperoxide synthase 2 (Ptgs2), peroxisome proliferator-activated receptor $\gamma$ (Pparg), leptin (Lep), and C-X-C motif chemokine 1 (Cxcl1) (Fig. 3C). The MCC algorithm scores for these genes were 14,711, 14,574, 13,752, 13,138, 11,658, and 11,526, respectively.

Fig. 3.

Integration analysis identifies FGF21 deficiency-associated retinal/choroidal differentially expressed genes (DEGs) overlapping with ocular fundus disease signature genes and their functional interaction networks. (A) Venn diagram depicting the intersection between retinal/choroidal DEGs from FGF21 KO mice and ocular fundus disease-associated genes curated from the GeneCards database. (B) PPI network diagram. Nodes represent individual proteins, and edges indicate predicted functional interactions. (C) Identification of hub genes from the PPI network using the Maximal Clique Centrality (MCC) algorithm implemented in CytoHubba plugin of Cytoscape. The top-ranked core genes are displayed in descending order of MCC scores, with color intensity proportional to the centrality ranking. (D) KEGG versus GO enrichment analysis. KO: FGF21 knockout; WT: normal control.

Functional enrichment analysis of these hub genes highlighted their involvement in key pathways and biological processes. KEGG pathway analysis showed enrichment in the TNF signaling pathway, Non-alcoholic fatty liver disease, and Leishmaniasis. Gene Ontology (GO) analysis indicated associations with the Acute inflammatory response, Temperature homeostasis, Photoreceptor inner segment, Keratin filament, Cornified envelope, RNA polymerase II core promoter sequence-specific DNA binding (Fig. 3D).

To delineate the specific relevance of FGF21 deficiency to age-related macular degeneration (AMD) and diabetic retinopathy (DR), we intersected our DEG list with disease-associated genes from GeneCards, DisGeNET, and Therapeutic Target Database. After applying a relevance score filter ($>$10 for GeneCards), 33 genes were common to the FGF21 KO DEGs and the combined AMD/DR gene sets (Fig. 4A).

Fig. 4.

Gene enrichment analysis of age-related macular degeneration (AMD)/diabetic retinopathy (DR) associated DEGs between KO and WT groups. (A) Venn diagrams depicting DR and AMD-related genes with significantly different genes from three databases: GeneCards, DisGeNET, and Therapeutic Target Database. (B) Correlation between disease (AMD and DR) related genes, FGF receptors, and FGF19 subfamily members. (C) Chord plot illustrating GO enrichment analysis of DR-related differential genes. (D) Chord plot illustrating GO enrichment analysis of AMD-associated differential genes. Differential genes are presented on the left, with color changes indicating alterations in log2FC from small to large, with red indicating up-regulated and blue indicating down-regulated genes. Enrichment term results are displayed on the right, and the two are linked to indicate gene enrichment into that term. KO: FGF21 knockout; WT: normal control.

Pearson correlation analysis between these 33 candidate genes and members of the FGF19 subfamily (including FGF21) and their receptors (FGFRs) revealed several significant associations: Mlxipl and Aoc3 showed strong positive correlations with FGFR1, Fos, Retn, and Adipoq correlated positively with FGFR4, and Fos also correlated with FGF15. In contrast, Cd5l exhibited a strong negative correlation with Fgf23 (Fig. 4B–D). Functional enrichment analysis of these 33 overlapping genes highlighted acute inflammatory response as a key biological process linking the FGF21-deficient retinal phenotype to both AMD and DR pathogenesis (Fig. 5A,B).

Fig. 5.

Screening of core DEGs related to AMD/DR in the retina/chorioid of FGF21 KO mice. (A) STRING diagram illustrating KEGG enrichment analysis of DR-related differential genes. (B) STRING diagram illustrating KEGG enrichment analysis of AMD-related differential genes. On the left side are the differential genes, with color changes indicating the change in log2FC from small to large, where red represents up-regulation and blue represents down-regulation. On the right side are the enrichment term results, and the two are connected to indicate gene enrichment into that term. (C) Real-time PCR analysis for differential genes in indicated groups. (D) PPI network diagram. (E) Core differential genes identified by the MCC algorithm. (F) MCC scores and rankings of related genes. KO: FGF21 knockout; WT: normal control. Statistical significance is indicated as * p 0.05 and ** p 0.01.

In addition, we constructed a PPI network from these genes using STRING and Cytoscape, which comprised 27 nodes and 109 interactionsFunctional enrichment analysis of these 33 overlapping genes (Fig. 5D). Application of the MCC algorithm identified five central hub genes: tumor necrosis factor (Tnf), interleukin 1$\beta$ (Il1b), peroxisome proliferator-activated receptor $\gamma$ (Pparg), chemokine (Ccl5), and C-X-C motif chemokine 1 (Cxcl1) (Fig. 5E,F). Notably, Pparg, Il1b, and Cxcl1 were upregulated in FGF21 KO mice, while Tnf and Ccl5 were downregulated. The alterations in gene transcription were further validated through qPCR (Fig. 5C).

To investigate the effect of FGF21 deficiency on the immune microenvironment, we estimated immune cell abundances from the RNA-seq data using the ImmuCellAI-mouse. The results revealed a significant shift in immune cell composition in FGF21 KO retinal/choroidal tissues. Specifically, the abundance of M2 macrophages was markedly increased, whereas NKT cells and germinal center B cells were decreased compared to WT controls (Fig. 6A,B).

Fig. 6.

Immune infiltration analysis. (A) Relative percentages of 36 infiltrating immune cells in the KO group versus the WT group. (B) Box plot of 36 infiltrating immune cells in the KO group versus the WT group. (C) Bubble plots showing the correlation between differential genes and immune cells. Bubble colors represent the magnitude of correlation coefficients, indicating positive and negative correlations. The size of the bubbles represents the magnitude of the p-value, signifying the significance level of this sample. (D) Representative images of IHC staining for CD163. (E) Quantification of CD163 staining. Scale bar = 50 µm. *, p 0.05; **, p 0.01. KO: FGF21 knockout; WT: normal control.

Subsequently, Pearson correlation analysis between the abundance of these altered immune cells and key DEGs (Tnf, Il1b, Pparg, Ccl5, Cxcl1) showed a clear dichotomy. The increased M2 macrophage signature positively correlated with Il1b, Pparg, and Cxcl1, but negatively correlated with Tnf and Ccl5. The opposite correlation pattern was observed for the diminished NKT and Germinal center B cell populations (Fig. 6C). Together, these data indicate that FGF21 deficiency dramatically alters the retinal immune milieu, primarily by disrupting macrophage homeostasis.