Co-expression network analysis identified hub genes critical to triglyceride and free fatty acid metabolism as key regulators of age-related vascular dysfunction in mice

Network construction and module detection

After microarray data preprocessing, a total of 2,256 differentially expressed Agilent transcripts were left. We further removed ambiguous probes and duplicated genes. The ambiguous probes included positive and negative controls in the microarray, and the duplicated genes were removed according to their gene symbols. Finally, the expression values of the 469 differential genes were selected for subsequent analysis. In graph theory, the nodes with higher connections are more likely to be pivotal connectors. Likewise, critical points controlling important dynamic components in biological networks [21], removal of them may cause biological systems failure in saving their coherence. In order to identify such interconnected genes from aging-associated co-expression network abstractly, adjacency matrix was obtained from the Pearson correlation matrix with a power β=12 based on the scale-free topology criterion [21]. Subsequently, the represented matrix was transformed to similarity matrix in clustering process during which genes grouped in a cluster were likely to be biologically relevant to the same pathway. The value of power β=12 could emphasize robust correlations and remove unreliable correlations between genes on an exponential level. Figure 1A shows the determination of β parameter based on the description in the WGCNA manual. Briefly, given Agilent series GSE50833, we used the WGCNA to establish a number of modules in co-expression network, and 6 modules were obtained (Figure 1B). As illustrated in Figure 1B, modules in gene dendrogram are shown in different colors and based on dynamic branch cutting algorithm underneath row color assigns the modules membership.

Figure 1. Parameter analysis of inferred co-expression network and modules. (A) A scaling factor beta determination based on the scale-free topology criterion. (B) Hierarchical clustering of genes in significant modules. The colors are assigned to each module by the Dynamic Tree Cut algorithm.

Modules associated with vascular aging and lipid levels

We aimed to gain new insights into the roles of specific modules that might be involved in vascular aging, module eigengenes associated with phenotypic traits were assessed in this study. An eigengene is referred to as gene expression profiles inside each of the modules summarized by the first principle component. The turquoise and blue modules that are positively correlated with triglyceride, leptin, and free fatty acid (FFA) respectively, were established to be significantly associated to vascular aging stage with an absolute correlation coefficient of > 0.6 and p-value of < 0.05. Figure 2A shows a highly significant correlation between gene significant (GS) versus module membership (MM) in the turquoise module with Triglyceride, and Figure 2B illustrates the correlation between turquoise module size and triglyceride. Next, we examined whether significant modules (turquoise and blue) were enriched for GO terms relevant to age-related vascular disorders. Gene ontology enrichment analysis for the turquoise module is presented in Table 1 and the GO analysis of other modules are presented in Supplementary Table 5. From Table 1, we noticed that 44 genes are involved in the proteolysis biological process, and the products of 46 genes are located under membrane-bounded vesicle.

Figure 2. Module-traits and module_membership-gene_significance correlation analyses. (A) Scatterplot shows a highly significant correlation between gene significant (GS) versus module membership (MM) in the turquoise module with Triglyceride. (B) Heatmap shows correlation between assayed traits and module eigengene values. Green and red colors represent the negative and positive correlation respectively. Decimals outside of round brackets are correlation, and decimals inside of round brackets stand for gene significance level.

Inference of miRNA–mRNA-transcription factor interaction

The top 30 interconnected genes within each of the turquoise, blue and red modules were graphically depicted by utilizing Cytoscape (Supplementary Table 17, Supplementary Table 19, Supplementary Table 21). Interestingly, miRNA-target analysis showed that in the turquoise module, gene F3 is putatively targeted by 23 and 38 predicted miRNAs from TargetScan and MicroCosm respectively (Figure 3A). Additionally, in the blue module, gene H2-Q7 illustrated to be targeted by 78 miRNAs in MicroCosm (Figure 3B). No genes were found to be targeted with miRNAs in the red module. Among the miRNAs, mmu-miR-449a, mmu-miR-449c, mmu-miR-34c, mmu-miR-34b-5p, mmu-miR-15a, and mmu-let-7 were listed as common between networks built by turquoise and blue modules. Due to the inconvenience of considering a highly dimensional set of transcription factors, we only investigated transcription factors in a network of the top 90 connected hub nodes depicted for mRNA-miRNA interaction analysis (Supplementary Table 9), when transcription factors Pax8 and Hsf1 were ranked as the most significant based on the normalized enrichment score test (NES) (Table 2).

Figure 3. The network of top 30 interconnected genes in the turquoise (A) and blue (B) modules and predicted miRNAs from TargetScan and Microcosm. The hub genes are shown in red and miRNAs in light yellow. mRNA-miRNA interactions from TargetScan and Microcosm are shown in green and blue lines respectively. Yellow circles show the genes that are putatively regulated by miRNAs.

Module preservation test

To test if the age-dependent modules are preserved between the constructed network and a reference network, we built a reference network based on 1036 significantly differentially expressed transcripts with P < 0.05 from 6-month old wild-type mice (n=3) and 20-month old wild-type mice (n=3) in the GSE10000 dataset (Supplementary Table 3). Among the 2,256 transcripts from GSE50833 datasets, 125 genes (according to the gene symbols) were matched with modules defined by Ghazalpour et al. [22] (Supplementary Table 4), by which we had already built a reference network. Median rank and Z-summary values were calculated to conduct module preservation test, and both measures revealed high preservation of green and blue modules between the two networks (Supplementary Figure 5A and 5B). Largely overlapped interactions between the two networks were found in this study, indicating the robustness of the built modules is strong. The strong module preservation could be due to the overlapped interactions between the two networks, which indicate that the gene interactions within a module are conserved.

qPCR validation of the key genes

To verify the main conclusion drawn from the microarray results, the relative expression levels of the five key genes (Enpp5, Fez1, Kif1a, F3, and H2-Q7) and their interacting microRNAs (mmu-miR-449a, mmu-miR-449c, mmu-miR-34c, mmu-miR-34b-5p, mmu-miR-15a, and mmu-let-7) were determined using qPCR. No significant difference in body weight was detected between the two groups of male mice. The RT-qPCR analysis results indicated that Enpp5, Fez1, Kif1a, F3, mmu-miR-34c, mmu-miR-34b-5p, and mmu-let-7 were significantly up-regulated in samples from 20-month old mice group compared with the 6-month old mice group, whereas mmu-miR-449a and mmu-miR-449c were significantly down-regulated. Meanwhile, there were no significant differences observed in the H2-Q7 and mmu-miR-15a genes between the two groups, but the expression level of H2-Q7 is a little higher in the group of old mice compared to the young mice.

Data collection

We carefully selected vascular aging-associated dataset from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, and identified two time-series microarray datasets regarding vascular aging in mice aorta (GEO access number: GSE50833; GSE10000). The Agilent TXT files for GSE50833 and Affymetrix TXT files for GSE10000 were downloaded from the GEO database respectively. The data GSE50833 [16] was used as first-stage discovery set for co-expression network construction, and GSE10000 [17] as a second-stage independent verification dataset for module preservation test. This dataset GSE50833 consists of a total of 12 samples with data generation platform of GPL10787, and these samples correspond to 6-month old mice (n=6) and 20-month old mice (n=6). The data GSE10000 consists of a total of 18 samples with data generation platform of GPL1261 and GPL8321 respectively, and these samples correspond to three pairs of 6-week old wild-type and ApoE-/- mice, three pairs of 32-week old wild-type and ApoE-/- mice, and three pairs of 78-week old wild-type and ApoE-/- mice.

Data preprocessing and identification of DEGs

Statistical software R (version 3.5.1) and its related packages were used for data preprocessing and identifying nominal DEGs between the two aging stages. The affyPLM package could fit the original data of the microarray to generate the weights and residuals diagram, the relative log expression (RLE), and the relative standard deviation (NUSE, Normalized unscaled standard errors) box diagram. Before we analyzed the data, we conducted quality testing of microarray data. In this process, we used some powerful and accurate R packages, such as affyPLM, affy, and RColorBrewer. After conducting the procedure of removing unqualified samples or probes, we got a reasonable and useful sample set. Raw files were normalized with quantile method (Supplementary Table 1) and the limma package (Linear Model for Microarray Data) was applied to extract genes with P < 0.05 by comparing the gene expression values between the 6- and 20-month old mice. These genes were considered as nominal DEGs without multiple test correction. Ultimately, the extracted Agilent probe IDs of the identified DEGs for the data GSE50833 were transformed into Agilent MIT IDs (Supplementary Table 2) and used for follow-up analysis.

Network construction

After removing the redundant probes that were not able to contribute greatly to the construction of modules, the nominal DEGs were used to build weighted gene co-expression network through the WGCNA R package, by which the calculated absolute value of the Pearson correlation coefficient for all pair-wise comparisons of gene expression values was transformed into a similarity matrix. Next, the similarity matrix was applied to stepBystep, a network construction, and module detection function. We carried out this step by building a weighted matrix with a scaling factor-beta based on the assumption that biological networks are scale-free [43]. The modules were computed by assigning a minimum of 30 genes per module and keeping the default value of SplitDepth in two for a medium sensitivity of cluster splitting. These parameters would optimize scale-free topology and robust node connectivity criteria when any two genes were connected, and the edge weight was determined by the topological overlap matrix (TOM). Furthermore, genes were clustered into modules by utilizing average linkage hierarchical clustering using topological overlap dissimilarity matrix (1-TOM) as the distance measure, and modules were determined by the dynamic hybrid tree cut algorithm. Finally, similar modules whose eigengenes were highly correlated were merged to trim genes whose correlation with module eigengene was less than the defined threshold (hereby we used 0.25 as a threshold for merging similar modules). WGCNA determines highly inter-connected nodes as modules designated with different colors. Moreover, this algorithm is able to relate the identified modules to external clinical outcomes, and to export the network information to external software for visualization, e.g. VisANT (http://visant.bu.edu/) [44] or Cytoscape [45] that visualizes biological interactions between subnetworks.

Gene ontology analysis and visualization

Gene Ontology (GO) analysis was conducted to reveal biological functions of gene products within significant modules from three scenarios: biological process (BP), cellular component (CC), and molecular function (MF) by using the software named Database for Annotation, Visualization, and Integrated Discovery (DAVID) [46]. DAVID facilitates high throughput gene functional analysis by accepting the input of a list of genes or an individual gene. Potential enriched functions of the DEGs in the determined modules were analyzed by the DAVID tool.

The top 30 high-rank intramedullary hub genes were prepared for visualization with VisANT (Supplementary Table 16, Supplementary Table 17, Supplementary Table 18, Supplementary Table 19, Supplementary Table 20, Supplementary Table 21).

Modeling miRNA–mRNA-transcription factor interaction

To explore transcription factor and miRNA–mRNA interaction in particular modules, the database TargetScan Mouse (Release 7.1) (http://www.targetscan.org/vert_71/) was used for searching the predicted microRNA targets, and MicroCosm Mouse (Release 5) (http://www.ebi.ac.uk/enright-srv/microcosm/htdocs/targets/v5/) was used to build a miRNA–mRNA interaction network and the network was visualized with Cytoscape 3.3.8 software [46]. Putative miRNAs were generated by the intersection of the miRNA-targets from TargetScan and Microcosm. Transcription factors may control different gene modules. Therefore, we further used iRegulon Cytoscape plugin [47] to identify the possible transcription factors and their co-factors that are associated with a collection of intramedullary hub genes extracted from prognosticate modules (Supplementary Table 6, Supplementary Table 7, Supplementary Table 8, 9). To expand our network, we employed CyTargetLinker plugin [48]. CyTargetLinker is able to enhance biological networks using the provided information in frame of Regulatory Interaction Networks or RegIN. RegIN is a network in xgmml format containing regulatory interactions.

Module preservation

To determine the reliability of the identified vascular aging-related modules and compare the modular structure in inferred co-expression network with a reference network generated based on an independent dataset, we used a previous study with GEO access number GSE10000 [17] to identify significantly differentially expressed genes (P < 0.05). In the dataset GSE10000, the two groups of samples corresponding to 32-week old wild-type mice (around 6-month old) (n=3) and 78-week old wild-type mice (around 20-month old) (n=3), were comparable with the first-stage samples corresponding to 6-month old mice (n=6) and 20-month old mice (n=6). Therefore, we removed the samples corresponding to 6-week old mice (n=3) and the ApoE-/- mice. Based on the significantly differentially expressed probe sets, we reconstructed a test co-expression network and contrasted the preservation of co-expression network across testing and reference datasets to detect the conservation of gene pairs between the two networks. Briefly, a Z-summary < 5 and lower median rank indicates weak preservation between the testing and reference networks.

Animal breeding, mouse aorta isolation, and real-time polymerase chain reaction

All animal experiments were approved by the Animal Ethics Committee of Shanghai Jiao Tong University. Male mice (C57BL/6J) were obtained from the Shanghai Laboratory Animal Center, Chinese Academy of Sciences (SLAC, CAS), and kept for 7 days in the local animal house for acclimatization. The mice were 6 and 20 months of age old, housed at a constant ambient temperature under a 12h/12h light-dark cycle and supplied with distilled water, and pelleted AIN-76A (Research Diets, New Brunswick, NJ) chow ad libitum. Quantitative real-time polymerase chain reaction (qPCR) was used to validate the expression of genes changed in aortas from mice aged 6-month (n = 20) and 20-month (n = 20).

A mouse aorta was perfused in situ with cautious to keep the adventitia intact. The vessel was flushed thoroughly with ice-cold phosphate-buffered saline (PBS), through the left ventricle of the heart, cleaned periadventitial fat and connective tissues, snap-frozen in liquid nitrogen, and stored at −80 °C.

Total RNA, including miRNA, was extracted from the tissue using Qiagen RNeasy Mini Kit (Qiagen, Hilden, Germany). Total RNA from aorta was treated with DNase I to remove genomic DNA. RNA integrity was checked by Agilent2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA) and RNA gel electrophoresis.

For the quantification of mRNAs, 1 μg of the total RNA was converted into cDNA using a reverse transcription kit (Promega) based on the manufacturer’s instruction. The gene expression of hub genes (Enpp5, Fez1, Kif1a, F3, and H2-Q7) were then detected by real-time qPCR. The qPCR were conducted using SYBR Premix Ex Taq (Takara, Japan) in 20 μl reaction solution containing 1μl cDNA, 10 μl SYBR Premix Ex Taq (2X), 0.4μl forward primer, 0.4 μl reverse primer, 0.4 μl ROX reference dye, and 7.8 μl ddH₂O. The PCR amplification procedure was carried out on an ABI 7500 fast real-time PCR system (Applied Biosystems, USA) at 95°C for10 s, followed by 40 cycles of 95°C for 5 s and 60°C for 35 s. The amplification reaction without the template was used as a negative control. β-actin gene was used as the internal reference. Primers used for gene amplification are available on request.

The miRNAs (mmu-miR-449a, mmu-miR-449c, mmu-miR-34c, mmu-miR-34b-5p, mmu-miR-15a, and mmu-let-7) were detected by qPCR. Briefly, 2 μg RNA were performed using the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). The qPCR was then conducted using the miRcute miRNA qPCR Detection Kit (SYBR) (Tiangen, Beijing, China) in 20 μl reaction solution containing 1 μl cDNA, 10 μl 2× miRcute miRNA premix, 0.4 μl 50× ROX Reference Dye, 0.4 μl forward primer, 0.4 μl reverse primer, and 7.8 μl ddH₂O. The PCR amplification procedure was carried out on an ABI Prism 7500 Fast Sequence Detection System (ABI, Carlsbad, CA, USA) at 94°C for 2 min, followed by 40 cycles of 95°C for 20 s and 60°C for 35 s. Small nucleolar RNA (RNU6B) was used as the housekeeping loading reference. Forward primers were designed based on mature miRNA sequences while the reverse primers were provided by the miRcute miRNA qPCR Detection Kit (SYBR) (Tiangen, Beijing, China).

The relative gene expression was calculated using the 2^-ΔΔCt method [49] (Song et al., 2018). The experiments were conducted three times independently. Statistical comparison of the levels was analyzed using two-tail unpaired Student’s t-test, and differences were considered significant if P < 0.05.