Identification of potential crucial genes and pathways associated with vein graft restenosis based on gene expression analysis in experimental rabbits

Occlusive artery disease (CAD) is the leading cause of death worldwide. Bypass graft surgery remains the most prevalently performed treatment for occlusive arterial disease, and veins are the most frequently used conduits for surgical revascularization. However, the clinical efficacy of bypass graft surgery is highly affected by the long-term potency rates of vein grafts, and no optimal treatments are available for the prevention of vein graft restenosis (VGR) at present. Hence, there is an urgent need to improve our understanding of the molecular mechanisms involved in mediating VGR. The past decade has seen the rapid development of genomic technologies, such as genome sequencing and microarray technologies, which will provide novel insights into potential molecular mechanisms involved in the VGR program. Ironically, high throughput data associated with VGR are extremely scarce. The main goal of the current study was to explore potential crucial genes and pathways associated with VGR and to provide valid biological information for further investigation of VGR. A comprehensive bioinformatics analysis was performed using high throughput gene expression data. Differentially expressed genes (DEGs) were identified using the R and Bioconductor packages. After functional enrichment analysis of the DEGs, protein–protein interaction (PPI) network and sub-PPI network analyses were performed. Finally, nine potential hub genes and fourteen pathways were identified. These hub genes may interact with each other and regulate the VGR program by modulating the cell cycle pathway. Future studies focusing on revealing the specific cellular and molecular mechanisms of these key genes and pathways involved in regulating the VGR program may provide novel therapeutic targets for VGR inhibition.


INTRODUCTION
Occlusive artery disease is a major cause of morbidity and mortality worldwide (Bansilal, Castellano & Fuster, 2015). Despite development of novel treatments in past decades, coronary artery bypass graft (CABG) surgery remains the standard of care for patients with left main coronary artery disease (CAD) and three-vessel CAD (Serruys et al., 2009). In contrast, most patients with late-stage peripheral artery occlusive disease are treated with peripheral artery bypass graft surgery (Weintraub et al., 2012). Due to their advantages in availability and length, veins are the most commonly used conduits in coronary and peripheral artery vascular surgeries (Goldman et al., 2004). However, one major issue that has dominated the field for many years is that autologous veins are especially prone to failure. Although data have shown that the patency rates of vein grafts diminish from 98% immediately after surgery to <88% within the first month post-surgery and to 60% at 10 years after surgery. To date, no unequivocally effective treatments are available for vein graft failure (Elmore et al., 2016). Thus, clarifying the cellular and molecular mechanisms involved in VGR and identifying potential novel therapeutic targets for the prevention of restenosis are significant goals.
At present, vein grafts are thought to undergo an adaptation process in new arterial environments during which intimal thickening and structural vessel wall remodeling appear (Buccheri et al., 2016;De Vries et al., 2016). However, the cellular and molecular mechanisms underlying intimal hyperplasia and the vascular remodeling process are largely unknown. Several studies have identified individual genes involved in restenosis; for instance, early growth response protein 1 (Egr-1) was suggested to promote endothelial cell(EC) proliferation and induce vein graft restenosis by up-regulating ICAM-1 expression (Zhang et al., 2013). Tumor necrosis factor alpha-stimulated gene 6 (TSG-6) was indicated to suppress restenosis of vein grafts in rats by inhibiting the inflammatory response . Despite the rapid development of genome sequencing and microarray technologies, very few studies have characterized genes and pathways associated with VGR at the transcriptome level, which may be partly due to the challenge of obtaining appropriate samples and establishing animal models. Kalish et al.(2004) identified the transcriptome responses of canine vein bypass grafts using transcriptional profiling and defined the individual contributions of ECs and VSMCs in a cell-specific vein graft transcriptome analysis (Bhasin et al., 2012). However, a main weakness of the former study was the unavailability of canine-specific gene arrays in 2004, as mentioned by the authors (Bhasin et al., 2012). The most important limitation of the latter study lies in the fact that only the role of ECs and VSMCs were taken into consideration. However, other cells and molecular components, such as inflammatory cells and immune cells, also play pivotal roles in mediating VGR (De Vries et al., 2016).
The purpose of the current study was to determine potential crucial genes and pathways involved in VGR. A bioinformatics analysis was performed with gene expression data from our previous study (Wang et al., 2016), which avoided the limitations of the aforementioned microarray data. After screening differentially expressed genes (DEGs) and performing a functional enrichment analysis, protein-protein interaction (PPI) networks of the DEGs were constructed and visualized. The present study makes notable contributions to our understanding of the molecular mechanisms involved in VGR and provides valuable biological information for further exploration of potential candidate biomarkers and therapeutic targets for the prevention of VGR.

Microarray data
The gene expression data used in this study were based on the Agilent GPL7083 platform (Agilent Rabbit 4x44K Gene Expression Microarrays). The whole rabbit genome oligo microarray provides a broad view by representing all known genes and transcripts in the rabbit genome. Sequences were compiled from a broad source survey and then verified and optimized by alignment to the assembled rabbit genome. As previously described, significant vascular wall thickening and a high level of cell proliferation were observed, whereas cell apoptosis was at the lowest level, in seven days after surgery group (Wang et al., 2016). Hence, seven days after surgery was considered a key time point for vein graft restenosis. Based on these observations, we compared the day 7 group with the control group as the best choice to identify DEGs involved in VGR. Raw microarray data (GSE110398) produced in our previous study were remodeled and further analyzed. A total of 10 samples were extracted from our previous dataset (Wang et al., 2016), including six vein graft samples (C1, C2, C3, C4, C5, and C6) removed seven days after surgery and four vein graft samples (O1, O2, O3, and O4) obtained from the sham surgery group (data can also be found in the Supplemental File). Herein, samples obtained seven days after surgery were defined as the surgery group, and samples obtained from the sham surgery group were defined as the control group.

Identification of DEGs
The raw signal intensities were normalized with the quantile method by GeneSpring GX v11.5, and low intensity genes were filtered. The probe quality control was determined using principal component analysis (PCA) in GeneSpring. The statistical software R (version 3.4.1; R Core Team, 2017) and the linear models for microarray data package (limma) (Ritchie et al., 2015) in Bioconductor (http://www.bioconductor.org/) were applied to identify DEGs by comparing expression values between samples in the surgery and control groups. An empirical Bayes method was used to select significant DEGs based on the ''limma'' package in Bioconductor. In this study, we were interested in determining which genes were expressed at different levels between the control and surgery groups. In our analysis, linear models are fitted to the data with the assumption that the underlying data are normally distributed. Therefore, a design matrix was set up with the vein graft tissue information.

Enrichment analyses of DEGs
The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 is an integrated functional annotation tool for investigators to examine the biological meanings underlying a large list of genes (Huang, Sherman & Lempicki, 2009). GO (gene ontology) (Gene Ontology Consortium, 2006) and KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa & Goto, 2000) pathway enrichment analyses of DEGs were performed using the authoritative online tool DAVID (https://david.ncifcrf.gov/). An adjusted p-value <0.05 and count ≥2 were considered to have achieved significant enrichment.

PPI network construction
The Search Tool for the Retrieval of Interacting Genes (STRING) database is an online tool designed to provide a critical assessment and integration of protein-protein interactions, including both direct (physical) and indirect (functional) associations (Szklarczyk et al., 2014). To evaluate the potential relationships among DEGs, we mapped the DEGs into STRING; only experimentally validated interactions with a combined score >0.4 were considered significant. Moreover, the Molecular Complex Detection (MCODE) app was used to identify modules of the PPI network in Cytoscape (Shannon et al., 2003) with a degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max. depth = 100. MCODE scores >3 and a number of edges >40 were set as the cut-off criteria. Pathway analysis of genes in each module was performed using the online tool DAVID, and p-values <0.05 were considered significant.

Identification of DEGs
The statistical analysis software R (version 3.4.1; R Core Team, 2017) was used to identify differentially expressed genes using an adjusted p-value <0.05 and |log2FC (fold change)|≥1 ( Fig. 1) as criteria. A PCA plot was generated using normalized data. The PCA plot shows clear separation between the surgery and control group samples (Fig. 2). A total of 858

GO term enrichment analysis
The GO functional enrichment analysis showed that the up-regulated genes were significantly enriched in biological processes (BPs), including immune response, inflammatory response and innate immune response (Table 1). The down-regulated DEGs were significantly enriched in positive regulation of the ERK1 and ERK2 cascades, positive regulation of the apoptotic process, and chondrocyte differentiation (Table 2). For the cellular component (CC) category, the up-regulated DEGs were significantly enriched in extracellular exosome, extracellular space and membrane (Table 1). Interestingly, the down-regulated DEGs were also enriched in extracellular exosome, cytoplasm, and extracellular space (Table 2). For the molecular function (MF) category, the up-regulated DEGs were enriched in ATP binding, transmembrane signaling receptor activity, and kinase activity (Table 1), whereas the down-regulated DEGs were mainly enriched in calcium ion binding, heparin binding, and protein tyrosine phosphatase activity (Table 2).

KEGG pathway analysis
In the KEGG pathway enrichment analysis, a total of 60 up-regulated pathways and 26 down-regulated pathways were obtained using the functional annotation tool DAVID. However, some enriched pathways evidently had no associations with VGR. For example, the most significant down-regulated pathway is pathways in cancer. Hence, a total of 14 potential key pathways were selected from these enriched pathways by literature review and general experience ( Table 3); most of these pathways are involved in regulating cell proliferation, apoptosis, vascular smooth muscle contraction, and cell adhesion, which are processes that may be highly associated with VGR. Moreover, the 164 DEGs contained in these 14 key pathways were extracted and defined as the ''KEGG genes'' gene list.

PPI network analysis
Based on the information in the STRING database, the PPI network consists of 695 nodes and 2,098 interactions. The top 34 hub nodes with more than 20 degrees were screened and defined as the ''PPI hub genes'' gene list. Moreover, nine hub genes obtained from the intersection of the ''KEGG genes'' and ''PPI hub genes'' gene lists (Fig. 4) are listed in Table 4. In addition, a comprehensive PPI network was constructed by merging the 9 PPI networks of the nine hub genes and visualized using the Cytoscape software (Fig. 5). A total of 695 nodes and 2,098 edges were analyzed using the plug-in MCODE, and the top three significant modules were screened (Fig. 6).
To further lend some more validity to the analysis results, the day 7 versus day 1 comparison was served as a control relative to day 7 versus sham surgery comparison. The expression of these nine hub genes in three groups (day 1, day 7 and control) is shown in Fig. 7. As expected, most hub DEGs identified in day 7 versus sham surgery group, are also differentially expressed in day 7 versus day 1 group, which may due to the influence of surgery hit, and all hub DEGs are also differentially expressed in day 1 versus sham surgery group.

DISCUSSION
Improving understanding of the cellular and molecular mechanisms underlying VGR is critically significant for reducing the rates of vein graft failure. In the present investigation, an integrated bioinformatics analysis was performed to identify potential key genes and pathways involved in regulating VGR. In summary, a total of 858 up-regulated and 817 down-regulated DEGs were screened. A total of 14 potential crucial pathways were selected from the functional annotation enrichment analysis. In the PPI network analysis, a comprehensive PPI network of the 9 hub genes was constructed. Furthermore, three significant modules were screened from the PPI network by module analysis. All nine hub genes showed higher degrees in the PPI network analysis and were involved in potential key pathways. Several previous studies demonstrated that fibroblast growth factor receptor 1 (FGFR1) activated the Akt/mTOR pathway (Chen & Friesel, 2009) and regulated phenotypic modulation of vascular smooth muscle cells (VSMCs) (Chen & Friesel, 2009;Chen, Simons & Friesel, 2009). As described in our previous study, the mTOR signaling pathway was indicated to play a major role in the arterialization of vein grafts (Wang et al., 2016). Moreover, a fibroblast growth factor receptor antagonist was validated to block growth factor-mediated VSMC proliferation by inhibiting binding of fibroblast growth factor 2 (FGF2) to VSMCs and soluble FGFR1 (Segev et al., 2002). Table 3 Potential key pathways selected out from KEGG pathway enrichment analysis of DEGs.

Category Term Count P Value Genes
Up-regulated Cytokine-cytokine receptor interaction  32  4.04E−06  CSF2, IL1R2, TNFRSF21, TNF, CCL2, CCR1, CXCL8,  IL15, IL7R, CCL4, IL17RA, CXCL10, TNFRSF11B, CXCR4,  IFNG, IL1B, FAS, XCR1, LTB, IFNGR2, IL1A, IL4, IL18R1,  IL6, TNFSF4, LTBR, TGFBR1, LOC100348776, CCR8,  TNFSF10, TNFSF13B,  VSMC proliferation was also thought to play a critical role in the pathophysiology of restenosis. The second hub gene PIK3CG (phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit gamma) was shown to be involved in three key pathways: the TNF signaling pathway, Toll-like receptor signaling pathway, and chemokine signaling pathway (Table 4). CDK1 (cyclin-dependent kinase 1) is a protein-coding gene that plays a key role  in regulating the eukaryotic cell cycle by modulating the centrosome cycle and mitotic onset and thus cell proliferation. Cell proliferation requires the expression of CDKs, and CDKs are considered a proliferation signaling cascade of VSMCs (Schad et al., 2011). Evidence indicated that radiation inhibited VSMC proliferation via cell cycle arrest by enhancing p21 expression and suppressing CDK1 and 2 (Kim et al., 2004). Consistently, a derivative of the CDK inhibitor roscovitine was validated to potently induce G1 phase arrest and therefore inhibit VSMC proliferation (Sroka et al., 2010). YWHAE (tyrosine 3monooxygenase/tryptophan 5-monooxygenase activation protein epsilon) has 34 degrees in the PPI network and is also included in the cell cycle pathway. Moreover, YWHAE belongs to the 14-3-3 family of proteins that induce signal transduction by binding to phosphoserine-containing proteins. Intriguingly, a recent study provided evidence that up-regulation of YWHAB in endothelial cells played a pivotal role in intimal hyperplasia Figure 5 A comprehensive PPI network were constructed by merging nine PPI networks of the nine hub genes respectively. Each node corresponds to a DEG, and edges represent the interactions between DEGs. DEGs with higher degrees appear larger size. The gradual color from blue to red represents the changing process from down-regulation to up-regulation. Full-size DOI: 10.7717/peerj.4704/ fig-5 following carotid artery injury via enhancing endothelial cell proliferation and migration (Feng et al., 2017). YWHAG knockdown in zebrafish was reported to reduce the brain size and increase the diameter of the heart tube (Komoike et al., 2010). CCNB2 (cyclin B2) is a member of the cyclin family, specifically the B-type cyclins. Evidence indicated that IL-18 promoted cell proliferation via NF-κB and the p38/ATF2 pathway by targeting CCNB2 (Zhang et al., 2015). Down-regulation of MAD2L1 was demonstrated to be involved in suppressing the proliferation, migration, invasion, apoptosis induction and cell cycle arrest of cancer cells (Li, Bai & Zhang, 2017). PCNA (proliferating cell nuclear antigen) is a molecular marker for proliferation due to its role in replication and is involved in both the cell cycle pathway and DNA replication; PCNA had 30 degrees in the PPI network. Previous studies showed that PCNA expression was correlated with stent-induced in-stent restenosis (Gao et al., 2013) and proposed that PCNA might be a potential future target (Wang, 2014). Guanine nucleotide-binding protein subunit beta-4 (GNB4) was involved in the chemokine signaling pathway and showed 20 degrees in the PPI network. Surprisingly, we noted that six of the nine hub genes (CDK1, YWHAE, YWHAG, CCNB2, MAD2L1, and PCNA) were involved in the cell cycle pathway, which further supported the enormous significance of this pathway in regulating VGR (Table 4). The key role of the cell cycle pathway in regulating VGR has been established (Sriram & Patterson, 2001); Hierarchical clustering analysis of DEGs in cell cycle pathways was shown in Fig. 8. However, the underlying molecular mechanisms remain largely unknown. Based on the abovementioned evidence, these hub genes may potentially interact with each other and regulate the VGR program by modulating the cell cycle pathway. Future studies focusing on revealing the specific cellular and molecular mechanisms involved in regulating the VGR program may provide novel therapeutic targets for VGR inhibition.
Over the past few decades, many powerful research tools in genomic technologies have advanced rapidly, such as genome sequencing and microarray technologies. Microarray technology has enhanced our ability to screen samples for thousands of genes at once, which may provide novel insights into potential molecular mechanisms involved in regulating the VGR program. Unfortunately, high throughput data associated with VGR are extremely scarce at present. ''Restenosis'' was used as a search term to search in the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/), and only 7 datasets were returned; this lack may be partly due to the challenge of obtaining appropriate samples in vivo and establishing a VGR animal model. Hence, we suggest that collaborations between clinicians and basic researchers may accelerate our understanding of the cellular and molecular mechanisms involved in VGR and therefore improve the clinical efficacy of occlusive artery disease.
Due to the extreme lack of high throughput data associated with VGR, more research is required to screen key genes and pathways involved in this process. Further experimental investigations are needed to estimate the role of these candidate genes and pathways. However, there are certain limitations in the present study, the aforementioned results, including the gene expression levels and their functions, were needed to be validated by experiments, which would be carried out in our further studies.

CONCLUSIONS
In conclusion, this study provides a set of potential therapeutic targets for future investigations into the molecular mechanisms involved in VGR. However, we acknowledge that there were certain limitations in this study. The aforementioned results, including the expression level of genes and the definite functions of potential hub genes and key pathways, were needed to be validated by experiments, which would be conducted in our further studies.