Transcriptome Analysis of Choroid and Retina From Tree Shrew With Choroidal Neovascularization Reveals Key Signaling Moieties

Pathological neovascularization in choroid, a leading cause of blindness, is a characteristic of many fundus diseases, such as diabetic retinopathy and age-related macular degeneration. The present study aimed to elucidate the key signaling pathways in choroidal neovascularization (CNV) by analyzing the mRNA profiles of choroid and retina in tree shrews with CNV. We induced choroidal angiogenesis by laser photocoagulation in 15 tree shrews and obtained mRNA profiles of their choroids and retinas by high-throughput transcriptome sequencing. Hierarchical cluster analysis, weighted gene co-expression network analysis (WGCNA), protein-protein interaction (PPI) network analysis, hematoxylin and eosin (HE) staining, CD31 immunohistochemistry (IHC), and reverse transcription quantitative PCR (RT-qPCR) were performed. After laser photocoagulation, we obtained a total of 350 differentially expressed genes (DEGs) in the choroid, including 59 genes in Module-FASN (“ME-FASN”) module and 28 genes in Module-RPL (“ME-RPL”) module. A total of 69 DEGs in retina, including 20 genes in Module-SLC (“ME-SLC”) module. Bioinformatics analysis demonstrated that DEGs in choroid were mainly involved in membrane transport; DEGs in “ME-RPL” were prominent in pathways associated with IgA production, antigen presentation, and cell adhesion molecules (CAMs) signaling. DEGs in “ME-FASN” were involved in fatty acid metabolism and PPAR signaling pathway, while DEGs in “ME-SLC” were involved in GABAergic synapse, neuroactive life receptor interaction, cholinergic synapse, and retrograde endocannabinoid signaling pathway. PPI network analysis demonstrated that the ribosomal protein family genes (RPL31, RPL7, RPL26L1, and RPL19) are key factors of “ME-RPL,” acyl-CoA superfamily genes (ACACA, ACAT1, ACAA2, and ACACB) and FASN are key factors of “ME-FASN” and superfamily of solid carrier genes (SLC17A6, SLC32A1, SLC12A5, and SLC6A1) and complement genes (C4A, C3, and C2) are key factors of “ME-SLC.” In conclusion, the present study discovered the important signal transductions (fatty acid metabolic pathway and CAMs signaling) and genes (ribosomal protein family and the complement system) in tree shrew CNV. We consider that our findings hold implications in unraveling molecular mechanisms that underlie occurrence and development of CNV.


INTRODUCTION
Choroidal neovascularization (CNV) is the formation of new blood vessels in the choroid. The blood vessels form/develop between the retinal pigment epithelium (RPE) and Bruch's membrane. They extend through the retinal neuroepithelial layer eventually forming a fibrous vascular tissue. CNV often occurs in the macular area, thereby causing macular hemorrhage and serous exudation under the retina. CNV is a common characteristic of many fundus diseases, such as age-related macular degeneration, high myopic maculopathy, central exudative chorioretinopathy, and diabetic retinopathy; it renders grave vision impairment in the affected individuals.
Although the pathogenesis of CNV remains poorly understood, it is believed to involve a variety of cell growth factors, such as the vascular endothelial growth factor (VEGF), basic fibroblast growth factor (bFGF), and platelet-derived growth factor (PDGF), available in the local microenvironment (Ming et al., 2011). Moreover, inflammatory cells and cytokines are also involved in its pathogenesis, while tumor necrosis factor-α (TNF-α), interleukins 6 (IL6), and intercellular cell adhesion molecules 1 (ICAM-1) released by macrophages and neutrophils promote the occurrence and development of CNV (Jin et al., 2010).
The early stage of CNV is characterized by changes in the retinal microenvironment along with production of VEGF by the RPE and photoreceptors (PRs; Nagineni et al., 2012). Further, VEGF promotes the migration of macrophages to the Bruch's membrane, eventually leading to proteolytic degradation of the membrane (Grossniklaus et al., 2002). However, the key molecular mechanisms underlying CNV and its signaling between retina and choroid remain unclear.
Animal models are important tools in investigating pathogenesis of CNV. The most commonly used method to establish a CNV model is laser-induced selective destruction of photoreceptors, RPE cells, Bruch's membrane, and choroidal capillaries. Destruction of the Bruch's membrane stimulates a series of damage repair processes that leads to the formation of new blood vessels (Archer and Gardiner, 1981;ElDirini et al., 1991;Yang et al., 2006;Hoerster et al., 2012;Liu et al., 2013;Kunbei et al., 2014;Poor et al., 2014). Pathogenesis of CNV in humans is same as that in animals, wherein an injury to the Bruch's membrane-RPE-choroidal complex leads to an imbalance between angiogenic and inhibitory factors, resulting in a series of neovascularization processes, such as endothelial cell proliferation and migration, and lumen formation.
Animals that are most often used for CNV experiments are monkeys, rabbits, and mice. Although the monkey model manifests changes/characteristics similar to that by humans, their application in CNV studies is limited owing to their high cost of maintenance and low incidence of CNV (Archer and Gardiner, 1981;Kunbei et al., 2014). Further, the CNV rabbit model fails to exhibit typical leakage characteristics (less than 30%) as observed by FFA. In addition, retinal vascular circulation of mice and rabbits is disparate from that of humans; therefore, mice and rabbits are not ideal models for CNV studies (Zhu et al., 1989;ElDirini et al., 1991;Tamai et al., 2002;Yang et al., 2006;Hoerster et al., 2012;Liu et al., 2013;Poor et al., 2014). Tree shrew is a small mammal sharing evolutionary and anatomical similarities with primates. Tree shrews have a well-developed visual system, with cone cells accounting for 96% of all the photosensitive cells. These mammals have good color vision as well as stereoscopic vision. In recent years, tree shrews have been used for investigating visual development and ophthalmic diseases (Lin et al., 2013;Jia and Dai, 2019). The retina of tree shrews is mainly composed of cone cells, as mentioned above, and exhibits functions similar to that of human retina. The unique retinal structure of tree shrew provides a good basis for investigating the disease process of human CNV.
In the present study, we aimed to investigate the signal transduction in choroid and retina in a tree shrew CNV model. We found that the genes of ribosomal protein family, superfamily of acyl-CoA, and solute carrier family were differentially expressed. Our findings imply that fatty acid metabolic pathway and CAM pathway play key roles in tree shrew CNV.

Experimental Animals
Fifteen adult tree shrews (Tupaia belangeri chinensis) with healthy eyes (seven females and eight males, aged 2-3 years, weighing 110-130 g) were collected from the Tree Shrew Germplasm Resource Center, Institute of Medical Biology, Chinese Academy of Medical Sciences, Kunming, China. The experimental animal production licenses were SCXK (Dian) K2018-0002 and use license SYXK (Dian) K2018-0002. All animal experiments were approved by the Animal Ethics Committee of the Institute of Medical Biology, Chinese Academy of Medical Science (approval number: DWSP201803019). The tree shrews were randomly divided into three groups (7, 21, and 30 days after laser photocoagulation) such that each group comprised five tree shrews (10 eyes per group, all left eyes were used as the experimental group and all right eyes as the control group).

Laser-Induced CNV Model
The animals were anesthetized by intraperitoneal injection of 0.3% pentobarbital sodium (0.6 mg/kg), and the body temperature was maintained by a heating pad. Five minutes before laser photocoagulation, compound tropicamide eye drops were used to disperse the pupil of both eyes. Eyes were anesthetized with two drops of 4 mg/ml oxybuprocaine hydrochloride eye drops (Santen Pharmaceutical Co., Ltd., China) and carbomer eye drops (Dr. Gerhard Mann, Chem.-Pharm. Fabrik GmbH, Germany) was used to prevent corneal dryness. After the animals were fixed, laser photocoagulation (Lumenis, United States) was performed at 1 PD around the optic disc, 15 spots in total. The laser wavelength was 647.1 nm, diameter of the spots was 50 μm, and the exposure time was 0.01-0.05 s. Energy of the laser was adjusted according to the retinal reaction. The effective spots were marked by the formation of bubbles without bleeding after photocoagulation. Eyes with ruptured capillaries small blood vessels were not used for subsequent experiments (Fundus camera: TOPCON, United States).

Hematoxylin and Eosin Staining
Eyeballs of tree shrew corneas were fixed with FAS (Servicebio, China) for 24 h, and 3-um-thick of were cut. The eyeball sections were dehydrated, cleared, embedded, and made into sections. The sections were used for hematoxylin and eosin (HE) staining and CD31 immunohistochemistry (IHC). The pathological changes were described by a section scanner (Mantra, United States).

CD31 Immunohistochemistry
In order to confirm the occurrence of CNV after photocoagulation, IHC for CD31 was conducted. Antigen repair was performed in citric acid sodium buffer solution (ph6.0; Servicebio, China). Sections were incubated in 3% hydrogen peroxide at room temperature for 25 min in dark to block endogenous peroxidase, and then blocked with 3% BSA.
Sections were incubated in 1:300 dilution of Anti-CD31 Rabbit pAb (Servicebio, China) at 4°C overnight, and then incubated in 1:200 dilution of HRP conjugated Goat Anti-Rabbit IgG (H + L; Servicebio, China) for 50 min at room temperature. Color development was performed using DAB (Servicebio, China). The result of CD31 IHC was observed under microscope (Mantra, United States), and inform was used for processing the image. Microvessel density (MVD) were counted refer to Weidner et al. (1993).

RNA Sequencing of Retina and Choroid of Tree Shrew CNV Model
Total RNA Extraction Tree Shrews were euthanized by injection of 2% pentobarbital sodium (2 mg/kg). Tree shrews were disinfected and the eyeballs were rinsed with iodophor. The canthus and the conjunctiva were cut. The conjunctiva and fascia were separated so that the sclera was fully exposed. The optic nerve was cut and then the eyeball was removed.
The eyeballs were suspended in PBS in clean bench. After the cornea was cut, the lens and vitreum were separated. The retina-choroid-sclera was transferred to RNase-free water (TaKaRa, Japan) to isolate the retina (soft and transparent). Since the choroid was fragile and easy to fall off, the choroidsclera was transferred to TRIzol and choroid (brown) was isolated.
Total RNA was extracted from the collected samples using TRIzol according to the standard protocol (Invitrogen, United States). Purity of the extracted RNA was spectrophotometrically quantified with NanoPhotometer (Implen, United States). RNA integrity was determined using Bioanalyzer 2100 system (Agilent Technologies, United States). In case the RNA content was insufficient (RNA < 2 μg), the sample was mixed with other samples of the same group.

cDNA Library Preparation and Sequencing
First-strand DNA was synthesized using random primers and M-MuLV reverse transcriptase (RNase H-). Second-strand DNA was synthesized by DNA polymerase I and RNase H, wherein dTTPs were replaced by dUTPs. Further, the dU-containing second-strand cDNA was degraded by the USER enzyme, and the cDNA was PCR-amplified to obtain a library.
The cDNA was quantified using Qubit 2.0, and insert size of the library was determined by diluting the library to 1 ng/ μl and subjecting to Agilent 2100 system. Effective concentration of the cDNA library was quantified accurately by qPCR (effective library concentration > 2 nM) to ensure quality of the cDNA library. Twenty-six DNA libraries were constructed. The libraries were sequenced on an Illumina HiSeq 4000 platform at the Novogene Bioinformatics Institute (Beijing, China).

GO and KEGG Analyses
In order to determine the functions of the DEGs, clusterProfiler was used to analyze gene ontoloagy (GO) function and kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment. GO analyzes the biological processes, cellular components, and molecular functions. KEGG analysis determines the signal pathways enriched by the DEGs, and it is expected to find the biological pathways that play key roles in tree shrew CNV. Threshold of GO/KEGG enrichment was considered significant at Padj < 0.05.

Protein-Protein Interaction Network Construction
STRING database 1 provides comprehensive information about interactions between proteins. In the present study, STRING was used to generate PPI networks among the differentially expressed mRNAs based on interactions with combined scores > 0.4. Additionally, Cytoscape was used to visualize the network, and PPI scores > 800 were considered for functional enrichment analysis of the modules. Weighted Gene Co-expression Network Analysis Hierarchical cluster tree was established using WGCNA based on the correlation of gene expressions at the Pearson correlation coefficient 0.8 and soft threshold (power) 9. Functions were considered relevant when the genes demonstrated a high degree of co-expression correlation in a module. The modules were clustered and the correlation heatmap among modules was constructed to evaluate the connectivity between two genes within the module according to the module eigenvalues. We analyzed the expression trends of DEGs in the modules during different courses of CNV.
The reaction was performed at the following conditions: initial cDNA synthesis at 45°C for 5 min; initial denaturation at 95°C for 10 s; 40 cycles of denaturation at 95°C for 5 s, annealing and extension at 60-61°C for 30 s; 95°C for 15 s, 60°C for 60 s, and 95°C for 15 s for the dissolution curve. The internal reference gene was GAPDH. Each experiment was performed in triplicate. PCR primers and amplification conditions used in this study are shown in Table 1.

Statistical Analysis
Statistical analyses were performed using SPSS 20.0. All data are expressed as mean ± SEM. p < 0.05 was considered statistically significant. Student's t-test was used for analyzing RT-qPCR and MVD results. GraphPad Prism 8.0 was used to prepare the graphs.

Pathology of the Tree Shrews CNV
In the control groups, tree shrews have intact retina ( Figure 1L). After 7 days of photocoagulation, retinal edema (Figures 1B,F), damaged haller's layer in choroid (Figures 1A,E) was observed. After 14 days of photocoagulation, neovascularization in retina leads to destruction of ganglion cell (GCL) layer and disorder of nerve fiber layer (NFL) layer. CNV was developed, and neovascularization surrounded by neutrophils, macrophages, and fibroblasts was observed found in choroid (Figures 1C,G). Haller's layer and sattler's layer was loose in choroid (Figures 1D,H). After 30 days of photocoagulation, retinal neovascularization was developed (Figure 1I). Neutrophils and macrophages were seen around the neovascularization in choroid. Damaged choroidal structure was not recover (Figures 1J,K).

IHC for CD31
The normal tree shrews have complete retina, choroid and sclera with distinct layers, and no neovascularization was observed (Figure 2A). After 7 days of photocoagulation, retinal edema was observed in the photocoagulation sites, and the choroid structure was destroyed. Neovascularization was observed in the subchoroid ( Figure 2C) and scleral ( Figure 2B), but no neovascularization was observed in the choroid and retina. After 14 days of photocoagulation, retinal edema was still observed in the photocoagulation sites. Neovascularization with small size was obvious in the retina (Figures 2D,F) and choroid (Figure 2D), while the neovascularization in the sclera was larger than that in 7 days of photocoagulation group (Figure 2E). After 21 days of photocoagulation, increased retinal ( Figure 2I) and choroid ( Figure 2H) neovascularization in the photocoagulation sites. The lumen of retinal neovascularization was larger than that in 7 days of photocoagulation group ( Figure 2G).
Microvessel density in laser photocoagulation group was higher than that in control group (p < 0.05). MVD in 14 days of photocoagulation group (7.00 ± 1.00) were higher than that in 7 days of photocoagulation group (3.33 ± 0.58). However, MVD were no significant difference between 14 days of photocoagulation group and 21 days of photocoagulation group (5.33 ± 2.31).

RNA Sequencing
Raw reads were generated from 36 samples, including 18 choroids and 18 retinas. The average clean data of each sample were not less than 7 GB. The average number of clean reads was 46,719,819; Q20 was higher than 97%; Q30 was higher than 94%.
The expression of ACAT1 was upregulated after 7 days (p = 0.039) and 21 days (p = 0.039) of photocoagulation in the choroid (Figure 3E), but downregulated after 21 days (p = 0.001) of photocoagulation in the retina ( Figure 3J). The qPCR results were consistent with the RNA sequencing data.

Cluster Analysis of DEGs
The mRNA expression levels were estimated with FPKM. In the choroid, 335 mRNAs, 9 mRNAs, and 6 mRNAs were differentially expressed after 7, 21, and 30 days of laser photocoagulation, respectively. In the retina, 6, 56, and 7 mRNAs were differentially expressed after 7, 21, and 30 days of photocoagulation, respectively. A cluster analysis of the differentially expressed mRNAs was conducted and the results are shown as heatmaps in Figure 4.

Heatmap of Inter-Module Correlation
Correlation analysis between the modules and samples showed that "ME-RPL, " "ME-FASN, " and choroid are highly (and positively) correlated, while "ME-SLC" and retina are highly (and positively) correlated, as depicted in Figure 6.

Cluster Heatmap of the Gene Modules
Analysis of interactions among the gene modules revealed that when the color of the region within and between the three modules is evident, the genes in these modules have a high degree of association. It is suggested that "ME-RPL, " "ME-SLC, " and "ME-FASN" modules interact closely with each other. The results in the form of heatmap visualization are shown in Figure 7.

Expression Patterns of the Gene Modules
In order to analyze the expression trend of the genes in different modules, the module gene expression pattern was constructed. The results showed that "ME-RPL" genes were downregulated in the retina but upregulated in choroid; "ME-SLC" genes were upregulated in the retina but downregulated in choroid; and "ME-FASN" genes were downregulated in retina but upregulated in choroid, as shown in Figure 8.

GO and KEGG Analyses of DEGs
GO analysis showed that the differentially expressed mRNAs in the choroid were mainly involved in membrane transport processes, while there was no significant enrichment for differentially expressed mRNA in the retina (Figure 9). KEGG analysis showed that "ME-RPL" genes were enriched in IgA production, antigen presentation, and CAMs pathways. "ME-FASN" genes are enriched in fatty acid metabolism and PPAR signaling pathway, while "ME-SLC" genes were enriched in GABAergic synapse formation/development, neuroactive ligand-receptor interaction, cholinergic synapse, and retrograde endocannabinoid signaling pathway (Figure 10).
In the present study, we showed that the expression of IL18 increased significantly after 21 days of laser photocoagulation in the choroid, although it decreased significantly after 7 and 21 days of laser photocoagulation in the retina. In mouse model of CNV, the expression of IL18 decreased significantly in RPE-choroid tissue (Choudhary et al., 2015). It has been reported that IL18 regulates pathological retinal neovascularization. In addition, IL18 inhibits the formation of experimental CNV through the stimulation of interferon-γ and thrombospondin (Xiao and Wu, 2016). IL18 reduce CNV development in the nonhuman primate eye, and inhibits vascular leakage in a mouse model of spontaneous neovascularization (Doyle et al., 2015). However, the mechanism of IL18 inhibiting CNV formation remains unclear and needs to be further investigated.
In present study, genes in "ME-RPL" had a high correlation with choroid. Moreover, KEGG analysis showed that "ME-RPL" genes were significantly enriched in CAMs pathway. CAMs are associated with neovascularization (Dong et al., 2018), wherein ICAM-1 is closely related to the occurrence and development of CNV (Gao and He, 2015). However, the role of CAMs in CNV remains to be fully understood.
There were significant differences in the expression of ribosomal protein family genes in the "ME-RPL" module, and the interaction among proteins in "ME-RPL" was close. It has been reported that RPL17 inhibits the growth of vascular smooth muscle cells and silencing RPL17 promotes proliferation of these cells (Smolock et al., 2012). It is suggested that the ribosomal protein family may be involved in the regulation of vascular growth. Concurrent with previous studies, we found that the ribosomal protein family may be involved in the regulation of CNV. However, the underlying mechanism requires further in-depth research.
In the "ME-FASN" module, fatty acid synthase (FASN) was the center of PPI network. KEGG analysis showed that many "ME-FASN" DEGs were enriched in the fatty acid metabolic, GABAergic synapse, neuroactive life receptor interaction, cholinergic synapse, and retrograde endocannabinoid pathways. FASN regulates tumor angiogenesis by altering secretion and activity of VEGF (Seguin et al., 2012). This study suggests that the acylation involved in fatty acid metabolism may be related to CNV. The acyl-CoA superfamily genes (ACACA, ACACB, ACAA2, and ACADVL) in "ME-FASN" module were significantly differentially expressed and had strong interactions with FASN. In mouse model of CNV, the ACACB and ACADVL gene was significantly differentially expressed (Choudhary et al., 2015), which correlate with present study. The present study indicates that fatty acid metabolic pathway and acyl-CoA superfamily genes may play an important role in the occurrence and development of CNV. It is reported that FASN knockdown in endothelial cells elevated malonyl-CoA levels, and reduced pathological ocular neovascularization (Bruning et al., 2018). The expression levels of acyl-CoA oxidase 1 (Acox1), fatty acid binding protein 4 (Fabp4) is associated with the progression of retinal pathological neovascularization in a murine model (Tomita et al., 2019). SLC27A2 and C3 gene in "ME-SLC" were significantly differentially expressed in CNV tree shrew models and CNV mouse model (Choudhary et al., 2015). The abscissa is the sample and the ordinate is the module. The numbers represent correlation between the module and the sample. Closer the value is to 1, stronger the positive correlation between the module and the sample; closer to −1, stronger the negative correlation. The number in brackets represents value of p, smaller the value, and higher the significance. The red specimen is choroid of the tree shrews.
In conclusion, we constructed the interaction network of "ME-RPL" involved in CAMs pathway and "ME-FASN" genes involved in fatty acid metabolism in choroid, "ME-SLC" genes involved in GABAergic synapse in retina during CNV in tree shrews. Our findings hold implications in unraveling molecular mechanisms that underlie occurrence and FIGURE 7 | Network heatmap of the gene modules. The tree represents a module (top and left), and the branch represents a gene. Darker the color of the dot (white → yellow → red), stronger the connectivity between the two genes corresponding to the row and column. Yellow tree represents "ME-FASN" module; Blue tree represents "ME-RPL" module; Turquoise tree represents "ME-SLC" module. FIGURE 9 | GO enrichment of the significant differentially expressed genes (DEGs) in the choroid and retina. The abscissa shows the rich factors and the ordinate represents the pathways. Size of each point indicates the number of candidate target genes in the GO term, and the color of each point corresponds to the value of p after correction.
FIGURE 10 | KEGG enrichment of the significant differentially expressed mRNAs in "ME-RPL," "ME-FASN," and "ME-SLC." The abscissa shows the rich factors and the ordinate represents the pathways. Size of each point indicates the number of candidate target genes in the pathway, and the color of each point corresponds to the value of p after correction.
FIGURE 12 | PPI networks of the DEGs in "ME-RPL" and "ME-SLC." Grayscale value of each line correlates with the strength of the interaction between the proteins. development of CNV. Although with cone cells accounting for 96% of all the photosensitive cells, tree shrew not have discernible macula. The problem of lack established transgenic technology for tree shrews need to be solved, and the role of CAMs and fatty acid metabolic pathway needs further investigations to elucidate CNV pathogenesis at the molecular level.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: The European Molecular Biology Laboratory's European Bioinformatics Institute (EMBL-EBI) ArrayExpress, https://www.ebi.ac.uk/arrayexpress/, E-MTAB-10198.