key: cord-0312335-z83ik5rv
authors: Wang, J.-H.; Wong, R. C. B.; Liu, G.-S.
title: Retinal aging transcriptome and cellular landscape in association with the progression of age-related macular degeneration
date: 2022-04-04
journal: nan
DOI: 10.1101/2022.04.03.22273375
sha: 8069c582a3932dcf4be7e9a8c692678f3ceb66c2
doc_id: 312335
cord_uid: z83ik5rv

Age is the main risk factor for age-related macular degeneration (AMD), a leading cause of blindness in the elderly, with limited therapeutic options. Here we systematically analyzed the transcriptomic characteristics and cellular landscape of the aging retina from controls and patients with AMD. We identify the aging genes in the retina that are associated with innate immune response and inflammation. Deconvolution analysis reveals that the estimated proportion of M2 and M0 macrophages is increased and decreased, respectively with both age and AMD severity. Moreover, we find that Muller glia are increased with age but not with disease severity. Several genes associated with both age and disease severity in AMD, particularly C1s and MR1, are strong positively correlated with the proportions of Muller glia. Our studies expand the genetic and cellular landscape of AMD and provide avenues for further studies on the relationship between age and AMD.

CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. and M2 macrophages were significantly increased with MGS level, while M0 macrophages 116 were also decreased (Extended Data Figs. 4 and 5, Supplementary Data 4). It is known that 117 human peripheral monocytes were differentiated into inactivated M0 macrophages and then 118 polarized to M1 and M2 phenotypes under various stress or stimuli 16 

shown that numbers of both M1 and M2 macrophages were increased in the retina 12 or 120 aqueous humors 17 of a small number of neovascular AMD patients or laser-induced murine 121 model with choroidal neovascularization 17, 18 . In line with these reports, our results further 122 revealed that age alone as a predominant factor that triggers the polarization of M0 123 macrophages towards to pro-inflammatory M1 and pro-angiogenic M2 phenotypes 12 , while 124 disease severity alone seems majorly affecting the numbers of M2 macrophages that feature 125 in the pro-angiogenic process. Collectively, these results together showed that activated M2 126 macrophages caused by age and disease severity is the major immune cell type involved in 127 the immune response in the pathogenesis of AMD. 128 129

We next investigated whether retinal aging is related to MGS level in AMD patients. Using 132 the OLR model that controls sex and age, we correlated age-associated genes to MGS level, 133 Although none of the Age-MGS genes we identified has been previously reported in 148 genome association studies in AMD, many plays critical roles in immune responses and 149 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint 6 inflammation in the retina. For example, S100 calcium binding protein B (S100B), a well-150 known biomarker of active neural distress that is involved in many neurodegenerative 151 diseases 19 , are pro-inflammation and can induce retinal degeneration in animal models 20, 21, 22 . We found that MATN1 expression was decreased with the disease severity. This result is 169 consistent with a study that showed increased angiogenesis during fracture healing in a 170 MATN1 deficient mouse, and purified MATN1 proteins enabled inhibition of angiogenesis 171 both in vitro and in vivo 28 . Of note, VSTM4, FBLN5, MPZL3 and SCPEP1, have been 172 reported to promote sprouting angiogenesis by endothelial cells 29, 30, 31, 32 , while their role in 173 retinal angiogenesis remains unstudied. 174

We also found that a group of Age-MGS genes are related to the regulation of retinoic 175 acid metabolic process (Supplementary Data 6), including cytochrome P450 family 26 176 subfamily A member 1 (CYP26A1) or B member 1 (CYP26B1), dehydrogenase/reductase 3 177 (DHRS3) and serine carboxypeptidase 1 (SCPEP1), and their expression (except for SCPEP1) 178

were reduced with disease severity. Retinoic acid is an active metabolite of vitamin A 179 (retinol), an essential signaling molecule involved in the normal development of the ventral 180 retina and optic nerve 33 as well as maintenance of the blood-retinal barrier 34 . Dysregulation of 181 retinoic acid can cause degenerative retinal diseases. CYP26A1, a retinoic acid catabolizing 182 enzyme, was reported to be enriched in the fovea in the human retinas and functions for 183 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint photoreception in rod-free zone generation 35 . Reduction of CYP26A1 may affect degradation 184 of retinoic acid, ultimately causing dysfunctionality of photoreceptors. Similarly, CYP26B1 is 185 also responsible for retinoic acid metabolism 36 , however, its role in the retina has not been 186 investigated yet. DHRS3, a retinal reductase and a possible gene target of all-trans-retinoic 187 acid, has been reported to be reduced in the liver during inflammation likely resulting in the 188 perturbation of the whole-body vitamin A metabolism that occurred in conditions with 189 inflammatory stress 37 . Together, these results imply that reduction of genes associated with CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint

We next sought to understand the correlation between the expression of Age-MGS genes. 219 Correlation analysis revealed 3 clear clusters of Age-MGS genes with correlated expression 220 using the EyeGEx database with sex and age controlled. (Fig. 3a) . For ease of description, 221 these clusters were numbered as cluster A, B and C. Cluster A consists of zinc finger protein To assess the cell type-specific expression pattern of the age-associated genes and 248

Age-MGS genes, we further analyzed the neural retina single-cell RNA-seq (scRNA-seq) 249 data from Human Cell Atlas (HCA, Extended Data Fig. 7 ) 52 . We found that the expression 250 pattern of the age-associated genes is distinctive among 9 retinal cell types, and a higher 251 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. cell types in the condition of AMD is not clear. We also found that correlation between 265 AMD-associated genes identified through different methodologies, such as LRT combined 266 with OLR (present study), transcriptome-wide association study (TWAS using EyeGEx) 7 and 267 GWAS 4 , is clearly distinct (Extended Data Fig. 9a) , indicating that these is limited interaction 268 between these genes. However, genes identified by TWAS (9 genes) and GWAS (25 leading 

To determine the functional roles of these Age-MGS genes involved in the pathogenesis of 277 AMD, we applied functional enrichment analysis (Biological Process) by Metascape. The 278 results revealed significant enrichment most related to regulation of innate immune response 279 (Fig. 4a, Supplementary Data 7) . Among which, innate immune response, leukocyte 280 activation, and multicellular organismal homeostasis are centered in the cluster of GO terms 281 and have most interaction (n ≥5) with other GO terms (Fig. 4b) , suggesting their important 282 roles in regulation of Age-MGS genes. Therefore, we combined genes involved in these 283 biological processes (Supplementary Data 7) and focused on the 8 Age-MGS genes with 284 significant changes with both MGS levels (Fig 4c) and age (Fig 4d) ), including C1s, C7, MR1, 285 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. the scRNA-seq data from HCA as a reference (Extended Data Fig. 7, Supplementary Data 8) . 313

To increase the confidence of down-stream analyses, the cell types were retained if their 314 estimated proportions were >0 in all samples (Extended Data Fig. 10) . 315

To assess whether age as a risk factor affects cell type proportions, we used an OLR 316 model to control for sex and MGS level and identified age-associated alteration in the 317 proportions of retinal cell types. We observed that the proportions of Müller glia were 318 significantly increased with age, while the proportions of horizontal and rod cells were 319 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint significantly decreased with age (Fig. 5a, Extended Data Fig. 11, Supplementary Data 9) . 320

Next, we focused on the risk factor of MGS level to assess its effect on cell proportions by 321 controlling for sex and age using an OLR model. Although the proportions of Müller glia and 322 rod were increased and decreased, respectively, with MGS levels, the changes were not 323 significant (Fig. 5b, Extended Data Fig. 12, Supplementary Data 10 is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint relationship in this form of retinal degeneration. Indeed, our results revealed that there is a 354 strong negative correlation (r = -0.89) between rod photoreceptors and Müller glia based on 355 its estimated proportions after controlling for sex and age (Fig. 5d) . In the present study, we systematically analyzed the retinal aging transcriptome and cellular 360 landscape in association with the progression of age-related macular degeneration. We 361 revealed that the aging retina is characterized by a range of biological processes, in particular 362 innate immune response and inflammation that could lead to the progressively severe AMD. Future research to investigate the function of these Age-MGS genes in the retina would be 385 important to understand the mechanism underly AMD pathology. 386 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint

Since the transcriptome data from human donors are complexed with a variety of 387 latent factors, such as sex, age and clinical diagnosis, it is important to use batch correction to 388 adjust those factors for the downstream analysis. Throughout our analysis, all the raw data 389 were subjected to the batch correction before relevant analysis. We noted a study using 390 various deconvolution algorithms to analyze involvement of retinal cell types in AMD 391 patients with different databases, including EyeGEx 67 . However, all the raw data from this 392 study was not subjected to batch correction, thus the conclusion resulted from was less 393 stringent compared to the present study. A clear example of difference between adjusted and 394 unadjusted data analysis can be referred to the cell type proportion estimated by 395

CIBERSORTx (Extended Data Fig. 3, 5, 11 and 12 ). There were several limitations of this To identify age-MGS associated genes, the TPM values of age-associated genes from both 445 clusters (increasing or decreasing with age) were analyzed by R package MASS (v7.3-54) 446 and ordinal (v2019.12-10) using ordinal logistic regression (OLR) model 71 . Using 447

Cumulative Link Model function in ordinal 71 , we controlled sex and age to screen input genes 448 that are associated with MGS level. Age-MGS-associated genes were determined at a 449 significance threshold of p < 0.05. The resulting genes were represented in a forest plot and 450 ranked according to its beta coefficient. 451

For gene-gene correlation of Age-MGS genes, the raw counts were processed by 452

Variance Stabilizing Transformation function in DESeq2 8 , followed by statistical adjustment 453 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

The authors declare no competing interests. 542 543 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint

Risk factors for age-related macular degeneration

Risk factors and biomarkers of age-related macular degeneration

Age-related macular degeneration is the leading cause of blindness

Fibulin-5 antagonizes vascular endothelial growth factor 608 (VEGF) signaling and angiogenic sprouting by endothelial cells

Transcriptomics of post-stroke angiogenesis in the aged brain

Newly Identified Peptide, Peptide Lv

Mesenchymal Transition (Plasticity) in CIN and

Early Invasive Carcinoma of the Cervix: Exploring Putative Molecular Mechanisms 618 Involved in Early Tumor Invasion

Retinoic acid signaling in mammalian eye development

Retinoic acid signaling is essential for 622 maintenance of the blood-retinal barrier

Fgf8 Expression and Degradation of Retinoic Acid Are 624

Required for Patterning a High-Acuity Area in the Retina

Photopic light-640 mediated down-regulation of local alpha1A-adrenergic signaling protects blood-retina 641 barrier in experimental autoimmune uveoretinitis

The retina as a window to the brain-from eye 643 research to CNS disorders

Intermediate filament proteins and their 646 associated diseases

Differential expression of GFAP in 648 early v late AMD: a quantitative analysis

Complement C1s and C4d as Prognostic Biomarkers in Renal 650

Cancer: Emergence of Noncanonical Functions of C1s

Systems-level analysis of age-related macular degeneration 653 reveals global biomarkers and phenotype-specific functional networks

Retinal Macrophages Synthesize C3 and Activate Complement in 656 AMD and in Models of Focal Retinal Degeneration

MR1-restricted mucosal associated invariant T 659 (MAIT) cells in the immune response to Mycobacterium tuberculosis

Early Involvement of

Immune/Inflammatory Response Genes in Retinal Degeneration in DBA/2J Mice

SLC9A9 Co-expression 665 modules in autism-associated brain regions

A single-cell transcriptome atlas of the adult human retina

Single-cell transcriptomic atlas of the human retina identifies cell 669 types associated with age-related macular degeneration

Glia-neuron interactions 671 in the mammalian retina

International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

New functions of Muller cells

Aging Changes in Retinal Microglia and their Relevance to Age-675 related Retinal Disease

The pivotal role of the complement system in aging and age-677 related macular degeneration: hypothesis re-visited

Evidence for association between multiple complement 680 pathway genes and AMD

Molecular pathology of age-related macular degeneration

Pathway Analysis Integrating Genome-Wide and Functional 684

Data Identifies PLCG2 as a Candidate Gene for Age-Related Macular Degeneration

Photoreceptor loss in age-related macular 690 degeneration

Photoreceptor topography in ageing and age-related maculopathy

Ablation of retinal horizontal cells from adult mice leads to rod 694 degeneration and remodeling in the outer retina

A single-cell transcriptome atlas of the aging human and macaque retina

Drusen-associated degeneration in the retina

Implication of specific retinal cell-type involvement and gene expression 700 changes in AMD progression using integrative analysis of single-cell and bulk RNA-701 seq profiling

Relationship between 703 RPE and choriocapillaris in age-related macular degeneration

International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

Comprehensive Integration of Single-Cell Data

Impact of outdated gene 708 annotations on pathway enrichment analysis

Modeling continuous response variables 710 using ordinal regression

limma powers differential expression analyses for RNA-sequencing 712 and microarray studies

Visualization of a 714 correlation matrix. R package version 073 230

The aging transcriptome and cellular landscape of the 716 human lung in relation to SARS-CoV-2

the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

for sex and age using the Remove Batch Effect function in limma 72 

For visualization of adjusted gene expression levels, the raw counts were processed by 458Variance Stabilizing Transformation function in DESeq2 8 , followed by statistical adjustment 459 for sex and MGS level or sex and age using the Remove Batch Effect determined as per the cell type annotation resulting from above analyses. We further scaled 483 the expression frequencies in cell types in z-score using the NMF package (v0.23.0). 484 485

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint Interaction network of gene ontology annotations described in Fig. 4a . Annotations 757 highlighted in red indicates that the specific biological process has the most interactions with 758 other biological processes. c, Tukey boxplots (interquartile range (IQR) boxes with 1.5× IQR 759 whiskers) showing the expression of C1s, C7, MR1, PLCG2, TRIM22, TMEM98, CYP26B1,  760 and VSTM4. Gene expression value are shown as log-transformed, controlled for sex and age 761 (Fig. 4c) or sex and MGS level (Fig. 4d) . Statistical significance of age-associated or MGS-762 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 4, 2022. ; https://doi.org/10.1101/2022.04.03.22273375 doi: medRxiv preprint