Molecular Vision 2020; 26:459-471
Received 19 November 2019 | Accepted 17 June 2020 | Published 19 June 2020
The first two authors contributed equally to this work.
1Medicine and Pharmacy Research Center, Binzhou Medical University, Yantai, Shandong, China; 2Department of Genetics, Genomics and informatics, University of Tennessee Health Science Center, Memphis, TN; 3State Key Laboratory of Ophthalmology, Zhongshan Ophthalmic Center, Sun Yat-Sen University, Guangzhou, Guangdong, China; 4Analytical Chemistry and Neurochemistry, Department of Chemistry-BMC, Uppsala University, Uppsala, Sweden; 5Transformative Otology and Neuroscience Center, Case Western Reserve University School of Medicine, Cleveland, OH; 6Departments of Otolaryngology, Case Western Reserve University School of Medicine, Cleveland, OH.
Correspondence to: Geng Tian, Medicine and Pharmacy Research Center, Binzhou Medical University, Yantai, Shandong, China; Phone: +86-535-6913395; FAX: +86-535-6913034; email: firstname.lastname@example.org
Purpose: Platelet-derived growth factor (PDGF) signaling is well known to be involved in vascular retinopathies. Among the PDGF family, the subunit B (PDGFB) protein is considered a promising therapeutic target. This study aimed to identify the genes and potential pathways through which PDGFB affects retinal phenotypes by using a systems genetics approach.
Methods: Gene expression data had been previously generated in a laboratory for the retinas of 75 C57BL/6J(B6) X DBA/2J (BXD) recombinant inbred (RI) strains. Using this data, the genetic correlation method was used to identify genes correlated to Pdgfb. A correlation between intraocular pressure (IOP) and Pdgfb was calculated based on the Pearson correlation coefficient. A gene set enrichment analysis and the STRING database were used to evaluate gene function and to construct protein–protein interaction (PPI) networks.
Results: Pdgfb was a cis-regulated gene in the retina; its expression had a significant correlation with IOP (r = 0.305; p value = 0.012). The expression levels of 2,477 genes also had significant correlations with Pdgfb expressions (p<0.05), among which Atf4 was the most positively correlated (r = 0.628; p value = 1.29e-10). Thus, Atf4 was highly expressed in the retina and shared the transcription factor (TF) Hnf4a binding site with Pdgfb. Gene Ontology and a pathway analysis revealed that Pdgfb and its covariates were highly involved in mitogen-activated protein kinase (MAPK) and vascular endothelial growth factor (VEGF) pathways. A generated gene network indicated that Pdgfb was directly connected to and interacted with other genes with similar biologic functions.
Conclusions: A systems genetics analysis revealed that Pdgfb had significant interactions with Atf4 and other genes in MAPK and VEGF pathways, through which Pdgfb was important in maintaining retina function. These findings provided basic information regarding the Pdgfb regulation mechanism and potential therapy for vascular retinopathies.
Platelet-derived growth factor (PDGF) signaling has been implicated in a broad range of vascular retinopathies . Its family consists of four different polypeptide chains encoded by four genes that form five different disulfide-linked dimers: PDGF-AA, PDGF-AB, PDGF-BB, PDGF-CC, and PDGF-DD. They act using two receptor tyrosine kinases: the PDGF receptors (PDGFRs) alpha and beta. Due to their functions in retinopathies, targeted PDGF family members have exhibited potential therapeutic values regarding the diseases . Among these isoforms, PDGFB (Gene ID: 5155, OMIM: 190040) is considered a key regulator of retina diseases and a potential therapeutic target [3,4]. It is a paracrine signal from endothelial cells to pericytes, which are important in regards to capillary morphogenesis in the retina . Indeed, the PDGF-AB complex concentration has been shown to be significantly elevated in patients with diabetes and age-related macular degeneration (AMD) . Further, signaling induced by Pdgfb has been reported to be related to PDGFB/PDGFRβ signaling, which is critical to the formation and maturation of the blood–retinal barrier through the active recruitment of pericytes onto growing retinal vessels . Studies of mice have shown that Pdgfb can be a neuroprotector against light-induced photoreceptor damage  and that the overexpression of Pdgfb in mice has phenotypes similar to diabetic retinopathy . Recent evidence has suggested the ability to simultaneously inhibit VEGF and PDGF to ameliorate AMD; however, clinical trials have not been successful . This could be due to the propensity of PDGF and VEGF family members to form complex regulatory networks . Thus, it is critical to investigate Pdgfb correlated genes and the related co-expression network to understand the function of Pdgfb in the retina.
A systems genetics approach is a unique method for investigating a gene’s coexpression networks. With well-maintained recombinant inbred (RI) strains of mice, a combinatorial application of genomics, transcriptomics, proteomics, and metabolomics can provide detailed insights into gene transcription regulation and subsequent coexpression networks for a gene of interest. Currently, one of the largest panels of RI mouse strains is the BXD family, which consists of more than 150 strains derived from a cross between the C57BL/6J (B6) and DBA/2J (D2) strains. While each line is fully inbred, the entire collection is highly diverse. Many studies of genetic networks have used the BXD RI strains of mice [11-14]. Moreover, the BXD lines have been used extensively in studies that have aimed to identify the genetic bases of ocular phenotypes . Therefore, these strains have a unique foundation with which Pdgfb has been investigated using a systems genetics approach.
This study aimed to assess the expressions and resulting regulatory networks of genes belonging to the PDGF and VEGF families in the retina; it aimed to do so with a systems genetics approach. The intent was to identify the expression quantitative trait locus (eQTL) of Pdgfb, analyze its potential pathways, and construct a genetic network that could contribute to retinal diseases. The data offered several novel putative Pdgfb correlated genes and a novel PDGFB regulatory network in the retina.
The gene expression data generated in this study can be accessed through the “Full HEI Retina Illumina V6.2 (Apr10) RankInv” dataset at the website GeneNetwork . The sections below briefly describe the data generation process.
BXD mice and their parental strains B6 and D2 were bred at the University of Tennessee Health Science Center (Memphis, Tennessee). The animals were fed chow ad libitum and maintained on a 12 h:12 h light-dark cycle. The animal experiment protocols were approved by the Institutional Animal Care and Use Committee of the University of Tennessee Health Science Center. The mice were handled in conformance with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research and the Guide for the Care and Use of Laboratory Animals from the Institute of Laboratory Animal Resources (the Public Health Service policy on the humane care and use of laboratory animals).
Animals 60 to 90 days old were sacrificed using rapid cervical dislocation. After being anesthetized with an IP injection of ketamine, xylazine, and acepromazine at doses of 40, 5, and 1 mg/kg, respectively. From each mouse, two retinas were immediately isolated; they were placed in RNAlater (Applied Biosystems, Foster City, CA), and stored in a single tube overnight at 4 °C, followed by a minimum of 24 h of freezing at −20 °C and long-term storage at −80 °C. Total RNA was extracted from the retinas using the RNA STAT-60 method (Tel-Test Inc., Amsbio, Abingdon, UK) following the manufacturer’s instructions. Briefly, the tissues were homogenized in RNA STAT-60 (1 ml per 50 mg to 100 mg of tissue) with a syringe, followed by an addition of 0.2 ml of chloroform and centrifugation at 12,000 ×g for 1 h at 4 °C. The aqueous phase was transferred into a clean centrifuge tube and incubated with 0.5 ml of isopropanol for 1 h or overnight at −20 °C. After centrifugation at 12,000 ×g for 1 h, the RNA was washed with 75% ethanol and dissolved in 50 μl of nuclease free water (Qiagen, Hilden, Germany). The quality and purity of the RNA was assessed using the Agilent 2100 Bioanalyzer system (Agilent, Santa Clara, CA). RNA samples with an RNA integrity number of eight or higher were selected for further analyses.
The gene expression data used in this study were obtained using Illumina MouseWG-6 v2.0 arrays (Illumina, San Diego, CA) following the manufacturer’s protocol. Briefly, Total RNA was processed with the Illumina TotalPrep RNA Amplification Kit (Applied Biosystems, Foster City, CA) to produce biotinylated complementary-RNAs (cRNAs). This procedure included the reverse transcription of RNA to synthesize the first strand complementary-DNA (cDNA), second strand cDNA synthesis, cDNA purification, in vitro transcription to synthesize cRNA, cRNA amplification, and purification. The biotinylated cRNAs were hybridized to the Illumina Sentrix® Mouse Whole Genome-6 version 2.0 arrays (Illumina) for 19.5 h at 58 °C. The detailed info for RNA extraction and Array hybridization can be found in our previous publication . This dataset included data from 75 BXD strains, their parental strains B6 and D2, and both reciprocal F1 hybrids between B6 and D2, with an average of four samples for each strain (including a roughly equal number of males and females). Raw microarray data were normalized using a rank-invariant method and background subtraction protocols provided by Illumina as a part of the BeadStation software suite (Illumina). In short, the individual p-values was sorted in ascending order, then a rank was assigned for each p-values. The B-H critical value was calculated using the formula formula (i/m) Q, where:
• i = the individual p-value’s rank,
• m = total number of tests,
• Q = the false discovery rate (0.05)
The original p-value was compare to the critical B-H value, and find the largest p value that is smaller than the critical value. The expression data were then renormalized using a modified Z score described in a previous publication . Briefly, log base two was calculated for the normalized values, the Z scores for each array were calculated, the Z scores were multiplied by two, and an offset of eight units was added to each value. This transformation yielded a set of Z-like scores for each array with a mean of eight, a variance of four, and a standard deviation of two.
The IOP of each mouse was measured at the age of six to eight weeks, with an average of eight mice for each strain. The mice were mildly anesthetized with isoflurane (Alfa Aesa, Haverhill, MA), then the IOP was measured from 10:00 a.m. to 4:00 p.m. during the light cycle using an induced impact tonometer (Tonolab, Colonial Medical Supply, Franconia, New Hampshire). The IOP values for the left and right eyes were then averaged. The dataset can be retrieved from a database server with the record ID 12,300 (GeneNetwork).
eQTL mapping was conducted using a WebQTL module in GeneNetwork software (GeneNetwork) using published methods . Briefly, the methodology used regression analyses to determine relationships between differences in a trait and relationships between differences in alleles at markers across a genome. Simple interval mapping was performed to identify potential eQTLs that regulated Pdgfb expression levels and to estimate significance at each location using known genotype data for sites. In addition, composite interval mapping was used to map secondary genetic loci by controlling genetic variances associated with a major eQTL. Both analyses produced likelihood ratio statistics (LRSs); quantitative measurements of links between observed phenotypes, which were variations in expression levels of Pdgfb; and known genetic markers. A genome-wide significance (p<0.05) for each eQTL was established using a permutation test that compared the LRS values of the novel site with the LRS values of 1,000 genetic permutations; other parameters were left as defaults.
The genetic variants of Pdgfb, including single nucleotide polymorphisms (SNPs) and small insertions and deletions (indels) between the parental strains B6 and D2, were determined using the Mouse Genomes Project (Mouse_SnpViewer) variant querying site . Briefly, we searched the genetic variants for Pdgfb at Sanger Mouse SNP Viewer with default settings, except for the “strains” option, where only DBA2/J (D2) was selected (mouse-genomes-project).
The Pearson correlation coefficients of strain means between expressions of Pdgfb and all other probe sets were computed across the mouse genome to produce sets of genetically correlated genes. Next, a literature correlation analysis was performed to rank the list of genes using the software Semantic Gene Organizer, which automatically extracted gene correlations from the titles and abstracts of MEDLINE citations . All genetic and literature correlations were computed using GeneNetwork. Briefly, genes with significant genetic correlations (p<0.05) and literature correlations (r>0.3) with Pdgfb, as well as mean expression levels of >7.0, were selected and uploaded to the open source analysis tools on Webgestalt . This was done for Gene Ontology (GO; biologic processes) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses and transcription factor (TF) target over-representation analyses. The mouse genome was used as the reference gene set. The p values generated from hypergeometric tests were automatically adjusted to account for multiple comparisons using the Benjamini and Hochberg correction method  . This calculation was automatically obtained from Webgestalt  with the formula (i/m) Q, where i is the individual p-value’s rank, m is the total number of tests, and Q is the false discovery rate. Categories with adjusted p values of <0.05 indicated that the set of submitted genes was significantly over-represented in those categories.
PPI networks provided valuable information for understanding cellular functions and biologic processes. The gene set was submitted to the online tool STRING , a database of known and predicted PPIs. The direct PDGFB PPI network was then retrieved with an interaction score of >0.4 (a median confidence).
To identify the associations of related traits with growth factors from the PDGF family, a correlation screening was performed using the currently known members of the PDGF and VEGF families in regards to the IOP. Among the various members of the families, Pdgfb was the only member that was significantly correlated with IOP (r = 0.305; p = 0.012; Figure 1A). The IOP ranged between 11.643 and 19.833 mmHg, with a median of 15.500 mmHg (Figure 1B).
The Pdgfb expression profile was further investigated across the BXD strains. The Pdgfb expressions varied significantly across the BXD strains with a fold change of 1.73 (Figure 2). The average expression of Pdgfb across all the BXD strains was 7.66 ± 0.02 (log2 scale, mean ± SEM), such that BXD16 showed the highest expression (8.09 ± 0.14) and BXD71 showed the lowest expression (7.29 ± 0.12). This indicated that the genetic variability influenced the Pdgfb expressions and thereby facilitated the identification of loci that controlled the Pdgfb expressions among the BXD mice strains.
To identify the genetic loci that regulated Pdgfb expressions in the retina tissues, simple interval mapping was performed for Pdgfb across the mouse genome. The results showed a significant eQTL on chromosome 15 at 80.25 Mb with an LRS of 39.3, which was close to the physical location of Pdgfb on chromosome 15 at 80.00 Mb (Figure 3). Further composite interval mapping revealed no secondary genetic loci that modulated Pdgfb expressions. This indicated that Pdgfb was a cis-regulated gene, which suggested that a variation in Pdgfb expression was caused by genetic variants nearby or within the Pdgfb gene. Further, a search for mouse genome variants from the Mouse Genomes Project yielded a total of 165 genetic variants in Pdgfb between the BXD parental strains, including SNPs and indels. Among them, 18 variants were defined as 5′-UTR UTR-3′ synonymous or splice variants (Table 1), while the rest were all intron variants. The 5'- untranslated region (UTR) and 3'-UTR we mentioned in the text are genomic regions rather than primers. 5' UTR is the portion of an mRNA from the 5' end to the position of the first codon used in translation. The 3' UTR is the portion of an mRNA from the 3' end of the mRNA to the position of the last codon used in translation. Some of these variants could have been responsible for the variable Pdgfb expression in the BXD mice.
After filtering was conducted based on a mean expression value, genetic correlation, and literature correlation, a total of 2,477 genes (p<0.05) were identified as correlated with Pdgfb (Appendix 1). The top ten positively and negatively correlated genes are listed in Table 2. Among them, Atf4 was expressed significantly in the retinas, and it was the top gene positively correlated with Pdgfb (r = 0.628; p = 1.29e-10; Figure 4). To explore the functional relationship between Atf4 and Pdgfb, a TF target enrichment analysis was performed and three motifs were identified: CTGCAGY, RAAGNYNNCTTY, and TGAMCTTTGMMCYT. They were colocated in the transcription start sites of Atf4 and Pdgfb (FDR p<0.05). However, according to the TRANSFAC database, the first two motifs did not match any known TF binding sites, but the last motif did match the annotation for HNF4A. That is, the gene encoded hepatic nuclear factor 4 alpha, which was essential for maintaining normal vision and eye phenotypes. According to Mouse Genome Informatics and the International Mouse Phenotyping Consortium (IMPC), mice that carry a null allele of this gene have abnormal retinal vasculatures, retinal blood vessels, and lens morphologies and increased total retina thicknesses and cataracts.
To further evaluate the gene function of the Pdgfb correlated genes, 2,477 genes were submitted to Webgestalt to perform GO and KEGG pathway enrichment analyses. The gene set was found to be significantly associated with the development of blood vessels, eyes, and retinas (Table 3). Furthermore, according to the KEGG pathway enrichment analyses, several pathways involved in AMD were overrepresented, including the MAPK signaling pathway, the PI3K-Akt signaling pathway and the apoptosis pathway (Table 3).
To further identify the proteins that directly interacted with PDGFB and its network, the gene set of 2,477 genes was submitted to the STRING database. The results showed that 45 proteins directly interacted with PDGFB, among which 17 proteins were already implicated in eye and vision related traits, according to mammalian phenotype ontology. These included CTGF, EGF, EGFR, FGFR1, FGFR2, HIF1A, JAG1, KDR, KIT, NRP1, PDGFC, PR1R, PTEN, S1C20A2, TGFB1, TRP53, and VEGFB (Figure 5).
Correlations between members of the PDGF and VEGF families and IOP, a major risk factor for vascular retinopathy, were evaluated. Among the members of these families, Pdgfb was the only gene in the retina with an expression that correlated with IOP. Previous reports indicated elevated levels of Pdgfb in the plasma of neovascular AMD individuals . Moreover, several genetic markers for IOP were suggested in population GWAS studies . However, the genetic markers of Pdgfb had not yet been elucidated. In the present study of BXD mouse strains, a significant correlation was identified between the expression of Pdgfb and IOP. This result agreed with the neuroprotection function of Pdgfb in physiologic conditions ; it had been found that in animal models, a high retinal expression of Pdgfb resulted in a tractional retinal detachment due to a proliferation of vascular and nonvascular cells, similar to diabetic retinopathy in humans . Thus, the present study’s results indicated the importance of Pdgfb in regards to retinal vascular disease and thereby reinforced the possibility of targeting Pdgfb for the treatment of AMD and glaucoma.
In addition to the known Pdgfb correlated regulators, the results of the gene correlation analysis suggested a group of putative Pdgfb correlated genes. Among them, Atf4 was expressed significantly in the retinas and was the top gene positively correlated with Pdgfb. In addition, expressions of Atf4 and Pdgfb had significant correlations in many human tissues (Appendix 2). Atf4 was characterized as cAMP-response element binding protein 2 (CRBP2), a DNA binding protein widely expressed in mammals . The function of Atf4 in normal and malignant cells was related to stress-induced transcription activation . In smooth muscle cells, Pdgfb could mediate cell migration through Atf4 . However, the relationship between Atf4 and Pdgfb in the retina remained unclear.
One potential mechanism was that Atf4 and Pdgfb shared common transcriptional factor binding motifs and were thus coexpressed. In the TF binding enrichment analysis, Hnf4a was identified as a shared TF for both Atf4 and Pdgfb. Mice with null Hnf4a presented retina disorders, but the mechanism remained unclear. The results indicated a potential mechanism through the PDGFB pathway.
Another shared regulator could be noncoding RNAs. MicroRNA-214 (miR-214) was a microRNA that regulated ocular neovascularization . Both Atf4 and Pdgfb were known targets of miR-214. Moreover, miR-214 inhibited vascularization through a reduced Pdgfb release, and a suppression of miR-214 could increase an Atf4 level . However, the detailed mechanisms needed to be further illustrated. Augmented Atf4 signals that occurred during retinal degeneration had cytotoxic roles by triggering photoreceptor cell death . Atf4 was overexpressed in retinas that suffered from oxygen-induced retinopathy, and its inhibition had been reported to prevent pathologic neovascularization . Furthermore, Atf4 had been reported to bind to a VEGF promoter , which suggested a complicated regulatory network involving ATF4-PDGFB-VEGF. The results indicated a critical function of Atf4 in regards to Pdgfb and VEGF in the retina.
One of the top genes negatively correlated with Pdgfb was Hdgf. Although the cross-regulation of Pdgfb and Hdgf remained unclear, Hdgf had been found to display a functional similarity with Pdgfb. Hdgf had stimulated smooth muscle cell growth and had acted as a mitogen for various types of cells. It had also recently been identified as an angiogenic factor that regulated the retinal vasculature in physiological and pathological conditions . In injured retinal ganglion cells, HDGF could have a neuroprotective role . Thus, the negative correlation between Hdgf and Pdgfb indicated potential compensatory roles for the two growth factors. However, more studies were needed to elucidate the detailed mechanism.
A gene enrichment analysis was used to identify the GO categories in which Pdgfb and its correlated genes were overrepresented and to identify the pathways through which they had roles in regards to the retina. The GO analysis showed that Pdgfb correlated genes were significantly involved in blood vessel development, retina development, differentiation, and morphogenesis. This confirmed the major functions of Pdgfb and its correlated genes in the retina. Further, the gene pathway enrichment analysis indicated MAPK as the topmost signaling pathway, and this result agreed with a previous study that showed Pdgfb had the ability to activate the MAPK pathway . Moreover, in the present study, it was found that 98 other Pdgfb correlated genes were located in the MAPK pathway, and Pdgfb interacted with them to maintain various functions, such as angiogenesis and cell survival.
Finally, the VEGF signaling pathway was one of the top enriched pathways. The cross interaction between VEGF and PDGF signaling had been noted in previous studies , and PDGF signaling had been suggested to be partially responsible for resistance to VEGF targeting therapy. Thus, the present study’s results confirmed this cross-family interaction in the retina and therefore provided a helpful resource for further investigations into VEGF-PDGF signaling interactions.
To further illustrate the function of the Pdgfb correlated genes at the protein level, a PPI network analysis was performed. The PPI analysis provided protein interaction information for the gene products. As a result, 57 gene products were reported to have directly interacted with PDGFB, among which 17 proteins were involved in eye diseases. Intriguingly, most of these proteins were growth factors and corresponding receptors. These results indicated a complex growth factor cotranscription regulatory network. In particular, two members from the PDGF family, PDGFA and PDGFC, had significant correlations with PDGFB, which confirmed that the PDGF family was involved in the regulation of PDGFB expressions. Further, a correlation with the VEGF family member VEGFB was found. Growth factors from the VEGF and PDGF families were promising therapeutic targets for major vascular retinopathies. The co-expression network could provide novel leads for therapy and drug design, and this work confirmed that the complicated gene network of the VEGF-PDGF interaction had a critical role in the retina.
With a systematic analysis that involved more than 70 mouse strains, it was found that a Pdgfb expression in the retina was highly variable and correlated significantly with IOP. A Pdgfb coexpression network revealed a complicated network of growth factors in the retina, and Atf4 and Hdgf were potential novel genes that interacted with Pdgfb in the retina. Further gene ontology and pathway analyses indicated that many of these genes interacted with the MAPK and VEGF pathways and thus provided the putative Pdgfb action mechanism of mediating cell survival and providing neuroprotection. Thus, the study confirmed the importance of Pdgfb in retinopathies. However, it was noted that further investigations should be conducted to elucidate the nature and underlying mechanisms of these relationships and assist in the development of efficacious treatment strategies for retinal vascular diseases.
Appendix 1: The list of 2477 genes (p<0.05) correlated with Pdgfb.
Appendix 2. Scatterplot of correlation between the expression of ATF4 and PDGFB in the human tissues from Human GTEx data set.
Appendix 3: The full list of GTEx tissue name.
The research was supported by Taishan Scholars Construction Engineering (G.T.), National Natural Science Foundation of China (31771284, 81670855), Open Research Funds of the State Key Laboratory of Ophthalmology (G.T. and X.L.), Science and Technology Support Plan for Youth Innovation of Higher Education of Shandong (2019KJE013), the Key Research and Development Plan of Shandong Province (2016GSF201100, 2018GSF118230), Shandong Provincial Natural Science Foundation(ZR2016JL026), Swedish Research Council grant 2015-4870 (JB), A grant from National Institutes of Health of United States (R01DC015111). The authors would like to express their gratitude to EditSprings (https://www.editsprings.com/) for the expert linguistic services provided. Dr. Geng Tian (email@example.com) and Dr. Xuri Li (firstname.lastname@example.org) are co-corresponding authors for this paper.