Weighted Gene Co-Expression Network Analysis Identifies Specific Modules and Hub Genes Related to Hyperlipidemia

Cardiovascular disease (CVD) is becoming the major cause of death around the world and imposes a significant health-economic load [1]. We had required risk status estimation of a standard lipid profile to assess risk severity, just as total cholesterol (TC) [2], triglyceride (TG) [3], low-density lipoprotein cholesterol (LDL-C) [4], apolipoprotein (Apo) B [5], high-density lipoproteins cholesterol (HDL-C) [6], ApoA1 [7] and the ratio of ApoA1 to ApoB [8] are recommended methods for cardiovascular risk prediction. Among these indicators, elevated levels of serum TC and/or TG, which have also been defined as hyperlipidemic, could identify as strong contributors to CVD mortality [9]. Besides, potential risk factors of CVD have been reported as diabetes mellitus, hypertension, tobacco use, male gender, age, and family history of atherosclerotic arterial disease, family history of heart attacks, lack of physical activity, alcohol consumption, obesity, stress and personal genetic profile [10-16]. All these risk factors which have been taken to individual are important genetic components.

Over the last 10 years, genome-wide association studies (GWASes) have revolutionized the discovery of common genetic variants associated with a range of serum lipid and most of these studies were performed by microarray analysis [17]. Gene co-expression network, a widely used method, has been carried out in analyzing microarray data, especially for marking functional modules [18]. Weighted Gene Co-expression Network Analysis (WGCNA), as one of the most useful gene co-expression network, may identify potential therapeutic targets via WGCNA data analysis [19]. we detected the mRNA expression profile of hyperlipidemia samples to find out the highly connected hub genes and significant modules to show the potential molecular mechanisms.

China is a multi-ethnic country and Maonan is one of the 55 minorities with a long-time history. Recent phylogenetic and principal component analyses revealed that the Maonan people belong to the Southeastern Asian group and the genetic relationship between Maonan ethnic group and other minorities in Guangxi was much closer than that between Maonan and Han nationalities. In a previous study, we have found that the incidence of hyperlipidemia was distinct between Maonan and Han ethnic groups [20]. In this study, we performed WGCNA to construct the co-expression network and mark significant modules in the network. The outcomes may help us for further elucidating the innate character of hyperlipidemia, and provide new insights to potential gene biomarkers and signaling pathways to treat hyperlipidemia.

Microarray data of GSE3059 [21] was downloaded from the National Center For Biotechnology Information (NCBI) Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database [22], which was based on the platform of GPL570 Affymetrix Human Genome U133 Plus 2.0 Array. GSE3059 contains 13 peripheral blood samples from hyperlipidemia patients (mean age = 50.15, 7 males and 6 females) and 19 peripheral blood samples from control individuals who had never been diagnosed with a chronic illness (mean age = 40.13, 10 males and 9 females). The CEL files were transformed into the expression value matrix using the Affy package in R [23], and the probe information was then transformed into the gene name using Bioconductor in R [24]. If one gene had more than one probe, the mean expression value of this gene was selected.

WGCNA is a widely used system biology method, which is used to construct a scale-free network from gene expression data [18]. To ensure that the results of network construction were reliable, outlier samples were removed. An appropriate soft threshold power (soft power = 3) was selected in accordance with standard scale-free networks, with which adjacencies between all differential genes were calculated by a power function. Then, the adjacency was transformed into a topological overlap matrix (TOM), and the corresponding dissimilarity (1-TOM) was calculated. Module identification was accomplished with the dynamic tree cut method by hierarchically clustering genes using 1-TOM as the distance measure with a deepsplit value of 2 and a minimum size cutoff of 30 for the resulting dendrogram. Highly similar modules were identified by clustering and then merged together with a height cut-off of 0.4. To test the stability of each identified module, module preservation and quality statistics were computed with the module preservation function implemented in the WGCNA package [25].

The module membership (MM) was defined as the correlation of gene expression profile with module eigengenes (Mes). And the gene significance (GS) measure was defined as (the absolute value of) the correlation between gene and external traits. Genes with highest MM and highest GS in modules of interest were natural candidates for further research [26]. Thus, the intramodular hub genes were chosen by external traits based GS > 0.2 and MM > 0.6 with a threshold of P-value < 0.05 [18]. The gene-gene interaction network was constructed and visualized using Cytoscape software package [27] and Molecular Complex Detection (MCODE) [28] was used to analyze the most notable clustering module. MCODE score > 6 was a threshold for next analysis.

Two groups of study population including 833 unrelated participants of Maonan (248 males, 29.78% and 578 females, 70.22%) and 826 unrelated subjects of Han (252 males, 30.50% and 574 females, 69.5%) were randomly selected from our previous stratified randomized samples [29]. All participants were resided in the Huanjiang Maonan Autonomous County in the Northwestern of Guangxi Zhuang Autonomous Region, which is located in Southwestern China. The participants’ age ranged from 18 to 80 years with a mean age of 55.98 ± 12.78 years in Han and 56.87 ± 12.12 years in Maonan; respectively. The gender ratio and age distribution were matched between the two groups (Table 1). All participants were essentially healthy with no history of coronary artery disease, stroke, diabetes, hyper- or hypo-thyroids, and chronic renal disease. They were free from medications known to affect lipid profiles.

The epidemiological survey was carried out using internationally standardized method, following a common protocol [30]. Information on demographics, socioeconomic status, and lifestyle factors were collected with standardized questionnaires. Cigarette smoking status was categorized into groups of cigarettes per day: ≤ 20 and > 20. Alcohol consumption was categorized into groups of grams of alcohol per day: ≤ 25 and > 25 [31] . Several parameters such as blood pressure, height, weight and waist circumference were measured, while body mass index (BMI, kg/m2) was calculated.

Venous blood samples were obtained from all subjects after at least 12 h of fasting. The levels of serum TC, TG, HDL-C and LDL-C in samples were determined by enzymatic methods with commercially available kits, Tcho-1, TG-LH (RANDOX Laboratories Ltd., Ardmore, Diamond Road, Crumlin Co., Antrim, UK, BT29 4QY), Cholestest HDL, and Cholestest LDL (Daiichi Pure Chemicals Co., Ltd., Tokyo, Japan), respectively. Serum ApoA1 and ApoB levels were detected by the immunoturbidimetric immunoassay (RANDOX Laboratories Ltd.). All determinations were performed with an autoanalyzer (Type 7170A; Hitachi Ltd., Tokyo, Japan) in the Clinical Science Experiment Center of the First Affiliated Hospital, Guangxi Medical University [32].

The normal values of fasting serum TC, TG, HDL-C, LDL-C, ApoA1 and ApoB levels, as well as the ApoA1/ApoB ratio in our Clinical Science Experiment Center were 3.10–5.17, 0.56–1.70, 1.16–1.42, 2.70– 3.10 mmol/L, 1.20–1.60, 0.80–1.05 g/L, and 1.00–2.50; respectively. The individuals with TC > 5.17 mmol/L, and/or TG > 1.70 mmol/L were defined as hyperlipidemic [33].

According to the microarray results, the 5 most dysregulated mRNAs were chosen for further validation by RT-qPCR in hyperlipidemia patients vs. healthy controls. Blood samples were collected and PBMCs were isolated. Total RNAs were isolated from PBMCs using Trizol reagent (Invitrogen), and complementary DNA was synthesized using the TransScript R Frist-Strand cDNA Synthesis SuperMix (Transgen, China) according to the manufacturer’s instructions. Primer sets for selected genes were designed by Sangon Biotech (Shanghai, China); their sequences and reaction conditions are available in Table 2. Each sample was run in triplicate in 96-well plates using LightCycler R 96 and FastStart Essential DNA Green Master (Roche Diagnostics GmbH, Germany). Quantification cycles (Cq) were calculated using the fit point method (LightCycler R 96 Software, Version 1.1 provided by Roche). The expression data were normalized to the reference glyceraldehyde-3-phosphate dehydrogenase (GAPDH). All experiments (sample collection, preparation and storage, primer design, qPCR normalization) were performed according to the MIQE guidelines.

All statistical analyses were performed using the statistical software package SPSS 21.0 (SPSS Inc. Chicago, IL, USA). A chi-square analysis was used to evaluate the difference of the rate between the groups. Continuous data were presented as means ± SD. Comparisons between groups for continuous data were made using the Kruskal-Wallis and Mann-Whitney nonparametric tests. A raw P value of < 0.05 was considered nominally significant. The heart-map of correlation models was measured by R software (version 3.3.0).

We analyzed the GSE3059 and obtained 54560 expression probes separately from each gene expression profile. After data preprocessing, the expression matrices of 21194 genes were obtained from the 13 samples.

We selected a suitable weighted parameter of adjacency function, which is the soft-threshold β, before construction the weighted co-expression network. After the calculation, we selected the correlation coefficient close to 0.8, soft-threshold β = 3 to construct gene modules using the WGCNA package (Fig. 1). After determining the soft threshold, all of genes were used to construct weighted gene co-expression networks. According to the basic idea of WGCNA, we calculated the correlation matrix and adjacency matrix of the gene expression profile of the hyperlipidemia groups, and then transformed them into a topological overlap matrix (TOM) and obtained a system clustering tree of genes on the basis of gene-gene non-ω similarity (Fig. 2). Together with the TOM, we performed the hierarchical average linkage clustering method to identify the gene modules of each gene network (deepsplit = 2, cut height = 0.4) (Fig. 3). In hyperlipidemia groups, a total of six gene modules were recognized by the dynamic tree cut (Fig. 4). Genes that did not belong to any modules were housed in the gray module. The gray gene modules were ignored in this study.

It is a hugely valued biological significance to find out modules most significantly associated with clinical features. The highest association in the Module-feature relationship was found between midnight blue module and glucose (r² = 0.67, P = 0.01), salmon module and WBC (r² = 0.60, P = 0.03) and salmon module and neutrophil (r² = 0.78, P = 0.002), which were selected as module of interest and clinical feature to be studied in subsequent analyses (Fig. 5). The other modules without enough relationship or statistical significance were not further considered. In order to explore biological relevance of midnight blue module and salmon module, 280 genes in midnight blue module and 312 genes in salmon module were respective subjected to Gene Ontology (GO) functional and KEGG pathway enrichment analyses by R cluster Profile package [34]. Biological processes of salmon module were found to focus on neutrophil mediated immunity (P = 4.77 × 10^–7), positive regulation of defense response (P = 1.19 × 10^–6), positive regulation of inflammatory response (P = 4.57 × 10^–6), and regulation of innate immune response (P = 6.26 × 10^–5). At the same time, biological processes of midnight blue module were found to focus on energy reserve metabolic process (P = 4.10×10^–5), and positive regulation of β cell activation (P = 5.67 × 10^–5). However, in KEGG pathway analysis, salmon module was found to focus on leukocyte trans-endothelial migration (P = 3.70 × 10^–5), and neurotrophin signaling pathway (P = 1.43 × 10^–4). At the same time, midnight blue module was found to focus on insulin signaling pathway (P = 3.64 × 10^–6), and insulin resistance (P = 3.64 × 10^–6; Fig. 6). The top 6 of the databases were shown in Table 3 and 4.

When the STRING database [35] was analyzed, a total of 48 nodes and 266 protein pairs were got with a combined weight score > 0.25 in salmon module (Fig. 7A). After analysis in sub-module, only one module with score > 6 was detected by MCODE. As shown in Fig. 7A-1, the top 10 high degree genes, including PACSIN2 (degree = 27), HDAC4 (degree = 24), C11orf21 (degree = 22), RAB35 (degree = 21), UPF1 (degree = 20), NADK (degree = 18), FFAR2 (degree = 17), NRARP (degree = 16), F2RL1 (degree = 11) and ABHD2 (degree = 10) were hub nodes with higher node degrees in this module. However, there were altogether 70 nodes and 232 protein pairs were got with a combined weight score > 0.42 in midnight blue module (Fig. 7B), after analysis in degree, he top 5 high degree genes, including PDZRN4 (degree = 45), TM4SF1 (degree = 42), SCG5 (degree = 36), FAM13A (degree = 34) and GOT1L1 (degree = 33), could be further considered. When combined with GO functional, KEGG pathway enrichment and PPI analysis, we hypothesized that HDAC4, F2RL1, ABHD2, TM4SF1 and FAM13A as the hub genes were closely relevant to hyperlipidemia.

We performed RT-qPCR to validate the hypothesized data. To verify the main conclusion drawn from the microarray results for peripheral blood samples, the relative expression levels of genes coding for HDAC4, F2RL1, ABHD2, TM4SF1 and FAM13A were determined. Overall, the RT-qPCR results were qualitatively consistent with the results of the microarray analysis. However, the RT-qPCR analysis has shown that the F2RL1, ABHD2, TM4SF1 and FAM13A relative expression tended to give higher up-regulation levels in hyperlipidemia group than those in normal controls, besides, F2RL1 and ABHD2 relative expression levels in Maonan ethnic group were higher than those in Han group (Fig. 8A). As shown in Fig. 8B and 8C, the relationship between genes and lipid parameter and HDAC4 was negative correlation with other genes and lipid parameter.

TC, LDL-C, HDL-C, and TG are correlated phenotypes inherited as complex traits and hyperlipidemia as a complex and multifactorial disease caused by an interaction between genetic and environmental factors [36]. In our study, we have identified 5 genes which may modify serum lipid levels, but, the reason for causing the change of serum lipid levels was slightly different.

Histone deacetylases (HDACs) are enzymes that regulate gene expression by removing acetyl residues from target histones or non-histone proteins and consequently modifying chromatin structure. HDACs are divided into four classes (class I: HDAC1, 2, 3 and 8; class II: HDAC4, 5, 6, 7, 9 and 10; class III: SIRT1-7; class IV: HDAC11) [37]. HDAC4, one member of the tissue-specific Class II HDACs, is highly expressed in neurons [38] and bone mass, and plays an essential role in maintaining neuronal survival [39] and chondrocyte hypertrophy [40]. When we analysis the other 4 genes in UCSC, we could only find that the F2RL1 and ABHD2’ concentrate peak appeared in H3K27AC before gene transaction, but nothing can be found in TM4SF1 and FAM13A (Fig. 9). Braun et al. found that DNA methylation was a key epigenetic mechanism that is suggested to be associated with blood lipid levels [41]. Maybe histone modification is another key epigenetic mechanism for lipid modified.

Besides, nuclear HDAC4 distribution was enhanced in Purkinje neurons from Atm-deficient mice after lipopolysaccharides (LPS) stimulation, and Atm was identified to be involved in ataxia-telangiectasia characterized by immune deficiency [42], indicating that HDAC4 may directly or indirectly regulate inflammation genes. Recently, novel data have emerged that immune responses and lipid metabolism interact in a unique metabolic pathway underlying atherosclerosis [43]. Mihaylova et al. have suggested that class I/II HDACs may be influenced blood glucose and would be potential therapeutics for metabolic syndrome [44]. In our current study, F2RL1, ABHD2, TM4SF1 and FAM13A expression tended to give higher up-regulation levels in hyperlipidemia group than those in normal controls and HDAC4 was negatively correlated with other genes. These results are consistent with the previous findings.

Except for genetic factors, environmental exposures were another aspect which cannot be ignored. Maonan people like to pickle sour meat, snails and vegetables. They get meat mainly from poultries and livestocks, such as pigs, oxen, chickens, ducks and so on. A typical food, Minglun Sliced Pig is a well-known dish of the Maonan ethnic group. Most of the Maonan people like to eat food which is cooked half ripe. In addition, they also like to eat beef, pork and /or animal offals in a hot pot which contain abundant saturated fatty acid [20]. Long-term high saturated fat diet is an important risk factor for obesity, dyslipidemia, atherosclerosis, and hypertension [45]. The major dietary saturated long-chain fatty acids such as myristic acid (14: 0) and palmitic acid (16: 0) have been associated with deleterious effects on blood lipid metabolism, especially due to their influence on plasma TC and TG levels [46].

Recently, HDAC inhibitors in turn regulate the activity of HDACs, and have been widely used as therapeutics in psychiatry and neurology, in which a number of adverse outcomes are associated with aberrant HDAC function. At the same time, dietary HDAC inhibitors have been shown to have a regulatory effect similar to that of pharmacological HDAC inhibitors without the possible side-effects [47]. Wood and Yang et al. have found that a small number of dietary compounds appear to have dual effects on histone acetylation status and excessive alcohol consumption and cigarette smoking have been associated with histone acetylation status modified [48, 49]. In the present study, we showed that the percent of alcohol consumption were higher in Maonan than in Han groups. Most of the local adult men of the Maonan people like to drink. They even have the custom that it would be considered to be impolite to treat their guests without wines. Some families made wines themselves using grain sorghums and corns. Excessive alcohol consumption, just as HDAC inhibitors, removal of acetyl moiety alters the histone configuration, once more returning the chromatin to the closed form and influence transaction. That would be an explanation for the reason why the relative expression of F2RL1 and ABHD2 were higher in Maonan than in Han groups.

In conclusion, thirteen hyperlipidemia microarray datasets from GEO series had been systematic analyzed. After finished weighted gene co-expression networks analysis, we found salmon and midnight blue module with the highest correlation. We performed GO functional, KEGG pathway enrichment and protein-protein interactions analyses for these two modules, five genes (HDAC4, F2RL1, ABHD2, TM4SF1 and FAM13A) were found with a significant meaning. After validation with RT-qPCR, we found that the relative expression of HDAC4 was lower, but F2RL1, ABHD2, TM4SF1 and FAM13A were higher in hyperlipidemia than in normal controls. Subgroup analysis showed that the relative expression of HDAC4 was lower, but, F2RL1 and ABHD2 expression was higher in Maonan than in Han ethnic groups. Except for genetic factors and environmental exposures, epigenetic influence was another mechanism which cannot be ignored, but, it was still needed to further confirm.

L.M. conceived the study, participated in the design, undertook genotyping, performed the statistical analyses, and drafted the manuscript. R.-X.Y. conceived the study, participated in the design, carried out the epidemiological survey, collected the samples, and helped to draft the manuscript. S.Y. collaborated to the genotyping. S.-L.P., D.-Z.Y. and W.-X.L. carried out the epidemiological survey and collected the samples. All authors read and approved the final manuscript. The authors acknowledge the essential role of the funding of the National Natural Science Foundation of China (No: 81460169) and the Innovation Project of Guangxi Graduate Education in this motif.

The authors have no potential conflicts of interest to report.