- Open Access
Bioinformatics analysis of vascular RNA-seq data revealed hub genes and pathways in a novel Tibetan minipig atherosclerosis model induced by a high fat/cholesterol diet
Lipids in Health and Disease volume 19, Article number: 54 (2020)
Atherosclerosis is a major contributor to cardiovascular events, however, its molecular mechanism remains poorly known. Animal models of atherosclerosis can be a valuable tool to provide insights into the etiology, pathophysiology, and complications of atherosclerosis. In particular, Tibetan minipigs are a feasible model for studying diet-related metabolic and atherosclerotic diseases.
We used vascular transcriptomics to identify differentially expressed genes (DEGs) in high fat/cholesterol (HFC) diet-fed Tibetan minipig atherosclerosis models, analyzed the DEGs gene ontology (GO) terms, pathways and protein-protein interactions (PPI) networks, and identified hub genes and key modules using molecular complex detection (MCODE), Centiscape and CytoHubba plugin. The identified genes were validated using the human carotid atherosclerosis database (GSEA 43292) and RT-PCR methods.
Our results showed that minipigs displayed obvious dyslipidemia, oxidative stress, inflammatory response, atherosclerotic plaques, as well as increased low-density lipoprotein (LDL) and leukocyte recruitment after 24 weeks of HFC diet feeding compared to those under a regular diet. Our RNA-seq results revealed 1716 DEGs in the atherosclerotic/NC group, of which 1468 genes were up-regulated and 248 genes were down-regulated. Functional enrichment analysis of DEGs showed that the HFC diet-induced changes are related to vascular immune-inflammatory responses, lipid metabolism and muscle contraction, indicating that hypercholesterolemia caused by HFC diet can activate innate and adaptive immune responses to drive atherosclerosis development. Furthermore, we identified four modules from the major PPI network, which are implicated in cell chemotaxis, myeloid leukocyte activation, cytokine production, and lymphocyte activation. Fifteen hub genes were discovered, including TNF, PTPRC, ITGB2, ITGAM, VCAM1, CXCR4, TYROBP, TLR4, LCP2, C5AR1, CD86, MMP9, PTPN6, C3, and CXCL10, as well as two transcription factors (TF), i.e. NF-ĸB1 and SPI1. These results are consistent with the expression patterns in human carotid plaque and were validated by RT-PCR.
The identified DEGs and their enriched pathways provide references for the development and progression mechanism of Tibetan minipig atherosclerosis model induced by the HFC diet.
Hypercholesterolemia and atherosclerosis are major lifestyle-related metabolic disorders and health problems in modern society . Their occurrences are known to be related to genetic and environmental factors [2, 3]. In specific, hypercholesterolemia is an important risk factor for atherosclerosis, which can lead to atherosclerosis by inducing a local inflammatory response, oxidative stress damage, and damages to vascular endothelial cell metabolism and function. Recently, studies have suggested that imbalances in nutrient intake may lead to hypercholesterolemia and atherosclerosis [4, 5]. However, a lack of appropriate animal models has made it difficult to assess such an association between hypercholesterolemia and atherosclerosis caused by dietary patterns. Rodents are not an appropriate model for atherosclerosis in humans, because they have such a different lipoprotein profile that atherosclerosis rarely occurs in rodents unless they are genetically modified or fed to an extreme diet . Moreover, in these models, genetically induced hypercholesterolemia usually leads to cholesterol levels that are much higher than those in humans. In particular, mouse models have been challenged recently for studying human inflammatory/immune responses . Therefore, the establishment of an ideal experimental animal disease model is essential for studying the pathophysiology of atherosclerosis.
Pigs, which can spontaneously form atherosclerosis, provide a more representative animal model for studying dyslipidemia and atherosclerosis [7, 8]. The physiological and anatomical aspects of pigs in lipid metabolism, cardiovascular metabolism and digestion are very similar to those in humans. Furthermore, the omnivorous nature of pigs makes them highly suitable for studying dietary-based patterns and for translating findings from pigs to humans . In particular, Tibetan minipigs have many advantages as a model for studying diet-related cardiovascular metabolic diseases [9, 10]. Our previous studies have shown that Tibetan minipigs are a viable model for studying atherosclerosis induced by dietary patterns. They can form atherosclerosis with 24 weeks of a high-fat/cholesterol diet, and have characteristics like lipid disorders, insulin resistance, hypertension and cardiac dysfunction . The lipid metabolism changes and atherosclerotic lesions in the Tibetan minipig model are also very close to those in humans, but the underlying molecular pathological mechanisms are not yet clear.
Recently, it has been recognized that the relationship between dietary patterns and cardiovascular diseases is complex, which involves many changes in biological systems . The increase in genomic information, especially after the emergence of high-throughput sequencing technologies, provides an opportunity to better understand the intricate biological processes of the disease. Previous studies have revealed molecular changes in rabbit and human atherosclerosis at the mRNA level [12, 13]. Recently, hypercholesterolemia has been found to stimulate immune responses and may promote inflammatory diseases [14, 15]. Herein, we utilized Tibetan minipigs fed with a high-fat/cholesterol diet to mimic human hypercholesterolemia and atherosclerosis associated with today’s poor dietary patterns. We used vascular transcriptome to identify the extent of changes in vascular gene expression, and performed bioinformatics analysis to discover potential hub genes and key pathways caused by atherosclerosis. This study therefore provided clues about the molecular mechanism of atherosclerosis, revealed new therapeutic targets and biomarkers, and can serve as a reference for future research and application of Tibetan minipigs in studying human atherosclerotic diseases.
Materials and methods
Animals and experimental protocol
Before the study, the experimental protocol was approved by the Animal Care and Use Committee at the Zhejiang Chinese Medical University (approval number: ZSLL-2016-0031). Ten male Tibetan minipigs (3–4 months old) were purchased from Guangdong Pearl Lab Animal Sci & Tech Co., Ltd. All minipigs were kept in the animal house of Experimental Animal Research Center of Zhejiang Chinese Medical University at 20–22 °C and with a 12 h light/dark cycle. After 4 weeks of environmental accommodation, the animals were divided into a normal control (NC) group and an atherosclerosis model group with 5 animals per group. The atherosclerosis model was established by feeding minipigs with a high-fat/cholesterol (HFC) diet (containing 15% shortening oil, 1.5% cholesterol, 10% egg yolk powder, and 73.5% basal feed) for 24 weeks, as we previously described . Minipigs were euthanized at the 24th week.
Changes in biochemical parameters and white blood cell (WBC) count
Following a 16 h fast at the 24th week of a HFC diet, the body weights of minipigs were measured. Heparin anticoagulant blood was collected from the anterior vena cava of minipigs, and was then centrifuged for plasma collection. Total cholesterol (TC), triglycerides (TG), high-density lipoprotein-cholesterol (HDL-C), and low-density lipoprotein-cholesterol (LDL-C) were measured (Hitachi 7020 automatic biochemical analyzer) as previously described . Atherogenic index (AI) was derived from the ratio of (TC–HDL-C) to HDL-C. Plasma C-reactive protein (CRP) and tumor necrosis factor-α (TNF-α) levels were detected using ELISA kits (Nanjing Jiancheng Bioengineering Institute, China). In addition, 1 mL of EDTA-K2 anticoagulant blood was taken, and the components of the blood cells including the WBC count were analyzed on an automated hematology analyzer (2120, Bayer Healthcare Ltd., USA).
Atherosclerosis lesion analysis
At the end of the study, all minipigs were sacrificed for autopsy. The blood vessels from the aortic arch to the abdomen of the descending aorta (flat 3–4 lumbar vertebrae) were fixed with 10% buffered formalin, stained with 0.5% Sudan IV, and then photographed. Image J software was used to quantify Sudan IV-stained area, expressed as a percentage of the total surface area. The arterial vessels at the end of the abdominal aorta were placed into 10% paraformaldehyde solution, embedded in paraffin, cut into 5-μm-thick sections, stained with hematoxylin-eosin (H&E) and Masson’s trichrome, and immunohistochemically stained with NF-ĸB p65 subunit, VCAM1, and MMP9 antibodies (1:300 dilution, NF-ĸB p65,SC-8008; VCAM1, SC-18854; MMP9, SC-21733; Santa Cruz Biotechnology, CA, USA). Image pro plus 6.0 software was used to measure the intima-media thickness (IMT) and the staining percentage of the positive area . In addition, immunofluorescent staining was performed with smooth muscle actin antibody (α-SMA) and CD68 antibody (1:200; α-SMA, SC-32251; CD68, SC-70761; Santa Cruz Biotechnology, CA, USA), which was then secondarily labeled with Alexa fluor 647-conjugated goat anti-rabbit IgG (ab150087) and Alexa fluor 488-conjugated goat anti-mouse IgG (ab150117) (Abcam, Cambridge, UK) respectively, and mounted with Fluorescent with DAPI. Slides were examined with a fluorescence microscope (VS120, OLYMPUS, Japan).
Construction and sequencing of messenger RNA (mRNA) library
We randomly picked three pigs in every group for deep sequence analysis. Trizol reagent was used for RNA isolation from blood vessels at the end of abdominal aorta (Invitrogen, CA, USA). RNA quantity and purity were assessed using a Bioanalyzer 2100 and RNA 6000 Nano Lab Chip Kit (Agilent, CA, USA) with RNA integrity number (RIN) > 7.0. Ten micrograms of total RNA was used to isolate Poly (A) mRNA with poly-T oligo attached magnetic beads (Invitrogen). After purification, the mRNA was cleaved into small fragments by divalent cation at high temperatures. According to the protocol of the mRNA Seq sample preparation kit (Illunina, San Diego, USA), the RNA fragments after cleavage were reverse-transcribed to generate the final cDNA library. The average insertion size of the paired-end library was 300 bps. We then performed the paired-end sequencing on Illumina Hiseq 4000 (LC Sciences, Hangzhou, China) as recommended by the supplier. RNA-Seq data have been submitted to the Gene Expression Omnibus (GEO) database (Accession no.GSE140412).
Transcriptome data analysis
Clean reads were obtained and mapped to the pig reference genome (Scrofa 11.1) using the HISAT package (version 2.0) with default parameters. Each sample’s mapped reads were assembled using StringTie (version 1.3.0). All the transcriptomes in the sample were merged to reconstruct a comprehensive transcriptome using perl scripts. The average FPKM (Reads per kilobase of exon per million) were calculated using StringTie. Differentially expressed genes (DEG) were selected with log2 (fold change) > 1 or log2 (fold change) < − 1 and with statistical significance (p-value < 0.05) using edgeR.
Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis of DEGs
To assess the functions of the identified DEGs, g: profiler (https://biit.cs.ut.ee/gprofiler) was used to identify the enriched GO biological processes . The results were integrated using the Enrichment Map , in which p-value < 0.01 was adjusted to indicate that the difference is statistically significant. KEGG Orthology Based Annotation System (KOBAS 3.0) (http://kobas.cbi.pku.edu.cn) was used to identify the enriched KEGG pathway based on adjusted P values, in which p-value < 0.05 was adjusted to indicate that the difference is statistically significant. In addition, Gene Set Enrichment Analysis (GSEA, version 4.0.3, the broad institute of MIT and Harvard) was used to assess whether a genetically defined genome is statistically significant between the two phenotypes . GSEA was performed to identify enriched hallmark gene sets pathways between atherosclerosis and NC (FDR < 0.05), where the hallmark gene sets were collected from the Molecular Signature Database (MSigDB).
Identification of modules and hub genes in PPI networks
DEGs were uploaded to the STRING database to build a PPI network , with a confidence score ≥ 0.4 as the cut-off criterion. The PPI network was then visualized in Cytoscape software 3.7.2 . The PPI network centrality parameters (degree, betweenness, eigenvector) was analyzed using CentiScaPe 2.2 similar to previous studies [21, 22]. In this work, the average values of degree, betweenness, and eigenvector were calculated. To identify nodes with higher centrality values for the three selected parameters, a Venn diagram was drawn using Venny 2.1 software (https://bioinfogp.cnb.csic.es/tools/venny/). Finally, the major PPI network was constructed based on the intersection of three centrality parameters in the Venn diagram. The most significant modules were identified from the major PPI network using the Molecular Complex Detection (MCODE) plugin. K-score = 5 was set as the criterion for MCODE analysis. In addition, we used the CytoHubba plugin and the degree topological analysis  to identify hub genes, where genes of the top 15 degrees were considered to be hub genes in the network. Transcription factor prediction was performed using the iRegulon plugin of Cytoscape 3.7.2 .
Hub genes validation
We validated the RNA-seq results of 9 genes by real-time PCR. Briefly, complementary cDNA was synthesized using a PrimeScript™ RT reagent kit (Takara, Dalian, China). Table 1 shows the sequences of the used primers. Real-time PCR was performed on Step one plus using SYBR Premix Ex Taq Reagent Kit (Takara, Dalian, China). GAPDH was used as a house-keeping reference to calculate the relative expression of target genes with the 2-ΔΔCt method. To confirm the expression patterns of the15 hub genes in human atherosclerosis, we used the GEO database of human carotid atherosclerotic plaques (GSE43292, containing 32 atheroma plaque samples and 32 adjacent carotid tissue) for validation.
All data are expressed as means ± SEM. Statistical analysis was performed with GraphPad Prism 7.0 software (GraphPad Software Inc., USA). The Mann-Whitney test was used to compare between the two groups. P < 0.05 indicates statistical significance.
Changes of clinical phenotypic characteristics and histopathological assessment in Tibetan minipigs with a HFC diet
We saw from the H&E staining that the abdominal aortic intima was thickened in the atherosclerosis group, with a large number of inflammatory cells infiltrating, calcifying and forming foam cells (Fig. 1b). In specific, our quantitative analysis revealed a significant increase of IMT increased in the atherosclerosis group (P < 0.01, Fig. 1c). Masson staining showed that the content of collagen fibers in the intima and the media grew significantly after 24 weeks of HFC diet induction (Fig. 1d). From Sudan IV staining, we saw a significant increase of aortic lipid deposition in the atherosclerosis group compared with the NC group (P < 0.01, Fig. 1e-f). Meanwhile, significant changes were also observed in clinical phenotypic characteristics, as shown in Table 2. The body weight, TC, LDL-C, HDL-C, AI, WBC counts, as well as CRP and TNF-α levels, were significantly increased in the atherosclerosis group compared to the NC group (P < 0.05, P < 0.01), except for TG.
Identification of differential genes and analysis of their functions and pathways
To observe the changes of vascular transcriptome in the Tibetan minipig atherosclerosis model, six cDNA libraries were constructed from high-throughput sequencing of three abdominal aortic samples from each group (N = 3). One thousand seven hundred sixteen DEGs were identified in the atherosclerosis/NC group, including 1468 up-regulated genes and 248 down-regulated genes (Fig. 2a-c).
To understand the biological function of atherosclerosis in Tibetan minipigs induced by the HFC diet, the g: profiler online tool was used to identify overrepresented biological pathways, and the Enrichment Map plugin was used for enrichment visualization (Fig. 2d). These results showed that the up-regulated DEGs are markedly enriched in many biological processes including immune-inflammatory response, cell proliferation activation, cell chemotaxis, metal ion calcium, cell differentiation, response to lipid, bone remodeling, reactive oxygen species metabolic process, and cell death, whereas the down-regulated DEGs are enriched in muscle contraction. In KEGG analysis (shown in Fig. 2e), we saw that the up-regulated DEGs are enriched in cytokine-cytokine receptor interaction, phagosome, lysosome, hematopoietic cell lineage, cell adhesion molecules (CAMs), natural killer cell mediated cytotoxicity, antigen processing and presentation, chemokine signaling pathway, complement and coagulation cascades, B cell receptor signaling pathway, NF-kappa B signaling pathway, T cell receptor signaling pathway, Fc gamma R-mediated phagocytosis, Toll-like receptor signaling pathway, TNF signaling pathway, PPAR signaling pathway, JAK-STAT signaling pathway and insulin resistance, while the down-regulated DEGs are enriched in steroid biosynthesis, vascular smooth muscle contraction, metabolic pathway, calcium signaling pathway, cGMP-PKG signaling pathway. In addition, gene set enrichment analysis (GSEA) method was used to identify the gene set hallmarks. The results showed that there were 17 pathways activated in the atherosclerosis/NC gene set, including interferon alpha response, interferon gamma response, inflammatory response, complement, coagulation, TNFα signaling via NFKB, reactive oxygen species pathway and apoptosis, while cholesterol homeostasis and myogenesis were significantly inhibited (Fig. 2f and g).
Hub genes and key network identification in the PPI network
To investigate the interaction of these 1716 DEGs in the atherosclerosis/NC gene set, an initial PPI network was created using STRING 11.0. The threshold of interaction score was set to > 0.4 (medium confidence), and the unconnected nodes were removed from the network. The initial PPI network consisted of 7137 edges and 1236 nodes (Fig. 3).
To obtain the major PPI network from the initial network, we performed the following analysis. First, each node’s topological connectivity was determined based on the degree, betweenness, and eigenvector measurement, as shown in Fig. 4a. Three hundred eighty-eight nodes were identified from the PPI network with a degree value above the average (−x = 11.549), and were therefore considered as hub nodes. Second, the betweenness of nodes exceeding the average value (−x = 3400.411) were considered as bottlenecks, which are the key connector proteins in the PPI network. In total, 360 bottlenecks nodes were identified. Finally, through further analysis of eigenvector parameters, a smaller and more important network was found from the PPI network. There were 297 nodes with an eigenvector value higher than the average (x = 0.013). According to the analysis of centrality parameters, a total of 195 nodes showed high centrality values, which is equivalent to 39.2% of the initial PPI network.
The resulting major PPI network consists of 195 nodes and 2147 edges. To analyze the major PPI network, we identified the top 15 hub genes with higher degree connectivity using the CytoHubba plugin (Fig. 4c), and sequentially ordered them as follows: TNF (Degree = 75), PTPRC (Degree = 71), ITGB2 (Degree = 58), ITGAM (Degree = 57), VCAM1 (Degree = 53, CXCR4 (Degree = 52), TYROBP (Degree = 52), TLR4 (Degree = 50), LCP2 (Degree = 50), C5AR1 (Degree = 48), CD86 (Degree = 47), MMP9 (Degree = 47), PTPN6 (Degree = 46), C3 (Degree = 46), CXCL10 (Degree = 45). The names and functions of these hub genes are shown in Table 3. The iRegulon plugin application in Cytoscape software was used to analyze the promoter binding motifs and to identify the transcription factors associated with each gene. The results show that SPI1 or NFĸB-1 are important regulators of the top 15 hub genes. These 15 hub genes mainly participate in the biological process of myeloid leukocyte activation (Fig. 5e), and are involved in KEGG pathways like natural killer cell mediated cytotoxicity, leukocyte transendothelial migration, cell adhesion molecules (CAMs), complement and coagulation cascades, Toll-like receptor signaling pathway, T cell receptor signaling pathway, TNF signaling pathway, phagosome, NF-kappa B signaling pathway, and osteoclast differentiation (Fig. 5e). Moreover, four significant modules were found via the MCODE plugin (Fig. 4b), with CXCR4, C5AR1, C3 and CXCL10 in module 1, ITGB2, ITGAM, VCAM1, TYROBP, LCP2 and PTPN6 in module 2, TNF, TLR4 and CD86 in module 3, and MMP9 in module 4. Enrichment GO biological processes and pathways of the four modules are displayed in Fig. 5. In these four modules, the most significantly enriched GO biological processes include cell chemotaxis, myeloid leukocyte activation, cytokine production and lymphocyte activation, while the main enriched KEGG pathways include chemokine signaling pathways, insulin resistance, leukocyte transendothelial migration, phagosome, complement and coagulation cascades, osteoclast differentiation, hematopoietic cell lineage, natural killer cell mediated cytotoxicity, cell adhesion molecules (CAMs), cytokine-cytokine receptor interaction, NOD-like receptor signaling pathway, TNF signaling pathway, adipocytokine signaling pathway, Toll-like receptor signaling pathway, lysosome, endocytosis, VEGF signaling pathway, Fc gamma R-mediated phagocytosis, NF-kappa B signaling pathway, and HIF-1 signaling pathway.
Hub genes validation
To verify the results of RNA-seq data, we selected 9 genes (8 up-regulated and 1 down-regulated) for real-time PCR analysis. We found that the mRNA levels of TNF, ITGB2, VCAM1, TYROBP, LCP2, SPI1, CTSS and PTPN6 in the atherosclerosis group were significantly increased (P < 0.05, P < 0.01) compared to the NC group, while the mRNA level of SQLE was significantly decreased (P < 0.01, Fig. 6a). The RT-PCR results are completely consistent with our RNA-seq results. Immunohistochemical analysis showed positive staining and significantly increased expression of NF-ĸB, VCAM1 and MMP9 proteins in the atherosclerosis group compared to the NC group (P < 0.05, P < 0.01, Fig. 6b and c). Immunofluorescence staining also showed that the expression of vascular SMCs and CD68 were significantly increased in the atherosclerosis group (P < 0.05, P < 0.01; Fig. 6d). Furthermore, to confirm whether these 15 hub genes have a similar expression pattern in human atherosclerotic diseases, we used the human carotid plaques (GSE43292) dataset for validation. As shown in Fig. 6e, the mRNA expression levels of the top 15 hub genes are increased in the human atheroma plaque compared to those in the adjacent carotid tissues.
Poor diet quality is a key factor leading to hypercholesterolemia and atherosclerosis. Current guidelines for reducing the risk of cardiovascular diseases also focused on dietary habits. Therefore, animal models that mimic the pathogenesis of human atherosclerosis have become increasingly important in preclinical studies for assessing dietary and drug interventions, as well as for elucidating their underlying mechanism. Our results showed that Tibetan minipigs developed dyslipidemia and immune inflammation accompanied by early-stage atherosclerosis after the consumption of a HFC diet relative to a regular diet. The main pathological changes are characterized by lipid deposition, hyperplasia of foam cells, smooth muscle cells and fiber-matrix, as well as calcification, which is similar to the pathogenesis of human atherosclerosis lesions. This provides a reliable model for studying the relationship between diet and atherosclerosis. Therefore, the pathogenesis of the Tibetan minipig atherosclerosis model deserves further study.
Atherosclerosis is a slowly progressive disease accompanied by inflammation, and inflammatory and immune responses are important links in its formation. High-throughput sequencing can provide a new way for studying the pathogenesis of atherosclerosis. In this study, 1716 DEGs were identified that include 1468 up-regulated DEGs and 248 down-regulated DEGs. GO biological process enrichment, KEGG pathway, PPI network, and module analysis were helpful to find key genes and pathways in Tibetan minipigs with atherosclerosis. GO biological process enrichment analysis showed that up-regulated DEGs are mainly related to immune-inflammatory responses, cell proliferation, differentiation and chemotaxis, response to lipids, reactive oxygen species metabolic process and cell death, which confirmed that cell proliferation and differentiation, cell adhesion, immune inflammation, and lipid metabolism are well associated with atherosclerosis . Meanwhile, our analysis showed that the down-regulated DEGs are involved in muscle contraction, which is in agreement with the results in human carotid atherosclerotic plaque from Liu et al. . These results supported that lipid metabolism, oxidative stress and immune-inflammatory responses drive the development of atherosclerosis, which is consistent with the enrichment analysis results of hallmark gene sets in GSEA.
After the HFC diet, the LDL-C level of Tibetan minipigs was significantly increased with increased leukocyte recruitment. It can be seen that the increase of lipids (especially LDL) and leukocytes are key causes of the chronic inflammation of atherosclerosis. There was increasing evidence that LDL stimulates innate and adaptive immunity in atherosclerosis [27, 28]. Similar results were also confirmed in KEGG pathways enrichment analysis, where phagosomes, lysosomes, cytokine-cytokine receptor interactions, chemokine signaling pathway, complement and coagulation cascades, and antigen processing and presentation were shown to participate in the activation of the adaptive immune response by innate immune cells . At the same time, T and B cell receptor signaling pathways associated with the adaptive immune responses were also enriched and involved in atherosclerosis and arterial thrombosis in humans . Natural killer cell-mediated cytotoxicity was also known to play a role in both innate and adaptive immune-mediated responses . In addition, the Toll-like receptor signaling pathway, NF-kappa B signaling pathway and JAK-STAT signaling pathway have been confirmed as atherosclerosis-related inflammation signaling pathways . Therefore, we believe that a HFC diet induces hypercholesterolemia in Tibetan minipigs and activates innate and adaptive immunity to drive the development of atherosclerosis. Interestingly, it was also found that insulin resistance pathways are significantly enriched in blood vessels, which is consistent with our early study showing that the HFC diet induces insulin resistance and promotes vascular insulin resistance in Tibetan minipigs .
To further understand the molecular mechanism of HFC diet-induced atherosclerosis in Tibetan minipigs, a representative PPI network of atherosclerosis was obtained. We first constructed an initial PPI network of 1716 DEGs. Based on the degree, betweenness and eigenvector analysis of the initial PPI network topology, 195 DEGs were selected to create a high centrality major PPI network, indicating the decisive roles of these proteins in atherosclerosis induced by HFC-diet in Tibetan minipigs. Four key modules were obtained by MCODE cluster analysis, in which we saw that the genes in these modules are mainly involved in the biological process of cell chemotaxis, myeloid leukocyte activation, cytokine production, and lymphocyte activation. Moreover, a total of 15 hub genes were extracted from the major PPI network through degree centrality analysis, including TNF, PTPRC, ITGB2, ITGAM, VCAM1, CXCR4, TYROBP, TLR4, LCP2, C5AR1, CD86, MMP9, PTPN6, C3, and CXCL10. GO biological process enrichment analysis revealed that these hub genes are mainly involved in myeloid leukocyte activation. Previous studies have shown that a variety of myeloid cells (monocytes, macrophages, dendritic cells, and neutrophils) in the aorta are involved in regulating the progression of inflammation and atherosclerosis , and suggested that the activation of myeloid leukocytes can promote vascular endothelial injury and inflammatory response, and accelerate the progression of atherosclerosis. For example, monocytes can differentiate into macrophages and swallow lipids, form foam cells, and secrete a variety of inflammatory factors such as TNF-α, interleukin-1 (IL-1) and IL-6, leading to cell adhesion, inflammatory cell infiltration, matrix degradation, and ultimately plaque rupture . Among the hub genes, TNF, PTPRC (also known as CD45), VCAM1, CD86 and MMP9 have been confirmed as inflammatory biomarkers of atherosclerosis [33, 34]. Chemokines are known to play a crucial role in atherosclerosis, especially in the recruitment of monocytes. CXCR4 is widely expressed not only on the surface of hematopoietic stem cells but also in atherosclerotic lesions. Animal experiments have shown that CXCR4 expression was significantly up-regulated in the thoracic aorta vessels of high-fat diet-fed ApoE−/− mice . Heller EA et al.  showed that lymphocytes indeed promote atherosclerosis, and identified CXCL10, a T cell chemokine which enhances this process. Complement and coagulation cascades are the main defense systems that mediate inflammation and thrombosis. Hypercholesterolemia can activate complement to produce C3 invertase through the classical pathway, lectin pathway or alternative pathway. At the same time, the downstream of C5 invertase was formed to lyse C5 and form C5a . Patzelt J et al.  found that the expression of C3aR and C5aR in patients with coronary heart diseases was positively correlated with platelet activation. Recently, NK cells, especially NKT cells, have been identified as important participants in immune metabolism because of their unique responses to lipid antigens and their mixed nature of innate and adaptive immune systems . ITGB2 (integrin subunit beta 2, CD18) is an adhesion molecule of the integrin family and a major marker of NK cell maturity, and is involved in adhesion between inflammatory cells and endothelial cells, inflammatory cells chemotaxis, and other processes that are active in the early stages of atherosclerosis . ITGAM encodes the integrin αM chain that is involved in CD40L-mediated inflammation during atherosclerosis . TYROBP (also known as DNAX activator protein 12, DAP12) is a kind of transmembrane receptors that is widely present in natural killer cells, neutrophils and monocytes/macrophages. Wang et al.  found that DAP12 was highly exhibited in plaques of APOE mice and promoted the formation of atherosclerosis through the TREM-1/DAP12 pathway. Type 6 protein tyrosine phosphatase non-receptor (PTPN6, also named as SHP-1) is a member of the protein tyrosine phosphatase (PTP) family located in the cytoplasm, which regulates cell growth, differentiation, mitotic cycle and oncogenic transformation . PTPN6 plays a vital role in negatively regulating insulin action and liver clearance, thereby regulating systemic glucose homeostasis . Nai et al.  found that PTPN6 was significantly up-regulated in atherosclerosis. LCP2, also known as SLP-76, is an adaptor protein specific to blood-derived cells, which is mainly present in the hematopoietic system and platelets, NK cells, neutrophils, mast cells, immature B cells and macrophages, and participates in T cell receptor-related signaling pathways and coagulation function . It can be seen that the activation of these hub genes is involved in promoting the adhesion between inflammatory cells and endothelial cells, and in accelerating the formation of atherosclerotic plaques.
In addition, TLR4, a member of pattern recognition receptors family, plays a key role in innate immunity and has recently received widespread attention. In monocytes, TLR4 and TNF-α can activate the classic NF-κB/AP-1 signaling pathway and promote the release of inflammatory cytokines . This is consistent with our findings about the enrichment of TNF, TLR, and NF-kappa B pathways. We found that 11 hub genes are regulated by the NF-ĸB-1 transcription factor, which was validated by the high expression of NF-ĸB in atherosclerotic vessels in Tibetan minipigs from immunohistochemical staining. On the other hand, 12 hub genes are controlled by the SPI1 transcription factor (spleen focus forming virus proviral integration oncogene), which is a member of the ETS transcription factor family. Previous studies have shown that the mRNA level of SPI1 is up-regulated during myeloid cell differentiation and is maintained at high levels in human B cells, mast cells, monocytes, peripheral blood neutrophils and other granulocytes . This is consistent with our observation of the significant increase of SPI1 expression in Tibetan minipigs induced by the HFC diet, further suggesting that the activation of myeloid cells plays a very important role in the formation of atherosclerosis.
To ensure the accuracy of candidate DEGs, we selected 9 genes for RT-PCR validation. The results confirmed the consistency between the gene expression and sequencing results. Furthermore, we used the human carotid atherosclerosis dataset (GSE 43292) to validate the expression pattern of 15 hub genes in human carotid plaques. We found that the expression patterns of these 15 hub genes in human carotid artery plaques agree with those in atherosclerotic vessels of Tibetan minipigs induced by the HFC diet, suggesting that the HFC diet-induced Tibetan minipigs atherosclerosis model is suitable for the study of atherosclerosis-related diseases in humans.
Our data showed that the HFC diet-induced Tibetan minipig atherosclerosis model exhibits significant dyslipidemia, oxidative stress, immune-inflammatory responses, and increased LDL and leukocyte recruitment. Our RNA-seq experiment revealed 1716 DEGs in the atherosclerosis/NC group, including 1468 up-regulated genes and 248 down-regulated genes. Functional enrichment analysis of DEGs showed that the HFC diet induced changes are related to vascular immune-inflammatory responses, lipid metabolism and muscle contraction, indicating that HFC diet-induced hypercholesterolemia activates innate and adaptive immunity to drive atherosclerosis development. Four modules, associated with cell chemotaxis, myeloid leukocyte activation, cytokine production, and lymphocyte activation, were identified from the major PPI network. Fifteen hub genes were discovered, including TNF, PTPRC, ITGB2, ITGAM, VCAM1, CXCR4, TYROBP, TLR4, LCP2, C5AR1, CD86, MMP9, PTPN6, C3, and CXCL10, as well as two transcription factors, i.e. NF-ĸB1 and SPI1. This bioinformatics study has identified the hub genes and key pathways for the development of atherosclerosis in Tibetan minipigs induced by a HFC diet.
Availability of data and materials
The datasets used and analysed during the current study available from the corresponding author on reasonable request.
Complement C5a receptor 1
C reactive protein
C-X-C motif chemokine ligand 10
C-X-C motif chemokine receptor 4
Differentially expressed genes
False discovery rate
Gene Set Enrichment Analysis
Hematoxylin and Eosin
High density lipoprotein cholesterol
Integrin subunit alpha M
Integrin Beta 2
Kyoto Encyclopedia of Genes and Genomes
KEGG Orthology Based Annotation System
Lymphocyte cytosolic protein 2
Low density lipoprotein cholesterol
Molecular complex detection
Matrix metallopeptidase 9
Molecular Signature Database
Normalized enrichment score
Protein tyrosine phosphatase non-receptor type 6
Recombinant Protein Tyrosine Phosphatase Receptor Type C
Smooth muscle cells
Spleen focus forming virus proviral integration oncogene
Search Tool for the Retrieval of Interacting Genes
Toll like receptor 4
Tumor necrosis factor
TYRO protein tyrosine kinase binding protein
Vascular cell adhesion molecule 1
White blood cell
Libby P, Ridker PM, Hansson GK. Progress and challenges in translating the biology of atherosclerosis. Nature. 2011;473:317–25.
Robledo JA, Siccardi LJ. Relationship between genetic and environmental factors and hypercholesterolemia in children. Arch Argent Pediatr. 2016;114:419–25.
Reschen ME, Lin D, Chalisey A, Soilleux EJ, O'Callaghan CA. Genetic and environmental risk factors for atherosclerosis regulate transcription of phosphatase and actin regulating gene PHACTR1. Atherosclerosis. 2016;250:95–105.
Torres N, Guevara-Cruz M, Velazquez-Villegas LA, Tovar AR. Nutrition and atherosclerosis. Arch Med Res. 2015;46:408–26.
Daugherty A, Tall AR, Daemen M, Falk E, Fisher EA, Garcia-Cardena G, Lusis AJ, Owens AP 3rd, Rosenfeld ME, Virmani R. Recommendation on design, execution, and reporting of animal atherosclerosis studies: a scientific statement from the American Heart Association. Circ Res. 2017;121:e53–79.
Seok J, Warren HS, Cuenca AG, Mindrinos MN, Baker HV, Xu W, Richards DR, McDonald-Smith GP, Gao H, Hennessy L, et al. Genomic responses in mouse models poorly mimic human inflammatory diseases. Proc Natl Acad Sci U S A. 2013;110:3507–12.
Priyadharsini RP. Animal models to evaluate anti-atherosclerotic drugs. Fundam Clin Pharmacol. 2015;29:329–40.
Vilahur G, Padro T, Badimon L. Atherosclerosis and thrombosis: insights from large animal models. J Biomed Biotechnol. 2011;2011:907575.
Yongming P, Zhaowei C, Yichao M, Keyan Z, Liang C, Fangming C, Xiaoping X, Quanxin M, Minli C. Involvement of peroxisome proliferator-activated receptors in cardiac and vascular remodeling in a novel minipig model of insulin resistance and atherosclerosis induced by consumption of a high-fat/cholesterol diet. Cardiovasc Diabetol. 2015;14:6.
Pan Y, Rong Y, Huang J, Zhu K, Chen J, Yu C, Chen M. Lower cardiovagal tone and baroreflex sensitivity associated with hepatic insulin resistance and promote cardiovascular disorders in Tibetan minipigs induced by a high fat and high cholesterol diet. J Diabetes Complicat. 2019;33:278–88.
Mozaffarian D. Dietary and policy priorities for cardiovascular disease, diabetes, and obesity: a comprehensive review. Circulation. 2016;133:187–225.
Seo D, Wang T, Dressman H, Herderick EE, Iversen ES, Dong C, Vata K, Milano CA, Rigat F, Pittman J, et al. Gene expression phenotypes of atherosclerosis. Arterioscler Thromb Vasc Biol. 2004;24:1922–7.
Tan L, Wang Z, Li Y. Rabbit models provide insights into bone formation related biological process in atherosclerotic vascular calcification. Biochem Biophys Res Commun. 2018;496:1369–75.
Proto JD, Doran AC, Subramanian M, Wang H, Zhang M, Sozen E, Rymond CC, Kuriakose G, D'Agati V, Winchester R, et al. Hypercholesterolemia induces T cell expansion in humanized immune mice. J Clin Invest. 2018;128:2370–5.
Mailer RKW, Gistera A, Polyzos KA, Ketelhuth DFJ, Hansson GK. Hypercholesterolemia induces differentiation of regulatory T cells in the liver. Circ Res. 2017;120:1740–53.
Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, Vilo J. g: Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019;47:W191–w198.
Merico D, Isserlin R, Stueker O, Emili A, Bader GD. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS One. 2010;5:e13984.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.
Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45:D362–d368.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Scardoni G, Petterlini M, Laudanna C. Analyzing biological network parameters with CentiScaPe. Bioinformatics. 2009;25:2857–9.
Vilela M, Chou IC, Vinga S, Vasconcelos ATR, Voit EO, Almeida JS. Parameter optimization in S-system models. BMC Syst Biol. 2008;2:35.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11.
Janky RS, Verfaillie A, Imrichová H, Van de Sande B, Standaert L, Christiaens V, Hulselmans G, Herten K, Naval Sanchez M, Potier D, et al. iRegulon: from a gene list to a gene regulatory network using large motif and track collections. PLoS Comput Biol. 2014;10:e1003731.
Van Assche T, Hendrickx J, Crauwels HM, Guns PJ, Martinet W, Fransen P, Raes M, Bult H. Transcription profiles of aortic smooth muscle cells from atherosclerosis-prone and -resistant regions in young apolipoprotein E-deficient mice before plaque development. J Vasc Res. 2011;48:31–42.
Liu W, Zhao Y, Wu J. Gene expression profile analysis of the progression of carotid atherosclerotic plaques. Mol Med Rep. 2018;17:5789–95.
Hansson GK, Hermansson A. The immune system in atherosclerosis. Nat Immunol. 2011;12:204–12.
Ammirati E, Moroni F, Magnoni M, Camici PG. The role of T and B cells in human atherosclerosis and atherothrombosis. Clin Exp Immunol. 2015;179:173–87.
Herrero-Fernandez B, Gomez-Bris R, Somovilla-Crespo B, Gonzalez-Granado JM. Immunobiology of atherosclerosis: a complex net of interactions. Int J Mol Sci. 2019;20:5293.
Zhu Y, Xian X, Wang Z, Bi Y, Chen Q, Han X, Tang D, Chen R. Research progress on the relationship between atherosclerosis and inflammation. Biomolecules. 2018;8:80.
Koltsova EK, Hedrick CC, Ley K. Myeloid cells in atherosclerosis: a delicate balance of anti-inflammatory and proinflammatory mechanisms. Curr Opin Lipidol. 2013;24:371–80.
Glass CK, Witztum JL. Atherosclerosis. The road ahead. Cell. 2001;104:503–16.
Soeki T, Sata MJMR. Inflammatory biomarkers and atherosclerosis. Int Heart J. 2016;57:134–9.
Woollard KJ, Lumsden NG, Andrews KL, Aprico A, Harris E, Irvine JC, Jefferis A-m, Fang L, Kanellakis P, Bobik A, Chin-Dusting JPF. Raised soluble P-selectin moderately accelerates atherosclerotic plaque progression. PLoS One. 2014;9:e97422.
Qin H, Liu P, Lin S. Effects of Astragaloside IV on the SDF-1/CXCR4 expression in atherosclerosis of apoE(−/−) mice induced by Hyperlipaemia. Evid Based Complement Alternat Med. 2015;2015:385154.
Heller EA, Liu E, Tager AM, Yuan Q, Lin AY, Ahluwalia N, Jones K, Koehn SL, Lok VM, Aikawa E, et al. Chemokine CXCL10 promotes atherogenesis by modulating the local balance of effector and regulatory T cells. Circulation. 2006;113:2301–12.
Gravastrand CS, Steinkjer B, Halvorsen B, Landsem A, Skjelland M, Jacobsen EA, Woodruff TM, Lambris JD, Mollnes TE, Brekke O-L, et al. Cholesterol crystals induce coagulation activation through complement-dependent expression of monocytic tissue factor. J Immunol. 2019;203:853–63.
Patzelt J, Mueller KA, Breuning S, Karathanos A, Schleicher R, Seizer P, Gawaz M, Langer HF, Geisler T. Expression of anaphylatoxin receptors on platelets in patients with coronary heart disease. Atherosclerosis. 2015;238:289–95.
Clerc G, Rouz PM. Lymphocyte subsets in severe atherosclerosis before revascularization. Ann Intern Med. 1997;126:1004–5.
Zhanglin M, Fen W, Yunfeng S. Identification of key genes and miRNAs associated with carotid atherosclerosis based on mRNA-seq data. Medicine (Baltimore). 2018;97:e9832.
Zirlik A, Maier C, Gerdes N, MacFarlane L, Soosairajah J, Bavendiek U, Ahrens I, Ernst S, Bassler N, Missiou A, et al. CD40 ligand mediates inflammation independently of CD40 by interaction with Mac-1. Circulation. 2007;115:1571–80.
Wang HM, Gao JH, Lu JL. Pravastatin improves atherosclerosis in mice with hyperlipidemia by inhibiting TREM-1/DAP12. Eur Rev Med Pharmacol Sci. 2018;22:4995–5003.
Kamata T, Yamashita M, Kimura M, Murata K, Inami M, Shimizu C, Sugaya K, Wang CR, Taniguchi M, Nakayama T. src homology 2 domain-containing tyrosine phosphatase SHP-1 controls the development of allergic airway inflammation. J Clin Invest. 2003;111:109–19.
Dubois MJ, Bergeron S, Kim HJ, Dombrowski L, Perreault M, Fournes B, Faure R, Olivier M, Beauchemin N, Shulman GI, et al. The SHP-1 protein tyrosine phosphatase negatively modulates glucose homeostasis. Nat Med. 2006;12:549–56.
Nai W, Threapleton D, Lu J, Zhang K, Wu H, Fu Y, Wang Y, Ou Z, Shan L, Ding Y, et al. Identification of novel genes and pathways in carotid atheroma using integrated bioinformatic methods. Sci Rep. 2016;6:18764.
Jeffrey CY, Sommers CL, Jordan MS, Hua G, Samelson LE, Koretzky GA, Hodes RJ. Inactivation of c-Cbl reverses neonatal lethality and T cell developmental arrest of SLP-76-deficient mice. J Exp Med. 2004;200:25–34.
Sindhu S, Al-Roub A, Koshy M, Thomas R, Ahmad R. Palmitate-induced MMP-9 expression in the human monocytic cells is mediated through the TLR4-MyD88 dependent mechanism. Cell Physiol Biochem. 2016;39:889–900.
Jonas W, Jacqueline M-J, Martin H. Functional polymorphism in ALOX15 results in increased allele-specific transcription in macrophages through binding of the transcription factor SPI1. Hum Mutat. 2006;27:78–87.
Financial support of National Science Foundation of China (No. 3157234 and No. 31970514) and Zhejiang Provincial Natural Science Foundation of China (No. LY16C040001) are gratefully acknowledged.
This work was supported by grants from the National Natural Science Foundation of China (No. 3157234 and No. 31970514), Zhejiang Provincial Natural Science Foundation of China (No. LY16C040001).
Animal experiments conform to internationally accepted standards and have been approved by the animal ethics committee of the Zhejiang Chinese Medical University (Approval Number: ZSLL-2016-0031).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Pan, Y., Yu, C., Huang, J. et al. Bioinformatics analysis of vascular RNA-seq data revealed hub genes and pathways in a novel Tibetan minipig atherosclerosis model induced by a high fat/cholesterol diet. Lipids Health Dis 19, 54 (2020). https://doi.org/10.1186/s12944-020-01222-w
- Tibetan minipig