This animal study was approved by the Institutional Animal Care and Use Committee (IACUC) at Iowa State University and was performed according to the Iowa State University Laboratory Animal Resources Guidelines. Male Sprague Dawley rats (n = 12) were obtained at 6 week of age (151–175 g) from Charles River Laboratories (Wilmington, MA). Current studies were limited to males, building on previous studies that provide the foundation for this study. Rats were individually housed in conventional cages in a temperature-controlled environment (22 ± 2°C) using a 12-h light-dark cycle. All rats were acclimated for 1 week on a standard rat chow diet, after which they were randomly assigned to one of two dietary intervention groups. There were no significant differences in baseline body weight between the two dietary groups (p = 0.62). Rats were placed on either a control CAS-or WE-based diet (Research Diets, New Brunswick, NJ; Table 1) matched for 20% protein (w/w). For 72 h, animals underwent a controlled fasted-refeeding protocol to train them to consume food ad libitum within a 4 h time period (Figure 1) for a serum collection. After training, 7-week old animals were fasted overnight for 10 h with water provided ad libitum followed by controlled feeding (4 g) of either the CAS- or WE-based diet. Serum was collected via the tail vein at 0, 2, 4, 6, and 8 h immediately following refeeding. Following the serum time curve collection, rats were given respective dietary treatments, ad libitum for 2 weeks. Food intake and body weight gain were monitored daily. Following the dietary intervention, 9-week old rats underwent a 12 h overnight fast with water provided ad libitum; rats were anesthetized with a ketamine:xylazine cocktail (90:10 mg/kg bw) via a single intraperitoneal injection. Following anesthesia, whole blood was collected via cardiac puncture for serum sample, followed by procurement of vital organs. Euthanasia was confirmed by performing a bilateral thoracotomy. The epididymal VAT, kidney and liver, were procured, weighed, and stored in RNAlater (Thermo Fisher, Waltham, MA). The fresh PFC was excised and rapidly micro dissected on a chilled platform. A flat edge razor was used to remove and discard the olfactory bulbs and make coronal slices through the brain. Anterior forceps of the corpus collosum (AFCC) was used as a visual landmark to identify the section from which to sub-dissect the PFC. Specifically, the darker tissue between the AFCC and the genus corpus callosum in roughly a diamond shape was dissected as the PFC. After sub-dissection the tissue was weighed and stored in RNAlater (Thermo Fisher, Waltham, MA).

Large- and smallRNAs were extracted using a SPLIT Total RNA Extraction Kit (Lexogen, Greenland, NH). RNA was quantified using a Qubit Fluorometer (ThermoFisher Scientific), and integrity was assessed using a Bioanalyzer 2100 (Agilent Technologies). Large RNA libraries were prepared using an automated protocol for the QuantSeq 3′ mRNA-Seq Library Prep Kit (Lexogen, Greenland, NH) and small RNA libraries were prepared using Small RNA Library Prep Kits (Lexogen, Greenland, NH). Total RNA samples were multiplexed across two lanes using an Illumina High-Seq 3000 resulting in an average of 7.5 million reads per sample prior to quality control. Small RNA libraries were also multiplexed and run on a separate lane on an Illumina High-Seq 3000 (San Diego, CA).

All large and small RNA reads were inspected using Fastqc and were trimmed to remove adapter sequences using BBDUK (BBMap: A Fast, Accurate, Splice-Aware Aligner)¹. Read segments matching common Illumina Truseq or Nextera adapter sequences were removed for the reverse-complement or forward sequence of the adapters during processing.

Rattus norvegicus reference genome (fasta) and genome annotation files (gtf) were obtained from the Ensembl genome browser². Reads were aligned to the RNO version 6 release 94 of the Ensembl genome using the STAR v2.5.2 aligner (19). Transcripts aligning to specific genes were counted using the STAR-quantMode geneCounts function to map transcripts to each genome. Files containing microRNA counts and gene counts are located in the GEO database. Small RNA samples were processed using the smallrnaseq python tool (20), which aligns samples using Bowtie and quantifies small RNA read counts using reference fasta and gtf files from RNAcentral.org.

Genes were filtered out if they were not expressed in any samples or had fewer than 10 counts in half of the samples for each gene. Data filtering and alignment settings were adapted from Lexogen's QuantSeq 3′ mRNA-Seq Kit and integrated Data Analysis Pipeline on Bluebee® platform according to the manufacturer's instructions. Following quality control, analysis included 10,101 liver genes, 9,217 adipose genes, 11,607 kidney genes, and 99,36 brain genes.

Read normalization was conducted using a weighted trimmed mean of the log expression ratios [trimmed mean of M values (TMM)] method to account for variable sequencing depth between samples. Differential expression analysis was conducted using the DESeq2 (21) package in the R programming language.

Principal Component Analysis (PCA) was conducted for the initial clustering and characterization of RNAseq data. Hierarchical clustering was used to create a dendrogram classifying samples according to similar transcriptomic profiles using Pearson correlation coefficients. Volcano plots were constructed to visualize genes which surpass a log-fold change of >1.0 (upregulated) or <-1.0 (downregulated) to assess biological significance at adjusted p < 0.05. PCA and volcano plots were all constructed using MatplotLib in python version 3.2.0rc1. Heatmaps and hierarchical clustering were constructed using “heatmapper” package in R software.

Pathway-based analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database that contains annotated biological functions for genes (22). All differentially expressed genes (DEGs) were additionally categorized into biological pathways using KEGG pathway analysis and Gene Ontology (GO) analysis based on cellular localization, function, and processes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 database through the online web application (DAVID Functional Annotation Bioinformatics Microarray Analysis)³. MicroRNA target genes were predicted through miRPath v.3 analysis tools (23). Genetic Interactome Analysis was completed on DEGs from the VAT and PFC using Cytoscape version 3.8.0 (24). STRING application—Protein Query generated a node network that allowed functional enrichment data of all DEGs to be generated. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was queried to identify the functional pathways corresponding to each gene in each species. All DEGs were cross-referenced to String—disease, specifically Diabetes, Parkinson's and Alzheimer's with a confidence interval of 0.95 and a maximum number of proteins of 500 (25).

Total RNA from each tissue was aliquoted and reverse transcribed into cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Catalog # 4368813). The cDNA was diluted to 100 ng/μL and qPCR reactions were performed using 200 ng of total cDNA with primers at 300 nM concentration in 10 μL FastStart Sybr Green Master (Roche) according to the manufacturer's instructions. All qPCR reactions were conducted in a Roche LightCycler 96 Real-Time PCR System. Primers sequences for qPCR (Supplementary Table 1) were obtained from Integrated DNA Technology (IDT) and 18s RNA was used as an internal control for normalization in each tissue.

Data analysis was conducted in Python within an IPython notebook unless otherwise specified. Body weight and food intake group means were analyzed using an unpaired t-test to analyze the difference at p < 0.05 between casein and whole egg treatment groups. When applying DESeq2 for differential gene expression analysis, DESeqDataSetFromMatrix was used to generate p-values and adjusted p-values were calculated using the Benjamini-Hochberg method (26) of false discovery rate (FDR) correction. For all analyses, FDR was controlled at 1% and all adjusted p < 0.05 were considered significant. For functional enrichment analysis, gene identifiers were cross-validated through our differentially expressed genes on the basis of significance at FDR correction of <0.05. All qRT-PCR samples and genes were analyzed in triplicate using the Livak Delta-Delta CT method (27).

Differential gene expression analyses of the tissues identified two DEGs in the liver, 22 in VAT, 52 in the PFC, and none in the kidney. Of the 76 DEGs that surpassed multiple testing corrections (adjusted p < 0.05), one gene was differentially upregulated across both the PFC (5.6-fold increase) and VAT (3.2-fold increase) tissue: indolethylamine N-methyltransferase (INMT). Supplementary Table 2 describes all of the DEGs in each tissue, whereas Supplementary Table 3 contains only those with a significant p-value after correction in all tissues.

To examine the functional pathways for the DEGs, mRNAs were mapped to KEGG/GO pathway terms, which are described in Supplementary Table 4. In the PFC, GO function analysis releveled genes modified by whole egg consumption were significantly enriched in Cholesterol biosynthetic pathways and lipid metabolism. Furthermore, a large genetic network of 16 genes in these two pathways were all down regulated (Supplementary Figure 1A). Of the 12 genes differentially regulated which were involved in Oxidation-reduction process, 9 were downregulated in response to whole egg consumption. In the VAT, Gsta3, Gstz1, and Gstp1, all members of the glutathione transferase activity pathway, were found to be up-regulated by whole egg consumption (Supplementary Figure 1B).

Using STRING: disease query analysis to examine disease related genes that were differentially expressed in response to WE diet, we identified one gene in VAT which is known to be associated with Diabetes. Specifically, Secreted Phosphoprotein 1 (SPP1) also known as Osteopontin was found to be down regulated in the VAT (−3.14 log2 fold change, Supplementary Table 3).

Differential expression analyses of microRNAs identified six upregulated microRNAs and 10 downregulated microRNAs across all four tissues based on non-adjusted p < 0.05 (Supplementary Table 5). No microRNAs survived multiple testing correction with FDR correction at 5% value (Supplementary Table 2). miRNAs largely changed in a tissue specific manner, but two miRNAs changed in multiple tissues including miR-10b-5p which was downregulated in both the VAT and PFC (non-adjusted; p = 0.03 and p = 0.02, respectively) and miR-192-5p which was downregulated in both the liver and PFC (non-adjusted; p = 0.02 and p = 0.05, respectively).

To gain a better understanding of the possible regulatory role of miRNAs, all 16 differentially expressed microRNAs (Supplementary Table 5) were mapped against their putative human genetic targets using miRPath, and cross-referenced against the DEGs in the corresponding tissue (Supplementary Table 3). Based on this mapping, three miRNAs were identified whose predicted target genes which were also differentially expressed in response to diet in the corresponding tissues (Table 2). These included miR-10b-5p which was downregulated in the PFC. Its corresponding predicted target, Arrestin Domain Containing 3 protein (ARRDC3), was upregulated. miR-192-5p was also suppressed in the PFC with nine putative targets which were upregulated. These included: SPC Membrane Recruitment Protein 1 (Amer 1), Fatty Acid Binding Protein 3 (Fabp3), Protocadherin 17 (Pcdh 17), FERM Domain Containing 4B (Frmd4b), Kinesin Family Member 1B (Kif1b), Collagen type V alpha 1 chain (Col5a1), zinc finger protein 36, C3H type-like 1 (Zfp3611), NIPA Like Domina Containing 1 (Nipal 1), and Sorting Nexin 33 (Snx33, Table 2). Additional gene targets of miR-192-5p in the PFC and miR-125b-5p in the VAT were predicted, but did not follow traditional patterns of miRNA regulation as they and the miRNAs were all down regulated (Table 2).

There was no significant effect of dietary treatment on serum microRNA at 0, 2, 4, 6, or 8 h immediately following refeeding based on multiple testing corrective measures using FDR adjusted p < 0.05 (data not shown).

To examine the similarity between transcriptomic profiles, low expression genes across all tissues were filtered and data illustrated in a PCA, as shown in Figure 2A. The PCA plots indicate that these rat samples cluster by diet (WE vs. CAS). Subsequently, the top DEGs were visualized using a hierarchical clustering heatmap for each tissue that displays distinct differential expression patterns according to each diet, as shown in Figures 2B–D. To visualize the magnitude of log-fold changes, all DEGs were plotted for each corresponding tissue using volcano plots. Generally an equal distribution of both up and down regulated genes were observed in each tissue in response to WE diet. Very few genes were differentially expressed at greater than 2 log₂ (Fold Change) in either a positive or negative direction (Figure 3A, Liver; Figure 3B, VAT; Figure 3C, PFC).

To validate the RNAseq data, five DEGs were randomly selected for qRT-PCR analysis across each tissue. The PCR gene expression validation data are contained in Supplementary Figure 2. suggesting that RNAseq and PCR results followed the same relative fold trend for each corresponding gene that was validated. A value of 1 indicates no fold change where as those values above the line indicate an increase in fold change and those below indicate a decrease.

Descriptive data indicated no significant effect of dietary treatment on food intake or total energy intake throughout the study (data not shown). Additionally, there was no statistically significant effect of dietary treatment on final body weight or cumulative body weight gain at 9 week of age after 2 weeks and 3 days on the WE diet. There were also no significant differences in organ weights following tissue collection except for the liver, where rats on the WE-based diets had 16% increased relative liver weight (p < 0.01).

In this study, we examined the effects of consuming a WE-based diet for 2 weeks on gene expression across multiple tissues in healthy rats. Our main results revealed 76 novel DEGs across the PFC, liver, VAT, and kidney transcriptomes. For VAT, GO analyses highlighted that glutathione metabolism was upregulated by WE consumption. Furthermore, SPP1 was identified to be highly down regulated in the VAT after WE consumption and should be further explored as it relates to outcomes of WE consumption and Diabetes. In the PFC, a large network of genes from the cholesterol biosynthetic and lipid metabolic process were all down regulated in the PFC. In addition, multiple genes involved in the oxidation-reduction process in the PFC were also altered by WE consumption. These results indicate that, even short-term, WE-based diets modulates the expression of genes. Together, these results may provide information for follow up in models where these genes and pathways are involved in the pathogenesis of disease.