A systems genetics approach to revealing the Pdgfb molecular network of the retina

Generating retinal gene expression data from BXD mice

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 [11]. The sections below briefly describe the data generation process.

Animals

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).

Tissue harvesting and RNA extractions

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.

Microarrays and data processing

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 [16]. 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 [17]. 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.

Intraocular pressure (IOP) measurements

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

eQTL mapping was conducted using a WebQTL module in GeneNetwork software (GeneNetwork) using published methods [17]. 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.

Sequence variant analyses

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 [18]. 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).

Gene set enrichment analyses

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 [19]. 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 [20]. 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 [21] . This calculation was automatically obtained from Webgestalt [20] 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.

Protein–protein interaction (PPI) network constructions

PPI networks provided valuable information for understanding cellular functions and biologic processes. The gene set was submitted to the online tool STRING [22], 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).

Significant correlation of Pdgfb with IOP

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).

Pdgfb expression levels in the retinas of BXD mice

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.

eQTL mapping and the sequence variants of Pdgfb

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.

Pdgfb correlated genes and gene function enrichment analyses

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).

PPI network analyses

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).

Conclusions

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.