Transcriptome-wide RNA sequencing analysis of rat skeletal muscle feed arteries. I. Impact of obesity

Animals and experimental design.

LETO (n = 12) and OLETF (n = 12) rats were obtained at age 4 wk (Japan SLC, Hamamatsu, Shizuoka, Japan). The OLETF rat, characterized by a mutated cholecystokinin-1 receptor that results in a hyperphagic phenotype, is an established model of obesity, insulin resistance, and type 2 diabetes (30). Animals were individually housed in a temperature-controlled (21°C) environment with light from 0600 to 1800 and dark from 1800 to 0600. All animals were given ad libitum access to standard chow with a macronutrient composition of 56% carbohydrate, 17% fat, and 27% protein (Formulab 5008, Purina Mills, St. Louis, MO). Rats were anesthetized at 30–32 wk of age with an intraperitoneal injection of pentobarbital sodium (50 mg/kg) between 0800 and 0930. Tissues were then harvested, and the animals were killed by exsanguination. Food was removed from the cages 12 h prior to death. All protocols were approved by the University of Missouri Animal Care and Use Committee.

Body weight, body composition, food intake, and citrate synthase.

Body weights and food intakes were monitored and recorded on a weekly basis. Weekly food intakes were averaged across ages 20–30 wk. Body composition was assessed by dual-energy X-ray absorptiometry (DXA; Hologic QDR-1000, calibrated for rodents) on the day of death. Omental and retroperitoneal adipose tissue depots were then removed and weighed to the nearest 0.01 g. Citrate synthase activity was measured from whole muscle homogenate of the vastus lateralis with the spectrophotometric method of Srere (70).

Isolation of skeletal muscle feed arteries and aortic endothelial cells.

Immediately after the gastrocnemius-plantaris-soleus muscle complex was harvested, the muscles were pinned down in a petri dish containing a cold RNA-stabilizing agent (RNAlater; Ambion, Austin, TX). Under the microscope, the SFA and GFA were then dissected clean of perivascular adipose tissue and excess adventitia as described previously (3, 28, 34, 39, 42, 82, 83). The single GFA supplying the medial head of the gastrocnemius muscle was used for the present study. In our experience, OLETF rats typically have one to three SFA. All SFA were dissected and pooled for the present RNA-Seq analysis. The reader is referred to our recent publication for a visual of the anatomic location and structure of GFA and SFA (28). Aortic endothelial cells were isolated by gentle scraping of longitudinally opened aortas as described previously (6, 50, 52, 66). This method of scraping the luminal surface yields an endothelial cell-enriched sample (50). The aorta was chosen to provide sufficient yield of endothelial cells. Indeed, we have established that we recover ∼200 ng of total RNA per rat aortic endothelial scrape, which is above the minimum required for RNA-Seq experiments. Samples were kept in RNAlater for 48 h at 4°C and then removed from the RNAlater solution and stored at −80°C until analysis.

Blood parameters.

Whole blood was collected on the day of euthanasia for analysis of glycosylated hemoglobin (HbA1c) by the boronate-affinity high-performance liquid chromatography method (Primus Diagnostics, Kansas City, MO) in the Diabetes Diagnostics Laboratory at the University of Missouri. Serum samples were prepared by centrifugation and stored at −80°C until analysis. Glucose, triglyceride (TG), and nonesterified fatty acid (NEFA) assays were performed by a commercial laboratory (Comparative Clinical Pathology Services, Columbia, MO) on an Olympus AU680 automated chemistry analyzer (Beckman-Coulter, Brea, CA) using commercially available assays according to manufacturers' guidelines. Plasma insulin concentrations were determined with a commercially available rat-specific enzyme-linked immunosorbent assay (Alpco Diagnostics, Salem, NH). Samples were run in duplicate, and manufacturers' controls and calibrators were used according to assay instructions.

RNA extraction.

Total RNA isolations were performed with the NucleoMag 96 RNA kit procedure (Clontech part no. 744350.1), which is a single-tube method based on reactive magnetic bead technology designed for automated small- or large-scale preparation of highly pure total RNA from tissue or cell samples. All liquid handling was optimized for use on a Beckman3000 robotic liquid handler housed within a laminar flow hood (with UV decontamination) designed to ensure a clean room environment for working with microtissues, which yield low (pg to ng) amounts of RNA. Briefly, groups of 24 sample vessels were removed from −80°C and immediately homogenized for 60–120 s in their own 2-ml microtube with the Mini-Beadbeater-24 (BioSpec Products) in the presence of NucleoMag lysis buffer and several miniature chrome-steel (RNase treated) BBs. Care was taken to get complete microvessel disruption without extending grinding times to prevent the generation of excess heat. The resulting homogenate was then loaded onto the robot deck, and a digital photo was taken before the sample was transferred into 96-deep well microplates. The photo allowed us to have a physical record of each sample ID prior to loading into the microplate for accurate tracking purposes. This process was repeated four times to completely fill a 96-well plate within 10–20 min. The combination of using stabilized tissue and immediate homogenization in chaotropic salt-containing lysis buffer ensured that the RNA was protected from RNase degradation during tissue disruption. After homogenization, the RNA was bound to RNA beads in the presence of alcohol, and a magnet was used to perform several wash and elution steps in a completely automated fashion. This method included a DNase digestion step ultimately yielding RNA of similar yield and quality from column-based procedures. Immediately after RNA isolation, pure RNA was transferred to a new 96-well plate and a 5-μl aliquot was taken into a second plate for RNA quality control. Both plates were stored at −80°C with cryogenic plate seals and placed in secondary containment to prevent frost buildup on the plates during storage.

RNA quality control (concentration and integrity).

For assessing total RNA yield and integrity, tandem Agilent Bioanalyzer 2100 instruments were used in combination with the Total RNA 6000 Pico Assay. At the time of this study, the RNA Pico LabChip Kit was the only platform to give unbiased assessment of RNA integrity (via RIN) and accurate results with extremely low RNA concentrations such as those provided by microvessels. The qualitative range for the total RNA assay is 200–5,000 pg/μl, and the most important advantage of this system is the small amount of sample used (1 μl), leaving the rest of the RNA for other applications. Typical yields from rat microvessels were ∼500–1,000 pg/μl. RNA quality control was performed with only the aliquot from each isolation plate.

Illumina library preparation (SMARTer amplification and RNA-Seq).

Because of the low yields of total RNA from microvessels, total RNA could not be used directly in traditional Illumina gene expression profiling methods (RNA-Seq) because of the low concentration of the samples (standard RNA-Seq kits during this project required 0.1–1 μg of total RNA). Thus the SMARTer Ultra Low RNA Kit for Illumina Sequencing (Clontech, catalog no. 634935) was utilized for generating full-length cDNA transcripts prior to Illumina RNA-Seq library preparation. Briefly, the technology involves SMARTer first-strand cDNA synthesis and purification, utilizing the SMARTer anchor sequence and poly(A) sequence that serve as universal priming sites for end-to-end generation of single-stranded cDNA, followed by cDNA amplification with LongDistance PCR (LD-PCR) using the manufacturer's recommended Advantage 2 PCR system (Clontech, catalog no. PT3281-1) containing a novel polymerase and ultrapure dNTPs specifically for Illumina sequencing. Using the concentration values from the Bioanalyzer RNA Pico Assay, we sought to use 100–1,000 pg of total RNA as input to the SMARTer 1st cDNA reaction.

After cDNA generation, validation was performed with the Bioanalyzer 2100 High Sensitivity DNA Assay (Agilent, catalog no. 5067-4626) for select samples from each plate of 96 samples in order to accurately size and quantitate DNA up to 12 kb in length, again consuming minimal sample volumes (1 μl). After 14 cycles of LD-PCR amplification the final cDNA yields were estimated at ∼1–10 ng for each microvessel, which is a suitable input amount for library preparation for cDNA/ChIP Seq library preparation. To generate Illumina Paired-End RNAseq libraries, cDNA was fragmented to ∼200 bp with the Q700 DNA fragmentation system (QSonica) and then used directly with the NextFlex DNA preparation kit (Bioo Scientific, catalog no. 5140-02) with some modifications. Briefly, fragment cDNA was end-repaired and purified with 1.8x SPRI beads to remove reaction components (Agencourt). The resulting blunt ends were A-tailed in preparation for cohesive ligation to the Illumina specific sequencing adapters diluted to 0.6 μM working concentration (NextFlex DNA Adapters, Bioo Scientific, catalog no. 514104). Ligated DNA was purified twice with 1.0x SPRI to remove adapter dimers and perform gel-free size selection and then amplified through 14 cycles of PCR. The final sequencing construct was purified with a 1.0x SPRI to remove low-molecular-weight adapter dimer artifacts (if any), and libraries were validated to contain ∼330-bp fragments with the Bioanalyzer 2100 High Sensitivity DNA Assay. Library quantitation was performed with the Qubit 2.0 fluorometer and the High Sensitivity DNA assay (Life Tech, catalog no. Q32851).

RNA sequencing.

By utilizing 48 unique adapter indexes during library preparation across each plate of 96 libraries created, we were able to overcome several common technical mistakes. First, it allowed us to account for technical biases through randomization of samples by vessel type (group) and treatment across each plate of samples. Second, by having a priori knowledge of the sequencing index used to identify each sample from the sequencing, we were able to use a single manufacturing lot of adapters that were uniformly diluted and preplated to ensure similar ligation efficiencies across several plates (hundreds of samples) used in the study. Third, this indexing scheme allowed us to standardize the pooling of several libraries by row within each plate, where equimolar volumes of each sample in a plate row were pooled to a final concentration of 5–10 nM. Altogether, this approach prevented inadvertent use of the wrong adapters during preparation, randomized the sample and index combinations, and allowed for reduced mixing up of libraries within each sequencing pool. The final pools (>50 total) were each loaded on a single lane of single-read 50-base sequencing on the Illumina HiSeq2000 and ultimately yielded ∼175–200 million useable reads per lane (14–17 million reads per RNAseq sample). It should be noted that the combination of the 48 adapters that resulted in 4 pools of 12 indexes was carefully designed and wet lab tested to be compatible with the HiSeq to ensure maximum sequence yields and to ensure that each sample was correctly identified by the HiSeq during the index identification steps.

Statistical analysis.

The primary analysis of the preprocessed RNA-Seq data was conducted with R/Bioconductor, with other free, open-source programs used for preprocessing. The preprocessing steps consisted of 1) quality assessment of raw sequenced library data (FASTQ files) with FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc), 2) trimming of adapters from each library with Cutadapt (40), 3) aligning libraries to the rn4 genome build (http://hgdownload-test.cse.ucsc.edu/goldenPath/rn4/bigZips/) with a seed-and-vote method (37) and saving in SAM format, and 4) forming a read-count matrix at the gene level by read summarization using featureCounts in Rsubread.

Nonspecific filtering of genes prior to statistical testing was carried out to increase detection power (5), based on the requirement that a gene have an expression level greater than 16 counts per million reads mapped (CPM) for at least 9 libraries across the GFA and SFA samples. For the nonspecific filtering of the aortic endothelial cells, the reader is referred to our companion article (49a). In either case, this CPM cutoff was established empirically based on the point at which the ERCC Spike-Ins at different concentrations were no longer distinguishable for that particular set of data. Library normalization to adjust for differences in library size was then made with the TMM method (59) in edgeR (58). This was used in conjunction with a flexible mean-variance modeling and transformation process known as Voom (35), as part of the limma package (69), a combination known as limmaVoom that was recently shown to perform among the best of available techniques in an comprehensive comparison of differential gene expression methods for RNA-Seq data (56). With this framework, normalized expression levels were compared between the groups of interest and significant differences identified by fold change and adjusted P values. Adjustment to the P values was made to account for multiple testing with the false discovery rate (FDR) method of Benjamini and Hochberg (4). For analysis of obesity effects (i.e., OLETF vs. LETO within GFA and within SFA) as well as analysis of vessel effects (i.e., GFA vs. SFA within LETO and within OLETF), we chose 10% as our FDR threshold for statistical significance. For analysis of interactive effects of obesity and vessel, computed as (OLETF GFA/LETO GFA)/(OLETF SFA/LETO SFA), we adopted a FDR of 20%, knowing that there is less power to detect an interactive effect than a main effect (e.g., obesity or vessel). For all comparisons, a blocking factor on the flow cell was used in the fitted model to adjust for systematic flow cell effects, and the robust version of limma was utilized during model fitting. For between-vessel comparisons, we accounted for the dependency due to multiple measurements coming from the same rat using the duplicateCorrelation option (69) in limma, but the correlation was negligible (∼6%) within the block [rat]. As an empirical measure of the FDR, we evaluated what proportion of the identical ERCC probe/concentration combinations (Set B) appeared in our list of differentially expressed genes. Similarly, we looked at a set of 13 putative housekeeping genes derived from a study of more than 13,000 rat samples (13) to have another estimate of our FDR. The set of genes was Actb, B2m, Gapdh, Gusb, Hprt1, Hmbs, Hsp90b1, Rpl13a, Rps29, Rplp0, Ppia, Tbp, and Tuba1. For both of these sets of controls, we also estimated the fold change of each of the genes as a measure of the accuracy of the fold change estimates.

Networks were generated through the use of Ingenuity Pathway Analysis (Ingenuity Systems, www.ingenuity.com), henceforth IPA, as previously described (52). The full lists of all differentially expressed genes were uploaded into IPA with Entrez GeneID, fold change, and adjusted P value. In the case of duplicate mapped IDs, the median values (fold change and adjusted P value) were used to represent the single results for that ID. Genes that were not detected, or those that were filtered out (as previously described), were assigned a fold change of 1 and an adjusted P value of 1. An adjusted P value <10% was set to identify molecules whose expression was significantly differentially regulated, which IPA terms as Network Eligible Molecules (NEMs). Our networks were built using only knowledge from experimentally observed relationships contained in the Ingenuity Knowledge Base (IKB) and were algorithmically generated to maximize their specific connectivity with each other. Higher network scores imply a lower probability of finding the observed number of NEMs in a given network by random chance. Specifically, the score is the negative log₁₀ of the P value from Fisher's exact test applied to a given network. For example, a score of 9 for a network implies a 1-in-a-billion chance of obtaining a network containing at least the same number of NEMs by chance when randomly picking 35 molecules from the IKB. A detailed description of the network-generating algorithm is provided by Calvano et al. (7).

Graphical representations of the networks were generated with IPA's Path Designer, which illustrates the relationships between molecules. Molecules are represented as nodes, and the biological relationship between two nodes is represented as an edge (connecting line). All edges are supported by at least one reference from the literature, from a textbook, or from canonical information stored in the IKB. These are rich, high-information content graphics, with details included in the figure legends.

Finally, between-group differences for all descriptive variables were determined by using an independent t-test, for which statistical significance was declared at P ≤ 0.05.

Impact of obesity on GFA and SFA.

Figure 1 displays the number of differentially expressed genes in GFA and SFA between OLETF and LETO rats. Obesity produced 282 upregulated and 133 downregulated genes in SFA and 163 upregulated and 77 downregulated genes in GFA (FDR < 10%), with an overlap of 93 genes between the arteries. Of the 93 genes altered in both SFA and GFA (overlapping portion of Venn diagram in Fig. 1), 68 were upregulated and 23 were downregulated in both vessels. Although the direction of change for most genes was the same in both feed arteries, two genes were oppositely responsive to obesity between the arteries, ceruloplasmin (Cp) and aspartyl aminopeptidase (Dnpep) (Fig. 1). One (Cp) was upregulated in SFA and downregulated in GFA, whereas the other (Dnpep) was downregulated in SFA and upregulated in GFA. Table 2 (GFA) and Table 3 (SFA) provide lists of the top 20 genes differentially expressed between LETO and OLETF groups (sorted by magnitude of fold change), and the full lists of differentially expressed genes are provided in Supplemental Tables S1 and S2.¹ Figure 2 and Figure 3 illustrate the top-scoring gene networks influenced by obesity in the GFA and the SFA, respectively. The scores for these gene networks were 76 in the GFA and 66 in the SFA.

Between-vessel differences in LETO and OLETF rats.

Figure 4 presents the number of genes differentially expressed between GFA and SFA in LETO and OLETF rats. In LETO rats, there were 89 upregulated and 114 downregulated genes in the GFA compared with the SFA. There were 244 upregulated and 275 downregulated within OLETF (FDR < 10%) in the GFA compared with the SFA. Of the 76 differentially expressed genes common to both LETO and OLETF, 25 were upregulated and 51 were downregulated in GFA compared with SFA. The directional differences were consistent for each of these genes (Fig. 4). Table 4 (LETO) and Table 5 (OLETF) provide lists of the top 20 genes differentially expressed between GFA and SFA, and the full lists of differentially expressed genes are provided in Supplemental Tables S3 and S4. Figure 5 and Figure 6 illustrate the top-scoring gene networks for between-vessel differences in LETO and OLETF rats, respectively. The scores for these gene networks were 45 in LETO (Fig. 5) and 43 in OLETF (Fig. 6) rats.

Obesity × vessel interaction analysis.

Table 6 presents the full list of genes differentially influenced by obesity between the vessels, i.e., the result of our interaction analysis. There were 16 genes identified as upregulated, i.e., for which the effect of obesity was greater in the GFA than in the SFA, and 20 downregulated genes for which the effect of obesity was greater in the SFA than in the GFA (FDR < 20%). The top-scoring gene network for the interaction analysis is presented in Fig. 7, and the score for this network was 36.

Impact of obesity on aortic endothelial cells.

A total of 396 transcripts were found to be differentially expressed between LETO and OLETF in aortic endothelial cell-enriched samples (FDR < 10%), of which 235 were upregulated in OLETF and 161 were downregulated in OLETF compared with LETO. Table 7 provides a list of the top 20 genes differentially expressed between the groups, and the full list is provided in Supplemental Table S5. The top-scoring gene network for the effects of obesity on aortic endothelial cells is illustrated in Fig. 8, and the score for this network was 38.

Between-vessel differences in LETO and OLETF.

We assume that the list of top genes differentially expressed between GFA and SFA within LETO (Table 4) reflects innate differences between these feed arteries in the healthy state. Several genes identified seem to reflect developmental differences between the vessels. For example, Mab21l2, Fhod3, Apod, and Srpk3 are all involved in developmental processes (29, 46, 60, 62). Periaxin is a myelin-specific protein, suggesting developmental differences between the arteries in the degree of innervation, consistent with previous evidence of heterogeneous innervation among skeletal muscle vascular beds (38). Similarly, pathway analysis revealed a network of genes related to developmental and structural differences between the vessels in LETO rats (Fig. 5). The network includes molecules related to cellular assembly (actins, alpha catenin), vascular development (Notch), and myosin family genes (Mylpf, Mylk, Myh14, collectively represented by the Mlc node in Fig. 5). These differentially expressed genes may reflect the greater smooth muscle content in the GFA compared with the SFA, which also may be of developmental origin. In addition to developmental differences, it is important to note that the GFA and SFA differ markedly with respect to local hemodynamic signals (e.g., shear stress) owing to differences in soleus and gastrocnemius muscle recruitment patterns during standing and walking (i.e., normal cage activity) in rats (32). As such, these local hemodynamic factors, which we and others have shown are important determinants of vascular function and phenotype (27, 28, 33, 49, 51, 53, 66, 74), could have contributed to differential transcriptomic profiles between GFA and SFA in LETO rats.

For the between-artery comparison in OLETF rats, we observed some of the same intrinsic differences as discussed above within LETO (e.g., Mab21l2 and Pck1). Thus it appears that some developmental/intrinsic differences between these arteries are unaltered by obesity. Moreover, a major finding of our study was that the number of genes differentially expressed between the vessels within OLETF was ∼2.5-fold greater than that within LETO, indicating greater heterogeneity between these feed arteries in the obese OLETF rat (Fig. 4). The known functions of selected genes listed in Table 5 suggest that the GFA takes on a more inflamed, proatherogenic phenotype compared with the SFA within the OLETF rat, including Nnat, a novel marker of endothelial inflammation in diabetes (45), Cidec, which is involved in ectopic fat accumulation (47), Gpihbp1, an atherosclerosis-related gene involved in the transport of lipoprotein lipase into capillaries (11, 79), and Ccl11, also called eotaxin, an emerging inflammatory marker of cardiovascular diseases in humans (17). In addition to these novel genes, we found some “classic” indications of a proatherogenic phenotype in the GFA compared with the SFA of OLETF rats, such as an upregulation of Vcam1, an established marker of inflamed vascular tissue, and a downregulation of Nos3, the endothelial isoform of nitric oxide synthase (Supplemental Table S4). Similarly, pathway analysis of genes differentially expressed between GFA and SFA of OLETF rats revealed a network involving Nos3 and Akt, atheroprotective genes that were downregulated in GFA compared with SFA, as well as upregulations in a number of mitochondrial NADH transporters linked to mitochondrial complex 1 in the GFA (Fig. 6). Considering the large body of work relating reduced expression and activity of Nos3 and increases in mitochondrially derived superoxide production to vascular pathologies (64), it seems reasonable to interpret this top-scoring network as evidence of a preferential shift of the GFA toward a dysfunctional, prooxidant state compared with the SFA in the obese OLETF rat.

Effects of obesity on GFA and SFA.

Both the GFA and the SFA displayed transcriptional alterations in response to hyperphagia-induced obesity in the OLETF rat. Examination of the list of obesity-altered transcripts in the GFA (Table 2 and Supplemental Table S1) revealed a number of genes previously implicated in the development of a dysfunctional/proatherogenic vascular cell phenotype. For example, several genes (Folr2, Hemgn, Pxmp4, Pla2g2a, Klra17, Ly49Q, and Oas1) have been reported to be involved in proinflammatory processes (20, 26, 36, 55, 76, 78), suggesting a more inflamed phenotype of GFA in OLETF compared with LETO. Notably, the gene Nnat, elevated by more than fivefold with obesity, was recently identified as a novel marker of endothelial inflammation in diabetic mice (45). The gene Srebf2, another novel proatherogenic marker, contains an intronic region that encodes miR33, which suppresses cholesterol transport, ultimately leading to increased risk for atherosclerosis (57). Regarding obesity-induced downregulated genes in the GFA (Table 2, bottom), several appear to be potent mediators of antiatherogenic processes. For example, Abra is highly upregulated in growing collateral vessels, is responsive to nitric oxide, and is increased in response to atheroprotective shear stress (75). Hspa1a (also known as HSP70) and Hspb1 have well-established antiatherogenic effects, characterized as inhibitory against inflammatory processes such as leukocyte adhesion (16, 54). Finally, the obesity-induced downregulation of complement component 2 (C2) in the GFA could be interpreted as a proatherogenic response to obesity, as mice deficient in C2 exhibit significant atherosclerosis (67). Together, this profile of top-responding genes provides evidence that the GFA is adversely affected by obesity at the transcriptomic level in the OLETF rat model.

Several of the top obesity-induced upregulated genes in the SFA were also upregulated in the GFA (Mis18a, Srebf2, Pxmp4, Fam151b, Folr2). Overall, the list of top differentially expressed genes presented in Table 3 suggests the presence of a proatherogenic effect on the SFA as well as the GFA. For example, Folr2, Vnn1, and Deaf1 have all been implicated in atherosclerosis and inflammation (10, 26, 31). Interestingly, the SFA and the GFA shared some of the largest-magnitude obesity-induced downregulated genes, including Cyp2e1 and C2.

Pathway analysis (Figs. 2 and 3) revealed in both GFA and SFA that obesity significantly downregulated expression of gene networks centered on Ubc. Ubc is a ubiquitin gene that encodes a polyubiquitin precursor. The ubiquitin-proteasome system is an intracellular proteolysis system and plays a multifaceted role in vascular cells (14, 71). In endothelial cells a reduction in ubiquitin proteasome system activity reduces the activity of endothelial nitric oxide synthase, and in smooth muscle cells ubiquitin proteasome inhibition increases sensitivity to endoplasmic reticulum stress-mediated cell death and suppresses the unfolded protein response (14, 71). On the other hand, inhibition of the ubiquitin-proteasome system has also been shown to reduce neointima thickness and reduce vascular inflammation in atherosclerosis models (14). Thus further investigation into the role of the ubiquitin system and the functional significance of reduced expression of this system in SFA and GFA of obese rats is warranted.

Obesity × vessel interaction analysis.

As the GFA and SFA have been shown in the OLETF rat model to be differentially susceptible to obesity-induced endothelial dysfunction (3) and display differences in insulin-induced vasomotor reactivity (28), we were particularly interested to know the genes for which there were differential effects of obesity between the GFA and the SFA in the present study. Our approach for examining this was to identify those transcripts that displayed significantly different fold differences between the groups in one vessel compared with the other (Table 6). Interestingly, genes involved in endothelial cell permeability [Pla2g2a (76)] and inflammation [Mgst1 (68), Svil (44), Retstat (61), Cd14 (72)] were found to be significantly upregulated by obesity in the GFA but not the SFA. In contrast, notable genes upregulated in the SFA but not the GFA included novel markers such as Cp, which is associated with incident mortality in heart failure patients (9) and nephropathy onset in type 2 diabetes patients (25), as well as Scd and Dnajc10, both of which have been implicated in endoplasmic reticulum stress (8, 41). Finally, the selective upregulation of three key neuronal genes (Plp1, Ncs1, and Prx) in the SFA suggests that neural innervation or perhaps neural control may be differentially altered between GFA and SFA with obesity, a hypothesis that warrants future research.

Pathway analysis of the gene list presented in Table 6 revealed a top-scoring network centered on CD14, with links to proinflammatory signaling molecules such as Jnk, interleukins 1 and 12, interferon alpha, P38MAPK, and nuclear factor-κB (Fig. 7). Interestingly, CD14 is widely recognized as a cell surface antigen on circulating monocytes but has also been detected in vascular tissues (72). As part of the receptor complex (along with toll-like receptor 4) for lipopolysaccharide, CD14 is upstream of a host of intracellular proinflammatory events. Our present findings lead us to propose that CD14 and the related downstream inflammatory signaling pathways identified in Fig. 7 may play roles in the differential effects of obesity between the GFA and the SFA in the OLETF rat model.

Effects of obesity on aortic endothelial cells.

Recognizing that the GFA and SFA analyses were performed on whole vessel homogenates, we also sought to determine the influence of obesity on a purified population of endothelial cells obtained from aortic scrapes. Overall, the genes identified as differentially expressed between LETO and OLETF aortic endothelial cells were consistent with the hypothesized proatherogenic effects of obesity. Some genes identified as being upregulated in GFA and/or SFA also appear in Table 7, e.g., Hemgn, Pxmp4, Pla2g2a, and Nnat. Additionally, a number of obesity-induced transcripts in aortic endothelial cells have been previously implicated in proatherogenic processes, including Scd1, previously reported to be elevated in mesenteric arteries in a mouse model of atherosclerosis (65); Clec1b, a ligand for activation of inflammatory cells (specifically natural killer cells) (81); Mmp12 and leptin, which have well-established roles in atherosclerosis and whose expression has been shown to be upregulated together in human abdominal aortic aneurism specimens (73); Irf7, previously found to be a central gene in an IPA network of particulate matter-induced inflammation in endothelial cells (18); Npy, a potent vasoconstrictor involved in pathological vascular remodeling (1); and Cct6a, a chaperone protein known to be induced by endoplasmic reticulum stress (15).

The top-scoring gene network in aortic endothelial cells was suggestive of an obesity-induced shift toward a proinflammatory phenotype, as reflected by alterations in a number of molecules that converge on the nuclear factor-κB complex (e.g., obesity-induced increases in Stat2, interferon signaling molecules, and Cxcl16 and a decrease in Sod3). This network therefore suggests that obesity confers a proatherogenic phenotype in endothelial cells of large arteries (i.e., the aorta) in addition to the deleterious effects on gene expression observed in whole vessel homogenates of skeletal muscle feed arteries.

Limitations.

A few limitations of our study warrant mentioning. Regarding our experimental model, the hyperphagic OLETF rat not only is obese but also displays a number of obesity-related comorbidities (insulin resistance, hypertension, fatty liver disease, etc.) that progressively increase in severity as the animals advance in age. We also did not determine whether OLETF and LETO rats differed in terms of daily energy expenditure or cage activity. Thus we cannot rule out that the differential gene expression profiles reported here might be partially attributable to comorbidities or physical activity. Additionally, our aortic endothelial cell gene expression results should not be interpreted as being representative of endothelial cells from all large and small arteries. The aorta was chosen because it is the largest artery and therefore the most convenient sample for the enrichment of arterial endothelial cells.

Conclusions.

In summary, our study indicates that both the GFA and the SFA are susceptible to unfavorable effects of obesity at the transcriptional level in the OLETF rat model. The SFA displayed a greater number of genes differentially expressed between LETO and OLETF rats, which, in light of previous evidence indicating that the SFA largely maintains its functional properties in obese animal models, suggests that the SFA may be more “plastic” in the face of a proatherogenic obese systemic environment. Despite differences in the number of genes found to be differentially expressed, examination of the known functions of genes differentially expressed between lean and obese rats indicated a more pronounced shift toward a proinflammatory/prooxidant phenotype in the GFA compared with the SFA in obese rats, a finding consistent with our previously published finding that the GFA displays greater susceptibility to obesity-induced endothelial dysfunction (3). Our study is also consistent with the large body of research indicating that obesity confers a proatherogenic phenotype in arterial endothelial cells; importantly, our data are the first RNA-Seq-based characterization of these effects. Finally, it is important to emphasize that the majority of differentially expressed genes reported here have little or no known function in vascular tissues. Therefore, our findings should aid in generating new hypotheses that will lead to novel mechanistic insights into the negative consequences of obesity on skeletal muscle resistance arteries.