Short-term high fat diet alters genes associated with metabolic and vascular dysfunction during adolescence in rats: a pilot study

Animal experimental protocol

The male Sprague-Dawley rats (140–160 g, 5–6 weeks of age) used in this study were part of a larger investigation that compared symptoms of metabolic dysfunction in rats on CHOW and HFD treatments (Sweazea, Lekic & Walker, 2010). The composition of each diet has been described previously (Crawford et al., 2019). Briefly, the CHOW diet is the regular maintenance diet of the animals and contained in calories: 18.9% protein, 57.33% carbohydrates (% sucrose NA), and 5% fat (Cat. No. 2018; Harlan Teklad). Half of the animals were switched to a HFD after arrival that modeled a Western diet high in animal-based ingredients and saturated fat. The macronutrient profile of the HFD was 20% protein, 20% carbohydrates (6.8% sucrose) and 60% fat (in %kcal), primarily from lard (Cat. No. D12492; Research Diets, Inc., New Brunswick, NJ, USA) as described by Crawford et al. (2019). This HFD was intentionally selected to be greater than 50% fat to induce and study more severe metabolic syndrome (Rodríguez-Correa et al., 2020). In contrast, the fat content of the chow diet was plant-derived (60% soybean oil and meal, 40% wheat and corn). Similarly, the protein composition of the HFD included a mixture of casein, lactic, and 30 mesh whereas the chow diet was comprised of 36% wheat, 34% corn and 26% soy proteins. The carbohydrate composition of the chow diet included 54% wheat products, 40% corn products, and 6% soybean meal whereas the high fat diet included 63% corn and 37% sucrose.

Rats were maintained on the respective diets for 6 weeks and the food was replaced every 2–3 days to prevent spoiling. All animals had access to food and water ad libitum and were pair-housed in identical cages in the same animal facility and exposed to a 12:12 h light dark cycle. Animal care procedures were reviewed and approved by the Institutional Animal Care and Use Committee of the University of New Mexico School of Medicine (AAALAC # 000222). The thoracic aorta (not including the aortic arch) of eight animals, four each from the normal CHOW and HFD, were removed, cleared of surrounding fat and connective tissue, and the lumen was cleared of blood. The aortic tissue was placed in RNA later (Qiagen, Valencia, CA, USA) and frozen (−80 °C). Total RNA was extracted with the Qiagen RNeasy kit. RNA from two rats was pooled for each library (CHOW1, CHOW2, HFD1, HFD2), so four animals are represented for each treatment and were contained within two libraries. The high reliability of Illumina sequencing eliminates the need for technical replicates (Marioni et al., 2008) so each library was sequenced once.

Library preparation and sequencing

Illumina cDNA libraries were constructed and sequenced at the National Center for Genome Resources (NCGR, Santa Fe, NW) using manufacturer’s standard operating protocols. One μg of total RNA for each sample was subjected to Illumina mRNA-SeqPrep Kit RS-100-0801 following manufacturer’s protocols. This involved two rounds of polyA selection using polyT beads, random fragmentation of RNA, random priming to make cDNA, ligation of bar-coded sequencing primers, and size selection. The bar-coded libraries were run in a single lane of a flow cell on an Illumina Genome Analyzer IIx, resulting in single-end 36 nucleotide (nt) read lengths. Raw data files (.fastq) and processed data files are available through the Gene Expression Omnibus (GEO), accession number GSE147335.

The post-sequencing bioinformatic workflow is diagramed in Fig. 1. Four raw data files, two each for the CHOW and HFD condition were uploaded into Lumenogix. Data quality parameters included FASTQC (Andrews, 2010, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and BLAST (Altschul et al., 1990) analysis of randomly selected reads. Reads were mapped to ENSEMBL release 79 of the rat genome (Rnor_5.0; GCA_000001895.3) using two methods, the Genomic Mapping and Alignment Program (GMAP) (Wu & Watanabe, 2005) and Tophat, which takes into account potential RNA splice sites, known as splice-site junction mapping (SJM) (Trapnell, Pachter & Salzberg, 2009). The results for each mapping program was quantitated per gene using HTSeq count (Anders, Pyl & Huber, 2015). The resulting files were archived in GEO (accession number GSE147335). Empirical Analysis of Digital Gene Expression data in R (EdgeR) (Robinson, McCarthy & Smyth, 2010) was used to evaluate statistical significance for the GMAP and SJM mapping methods, these results are also archived under GEO GSE147335. The resulting data files were downloaded and opened with JMP version 14 (SAS Institute, Cary, NC, USA) for data filtering and further analyses. Zipf’s power law is a graphic method employed to model complex phenomena (Adamic, 2011) that has been used to identify low count theshold in environmental metagenomic shotgun sequence data (Reiss, Guerra & Makhnin, 2016) and in a study of cell response to nanostructured titanium (Reiss et al., 2020). Transcription profiles of tissue qualifies as a complex system and a low count cutoff was determined by a Ziph graph, in which the log abundance per gene, as measured counts per million (CPM) for all samples, is graphed against the log rank of each gene (Fig. 2). The use of Ziph’s law to filter out low expressed genes replaces the application of an arbitrary threshold for gene expression. The requirement that genes be significant at p(adj) ≤ 0.05 in at least one mapping algorithm ensures that only those genes that are highly significant are passed onto the enrichment step.

Data-mining and visualization with cytoscape

Genes with p(adj) ≤ 0.05 in either GMAP or Tophat (SJM) results were passed to Cytoscape (v3.7.1) App ClueGO/CluePedia (v2.5.5) for enrichment and network analyses. Each gene is assigned GO biological process terms (if available) and this information is used to find terms that are enriched in the data. ClueGo parameters used included a two sided-hypergeometric test with an initial cut-off of p(adj) 0.025, corrected with the Bonferroni step down method. Terms were limited to GO levels between 5 and 9 and the kappa score threshold was 0.4. At least three genes per term were required. The GO term fusion GO groups options were selected; these account for the redundancy of GO terms. ClueGO selects GO groups overview terms with the highest significance level (lowest p(adj) value). The log₂ fold-change values were transferred into Cytoscape and displayed as gene node color. Adobe Illustrator (v24.0.1) was used for final mark-up of the networks.

Muscle cell remodeling

Our analysis identified 17 (~19%) genes related to muscle cell remodeling with a number of these genes associated with processes involving the regulation of smooth muscle proliferation and mesenchyme migration. While expression of the majority of these genes were consistent with the negative alterations in vascular dysfunction, these findings were notable given the brevity of the dietary intervention, age of the animals, and time-course for muscle cell remodeling in vascular diseases. Reduced expression of Actc1 parallels findings that occur during the onset of heart disease through induction of cardiomyocyte apoptosis (Jiang et al., 2010) and increased expression of Uchl1, part of the ubiquitin-proteasome system associated with muscle atrophy in diabetic rats (Reddy et al., 2018). Further, Nr4a3 was downregulated, which is a transcription factor connected to all five of the main biological functions included in our analysis. Nr4a3 is a transcriptional target of p53, triggering apoptosis and playing a suppressive role in cancer (Fedorova et al., 2019). While the aorta was analyzed in the present research, previous studies in rats suggest the decreased expression of Nr4a3 may be associated with reduced signaling for metabolic capacity in muscle (Oita et al., 2009; Pearen et al., 2012).

Vascular functions

Genes associated with vascular function were abundant and included terms in the GO biological process network related to the regulation of blood pressure, vasoconstriction, temperature homeostasis, regulation of apoptotic cell clearance, and oxygen carrier activity (Figs. 3 and 4). In addition, many of the genes associated with monovalent inorganic anion homeostasis were downregulated, including Slc4a1. This encoded protein mediates an electroneutral chloride/bicarbonate exchange in physiological conditions and is the main membrane protein in vertebrate erythrocytes (Barneaud-Rocca, Etchebest & Guizouarn, 2013). In addition, Slc4a1 plays a role in signaling aged cells for clearance through macrophage-mediated phagocytosis (Badior & Casey, 2018). Relatedly, genes involved in oxygen carrier activity displayed reduced expression, including Hba1, Hba-a2, Hbb-b1, and Loc689604, a novel hemoglobin alpha gene. Recent studies in human and mice tissue have found hemoglobin alpha subunit expression in blood vessel walls, where it appears to localize at the myoendothelial junction and plays a role in regulating nitric oxide signaling between endothelium and smooth muscle (Butcher et al., 2014). Decreases in associated gene expression in the current study is suggestive of vascular dysfunction, perhaps through reduced nitric oxide signaling. This finding is corroborated with previous physiological data from these HFD animals, that exhibited impaired endothelium-dependent vasodilation (Sweazea, Lekic & Walker, 2010; Sweazea & Walker, 2011).

Egr1 functions as a master transcriptional switch activated by ischemia to trigger expression of pivotal regulators of inflammation, coagulation, and vascular hyperpermeability (Silverman & Collins, 1999; Yan et al., 2000) and was upregulated. Although ischemic tissue damage was not directly assessed in the present study, Erg1 activation is a notable finding considering its central role in the pathogenesis of vascular dysfunction (Yan et al., 2000). Therefore, this finding may be suggestive of early tissue damage occurring in the aortic wall of the HFD animals. Another intriguing finding was the decreased expression of Nampt (aka Visfatin). Nampt is the rate limiting component in the mammalian nicotinamide adenine dinucleotide (NAD) biosynthesis pathway and connects NAD-dependent sirtuin signaling, constituting a strong intrinsic defense system against various stresses (Wang & Miao, 2015). In the present study, we show a significant decrease in Nampt expression in aortic tissue potentially implicating this gene as nutrient-sensitive to the HFD treatment. Together the expression profile of vascular function related genes highlight potential negative alterations occurring in aortic tissue during the early stages of metabolic dysfunction, which are consistent with the impaired endothelium-mediated vasodilation observed in prior studies of ex vivo small resistance arteries from these animals (Sweazea & Walker, 2011; Sweazea & Walker, 2012).

Insulin/peptide hormones

The adipokine encoding gene, Lep, displayed increased expression paralleling the elevated leptin levels in an obesogenic state (Pan, Guo & Su, 2014). In rats, leptin has been reported to stimulate the proliferation and migration of vascular smooth muscle cells in the aorta, implying leptin may play a role in the formation and development of vascular remodeling (Oda, Taniguchi & Yokoyama, 2001). Leptin appears to directly stimulate cellular hypertrophy via p38 MAP kinase in rat vascular smooth muscle cells (Shin et al., 2005). In fact, leptin receptors are present in endothelial cells providing a target for leptin action on the vascular wall (Fortuño et al., 2002).

Hspa1a and Hspa1b encode for heat shock proteins that are major orchestrators of the cellular stress response associated with several human diseases (Daugaard, Rohde & Jäättelä, 2007; McCallister et al., 2015; Richter, Haslbeck & Buchner, 2010). Both are down regulated in the HFD samples. Hspa1a overexpression protects many cell types against a wide variety of stresses, including hypoxia and reoxygenation injury (Hutter et al., 1994; Radford et al., 1996). Indeed, elevated Hspa1a has been reported to provide cardioprotection against ischemia-reperfusion induced myocardial infarction and apoptosis (Quindry et al., 2007). Whereas low levels of intracellular and circulating Hspa1a promote a proinflammatory state and increase the vulnerability of the arterial wall to the damaging action of vascular risk factors involved in endothelial dysfunction—the first stage in the development of the atherosclerotic plaque (Dulin, García-Barreno & Guisasola, 2012). The role of Hspa1b in the vasculature is less apparent and, although closely related to Hspa1a, it is distinct and has remained conserved in a lineage preceding placental mammals (Hess et al., 2018). Therefore, selection may be acting to preserve not only the primary function but also other less well-studied secondary functions of the molecules, like lipid binding and cell signaling (Hess et al., 2018). Gja1 (connexin) was downregulated and similar to Hspa1a and Hspa1b is associated with both vascular functions and insulin/peptide hormones. Gja1 encodes the protein connexin 43, a gap junction protein important in cell-to-cell communication and in vascular myoendothelial gap junctions in the heart and other tissues (Triggle et al., 2012; Xu et al., 2017). In diabetic vascular disease where blood glucose levels are both elevated and not well controlled, advanced glycation end products (AGEs) may also impact both the expression and function of connexins, and thus, contribute to the development of endothelial dysfunction and vascular disease (Wang et al., 2011).

Lipid processes

Similar to the pattern noted with insulin/peptide hormone process, there was an expected decrease in genes encoding for enzymes involved in fatty acid biosynthesis. Of these, Fasn and Acaca encode enzymes in the fatty acid synthesis pathway and are well-known to increase with rising insulin levels associated with elevated carbohydrate intake (Saltiel & Kahn, 2001; Wong & Sul, 2010). Beyond the high level of saturated fat and low carbohydrate consumed by the HFD rats in the present study, the downregulation of these genes may also be related to increased oxidative stress, which has been noted to down regulate intermediates in the insulin signaling pathway and promote insulin resistance (Bloch-Damti & Bashan, 2005; Fridlyand & Philipson, 2006; Hurrle & Hsu, 2017). In relation, the HFD rats showed significant increases in circulating thiobarbituric acid reactive substances (TBARS) as well as vascular oxidative stress as measured by DCF fluorescence (Sweazea & Walker, 2011). In addition, they developed significantly impaired vasodilation in the mesenteric arteries that was mediated by oxidative stress (Sweazea & Walker, 2012).

Another important finding was the decreased expression of genes associated with brown adipose tissue (BAT). Fndc5 gene was reduced in HFD rats. The encoded protein from Fndc5 is a precursor to irisin which reportedly acts on white adipose cells to stimulate uncoupling protein expression and brown-fat-like development (Boström et al., 2012). In relation, we report decreased expression Ucp1 and of LOC100909612, annotated as encoding a potential Upc1-Like protein. Also decreased was Acot11, which catalyze the hydrolysis of fatty acyl-CoA thioesters and is enriched in BAT (Cohen, 2013). Acot11 also appears to regulate the availability of substrates for beta-oxidation and uncoupling (Han & Cohen, 2012; Okada et al., 2016; Zhang et al., 2012). While there is a reduction in BAT in humans with increased obesity and metabolic dysfunction (Franssens et al., 2017), the decreased expression of these BAT genes in the aortic wall from the present study was an unexpected finding. This may have been due to remnants of perivascular adipose tissue (PVAT), the adipose tissue surrounding blood vessels, in our tissue samples. PVAT shares several of the defining characteristics of BAT, including its cellular morphology and expression of thermogenic genes characteristic of brown adipocytes (Hildebrand, Stümer & Pfeifer, 2018). Moreover, PVAT is the source of several paracrine factors that regulate vascular tone (Chang et al., 2012; Friederich-Persson et al., 2017).

Carbohydrate processes

There was a general decrease in genes related to the regulation of carbohydrate metabolic processes including Hk2, which phosphorylates glucose to produce glucose-6-phosphate and is regulated in part by insulin (Osawa et al., 1995; Okada et al., 2016). Development of diabetes is associated with a shift in cardiac glucose metabolism including significant decreases in Hk2 content (Calmettes et al., 2015). The general downregulation in the insulin signaling pathway is notable as HFD rats display several features of metabolic dysfunction including hyperglycemia, impaired glucose tolerance, oxidative stress as well as oxidative-stress mediated vascular dysfunction (Sweazea, Lekic & Walker, 2010; Sweazea & Walker, 2011; Sweazea & Walker, 2012). Because the vasculature is located at the interface between blood and tissue, vessel walls are particularly exposed to this oxidative stress and may experience negative alterations in insulin signaling before other tissues (Wiernsperger, 2003).

Other considerations

In the present study, the number of genes associated with circadian function is suggestive of biological crosstalk between circadian rhythm and metabolism in the vasculature. Specifically, there was a notable decrease in expression of genes related to circadian function, including Nr1da1 (a.k.a. Rev-ERBA), Per1, Noct, and Nampt. In addition, the increase in Lep mentioned previously is intriguing given that leptin has been reported to affect circadian rhythm in mammals (Prosser & Bergeron, 2003). The relationship between the circadian clock and these genes along with their related metabolic pathways may represent a novel target for studying the pathogenesis of vascular diseases. Disruption of daily rhythms with a chronic HFD in mice has been found to be rapidly reversible in liver tissue (via PERIOD2:LUCIFERASE bioluminescence rhythm) (Branecky, Niswender & Pendergast (2015)). Though, it is unclear if this would be similar in the vasculature. In relation, time-restricted feeding in alignment with light and dark cycles has been found to partially restore day-night differences in cardiac remodeling, triglyceride synthesis, lipidome, and metabolic flexibility in mice (Mia et al., 2021). Therefore, future work should investigate HFD cessation, as well as time-restricted feeding models on circadian-related genes in the vasculature. Though this would be complex as the number of animals required, the appropriate time series of sample collection, and selection of vasculature tissue location throughout the cardiovascular system is not fully apparent. Indeed, circadian clock components display differences in expression and localization throughout the cardiovascular system, which may confer nuances of circadian clock signaling in a cell-specific manner (Anea et al., 2018).

In rodents, high fat feeding has been reported to alter circadian synchronization to light (Mendoza, Pévet & Challet, 2008), and obesity reduces circadian gene expression in concert with impaired endothelial function in the aorta (Nernpermpisooth et al., 2015). Two core clock genes, Per1 and Nrd1d, were downregulated in the aorta of rats on the HFD; both are transcriptional repressors. Per1 is reportedly decreased in vascular smooth muscle in type 2 diabetic db/db mice (Su et al., 2012). Rev-ERBA can reversibly bind to heme, which changes its interaction with other clock components (Yin et al., 2007). In the present study, the enzyme responsible for the rate-limiting step in Heme production, Alas1, was downregulated. The role of heme availability in circadian signaling and metabolic dysfunction are unresolved. NAMPT catalyzes the rate limiting step in the NAD salvage pathway and is also a clock-controlled gene (Bass & Lazar, 2016). Nampt-dependent oscillatory production of NAD regulates oscillation of clock target gene expression by releasing the core clock component: Clock-Arntl/Bmal1 heterodimer from NAD-dependent Sirt1-mediated suppression (Ramsey et al., 2009). Noct is another mitochondrial enzyme also under the control of the circadian clock, it catalyzes the dephosphorylation of NADP+ and NADPH. In contrast, Egr1 was increased. Egr1 appears to play a potential role in the regulation of the hepatic circadian clock in response to feeding/fasting (Tao et al., 2015). The role of Th expression in peripheral clocks remains to be elucidated. This provides evidence for the disruption of circadian pathways in aorta are induced by HFD, but the nature of changes remains to be determined. Furthermore, it is not clear if females would respond in a similar fashion. In fact, the response may very well be different in females as the effect of HFDs on the circadian rhythm display some notable dependence on biological sex. For instance, in female mice estradiol has been reported to regulate the daily rhythms underlying diet-induced obesity (Omotola et al., 2019). Ovarian hormones in female mice provide a protective circadian mechanism, preventing disruption of daily rhythms from high-fat feeding (Palmisano, Stafford & Pendergast, 2017). Though the age of mice is an important determinant of this response since a HFD has been noted to alter circadian genes in mammary glands of female mice during puberty (Sundaram, Johnson & Yan, 2020). These mice were the same age as the adolescent male rats used in the present study (6-weeks), therefore, future work should consider female animals in the context of puberty, as well as menopause. Such research may better delineate the putative role of estrogen as a protective agent against HFD-induced circadian disruption.

The current study has some limitations. Due to the nature of tissue extraction and the limited number of samples, it was not possible to assess time-related changes in gene expression leading to vascular dysfunction. Therefore, significant changes may have occurred at earlier or later time points. Another important aspect relates to the HFD diet used in the present rodent model. These rats were not genetically prone to develop atherosclerotic lesions, therefore the diet is likely unable to induce significant vascular disease. However, our previous work showed impaired function in these animals (Sweazea, Lekic & Walker, 2010; Sweazea & Walker, 2011). A third limitation of the present study was the small sample size in each group, although this is typical of sequencing studies of this nature as evidenced by similar sample sizes reported in other research (Bao et al., 2016; Xi et al., 2017). As a pilot study, this work is intended to inform future research efforts in this area, with larger samples. We acknowledge that the generated statistics and described effects could be within normal biological variation and need to be replicated with larger scaled research. Such research has been performed in other animal models, including a 12-month HFD study in adult pigs (Milenkovic et al., 2020), however this study and others investigating the vasculature have used microarray gene assessments (Bao et al., 2016). The technological platform employed is not trivial as microarray analysis may induce a more biased view of gene expression and does not allow the whole transcriptome to be assayed in comparison to RNA-seq (Rao et al., 2019). Importantly, the study was able to demonstrate that RNA-seq has the sensitivity to detect differences in gene expression in the given the experimental protocol. Another point to emphasize is that this data validates several findings established in the physiological literature (e.g., involvement of circadian clocks) and experimentally (e.g., the increased circulating leptin levels, hypertension and impaired vasodilation previously observed in these animals). It is important to note that gene expression data is a proxy for active protein concentrations. As such, there are four principles of data mining using GO terms that are important to understand its potentials and weaknesses. First, differential gene expression data is a direct measure of events in the nucleus. Although an mRNA template is required to produce proteins, it is subjected to post-transcriptional control and subsequent proteins require post-translational modifications, so mRNA-based assays are a proxy for fully active protein. Second, transcripts from the same gene may code for multiple proteins, some of which act in a tissue-specific manner. An example of this in Figs. 3A and 4 are the inclusion of the GO terms ‘memory’ and ‘learning;’ these terms reflect the role in the nervous system. GO term enrichment alone does not take into account the tissue and experimental protocol. Third, by their nature, gene ontologies are complex, and this is reflected by the redundancy in networks of GO terms. For example, vasoconstriction is an overview term that 4.85% of the 86 genes share, and both the negative and positive regulation of vasoconstriction are included as GO BP terms since these processes involve the same genes. The inclusion of multiple levels of GO terms (in this case 5–9) is necessary to capture an optimized subset of genes. Thus, any conclusions regarding specific gene products will need to be empirically validated.