Diurnal Patterns of Gene Expression in the Dorsal Vagal Complex and the Central Nucleus of the Amygdala – Non-rhythm-generating Brain Regions

Genes that establish the circadian clock have differential expression with respect to solar time in central and peripheral tissues. Here, we find circadian-time-induced differential expression in a large number of genes not associated with circadian rhythms in two brain regions lacking overt circadian function: the dorsal vagal complex (DVC) and the central nucleus of the amygdala (CeA). These regions primarily engage in autonomic, homeostatic, and emotional regulation. However, we find striking diurnal shifts in gene expression in these regions of male Sprague Dawley rats with no obvious patterns that could be attributed to function or region. These findings have implications for the design of gene expression studies as well as for the potential effects of xenobiotics on these regions that regulate autonomic and emotional states.


INTRODUCTION
Diurnal rhythm is a major factor in many aspects of mammalian function, including the regulation of gene expression (Reppert and Weaver, 2002;Sukumaran et al., 2010). In mammals, such diurnal regulation of expression is attributed to the presence of an endogenous circadian clock located in the suprachiasmatic nucleus (SCN) of the anterior hypothalamus that responds to solar time (Dibner et al., 2010). This circadian clock establishes and maintains the daily rhythms of physiological processes (Bugge et al., 2012). Perturbations in the circadian clock can lead to multiple disorders including cancer (Cadenas et al., 2014), major depressive disorder (Li et al., 2013), and alcoholic liver disease (Forsyth et al., 2015).
Peripheral tissues including the heart, kidney, liver and gastrointestinal tract also demonstrate gene expression variations based on circadian time (Oishi et al., 1998;Shieh, 2003;Lim et al., 2006;Hoogerwerf et al., 2007;Sládek et al., 2007;Hughes et al., 2009;Froy, 2011). For example, the pharmacokinetics and dynamics of xenobiotic metabolism in the liver have been shown to vary based on circadian time, thereby altering the balance between therapeutic response and toxicity (Paschos et al., 2010;Froy, 2011;Gachon and Firsov, 2011). Within the brain, most studies of clock-modulated gene expression effects have focused on the SCN. However, there also have been studies that have shown spatiotemporally controlled diurnal patterns of clock-controlled genes in the brain (e.g., Harbour et al., 2014), suggesting that although the SCN serves as master regulator of circadian rhythms, non-SCN brain regions also respond to and intone a diurnal rhythm. By most estimates, the number of diurnallyassociated genes is modest; in mouse, it has been estimated that 5-25% of the transcriptome exhibits diurnal patterning (Panda et al., 2002;Li and Zhang, 2015). Here, we report widespread diurnal gene expression patterns in two non-rhythmgenerating brain regions.
The dorsal vagal complex (DVC) and the central nucleus of the amygdala (CeA) (Figure 1) are involved in receiving visceral inputs and generating autonomic outputs. We and others have found that the DVC, a major viscerosensory region in the brainstem, responds to a variety of inputs in order to maintain homeostasis in response to drugs, such as alcohol (Covarrubias et al., 2006;McDonald et al., 2008). Furthermore, the DVC has projections to the CeA that can modulate emotional state and behavior (Schwaber et al., 1982;Rosa et al., 2014). Moreover, the CeA is implicated in response to many drugs and plays an important role in addiction (Smith and Aston-Jones, 2008;Koob and Volkow, 2010;Roberto et al., 2012;de Guglielmo et al., 2019). We aimed to investigate the consequential influence these viscerosensory functions have on emotional state and behavior in the face of xenobiotic exposure. The focus herein stems from a multidimensional analysis of the resulting data, which revealed considerable patterns of diurnal variation in gene expression in the DVC and CeA. To the best of our knowledge, neither of these brain nuclei has been classified previously as rhythmgenerating or solar-time-responsive like the SCN. However, it is reasonable to hypothesize that state change in these brain nuclei with respect to circadian time could have substantial influence on the subjective emotional experience of xenobiotics -particularly drugs of abuse and withdrawal from them.
Of the 145 transcripts measured in this study, 114 were determined to have a diurnal expression pattern in the DVC and/or the CeA. Of those, only 67 had diurnal patterns in both the DVC and the CeA, and only 22 of that subset were determined to have the same pattern in both brain regions. Taken together, this suggests broad spatiotemporally controlled effects of diurnal time on gene expression in the brain. Furthermore, diurnal gene expression patterns in the DVC and the CeA may underlie diurnal rhythms in physical and emotional regulation, which presents a potential mode of modulation by exogenous stressors that affect these systems. Ultimately, these results suggest that diurnal time should be controlled for explicitly in experimental designs involving gene expression in any brain tissue, and likely any in vivo tissue source.

Animals
All protocols were approved by the Institutional Animal Care and Use Committee of Thomas Jefferson University. Forty young, male, Sprague Dawley rats (35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45)Harlan Laboratories,Indianapolis,IN,United States) were housed individually in the Thomas Jefferson University Alcohol Research Center Animal Core Facility and fed standard chow and water until their weight reached 120 g (approximately 10 days). Upon achievement of the weight criteria, the rats were given an alcohol-free, maltose-dextrin substituted, Lieber DeCarli liquid diet (BioServ, Frenchtown, NJ, United States) ad libitum for 3 days. This diet is a liquid control diet that is nutritionally balanced to provide all necessary nutrients without supplementation (Lieber and Decarli, 1994). On subsequent days, each rat received a restricted quantity of the liquid diet such that the 24 h delayed caloric intake equaled the previous day's intake of matched alcohol-fed animals. Only the calorie-restricted, alcohol-naïve animals are considered in this study. These animals and their alcohol-fed matches were also considered elsewhere (Freeman et al., 2012a(Freeman et al., ,b, 2013. Rat weights were recorded weekly and linearly interpolated for intermediate days. The animal facility was maintained at constant temperature and humidity with a regular 12/12 light cycle [lights on at 07:00 (ZT0)].
The animals were sacrificed via rapid decapitation at one of three times of the day. For clarity, we use the standard circadian nomenclature of Zeitgeber time (ZT) to define the three times of the day, however, we note that this is not a study of circadian rhythms; the light cycles of all animals are identical and undisrupted. The three time points are: ZT3 (n = 7), ZT5 (n = 18), and ZT9 (n = 15). There were no differences among groups in the number of days on the calorie-restricted diet (ZT3: 39.71 ± 1.83 days; ZT5: 39.67 ± 0.95 days; ZT9: 39.07 ± 0.92 days; p > 0.3), and the overall rate of weight gain did not differ between groups (ZT3: 3.77 ± 0.39 g/day; ZT5: 3.73 ± 0.24 g/day; ZT9: 3.70 ± 0.24 g/day; p > 0.7). Thus, although these caloricallyrestricted animals may be different from standard chow-fed control animals, we consider all of the alcohol-naïve animals examined in this study to be experimentally equivalent to each other except for the Zeitgeber time of sacrifice; any baseline differences in expression due to caloric restriction are not examined herein.

DVC and CeA Microdissection
Brains were excised immediately after decapitation and placed into ice-cold artificial cerebral spinal fluid (ACSF: 10 mM HEPES, pH 7.4; 140 mM NaCl; 5 mM KCl; 1 mM MgCl 2 ; 1 mM CaCl 2 ; 24 mM D-glucose). Brainstem and forebrain were separated and secured in individual agarose (4% UltraPure TM low melting point agarose (Invitrogen, Carlsbad, CA, United States) in ACSF) blocks for slicing. Using a McIlwain Tissue Chopper (Gamshall, United Kingdom), the brainstem was sliced into 275 µm transverse sections, and the forebrain was sliced into 625 µm coronal sections. The slices were floated in ice-cold ACSF, and slices containing the regions of interest were selected. The DVC was identified and micropunched as previously reported (Khan et al., 2008). In an analogous manner, the CeA was identified by neuroanatomical landmarks defined by Paxinos and Watson (2007) and micropunched using size-matched micropunches (≤1.25 mm; Stoelting, Wood Dale, IL, United States). Bilateral punches were treated as a single sample. Each sample (two per animal: DVC and CeA) was tested separately and analyzed individually without pooling.

High-Throughput qRT-PCR
RNA was extracted with either the RNeasy or the AllPrep DNA/RNA extraction kit (Qiagen, Valencia, CA, United States). Total RNA was DNAse treated (DNA-Free RNA kit, Zymo Research, Orange, CA, United States), RNA concentration and integrity were assessed with an ND-1000 (NanoDrop, Wilmington, DE, United States) and RNA nano 6000 chips on an Agilent 2100 Bioanalyzer, and the RNA was stored at −80 • C. Reverse transcription was performed from 100 ng total RNA using SuperScript II (Invitrogen, Carlsbad, CA, United States).
A panel of genes representing various cellular processes and cellular localizations were developed for investigation (Supplementary Table S1). The panel is not exhaustive, and the genes were not selected with regard to circadian rhythm regulation nor caloric restriction, but rather were selected to provide a wide perspective of cellular activity (primarily for investigating the role of xenobiotic exposure). Cellular localization and function of the 145 genes are illustrated in Figure 1C. Intron-spanning PCR primers and probes were designed using Roche's Universal Probe Library Assay Design Center 1 . The primers and probes for the 145 genes that passed quality control measures are listed in Supplementary Table S1. The standard BioMark TM protocol was used to pre-amplify cDNA samples for 14 cycles (TaqMan R PreAmp Master Mix: Applied Biosystems, Foster City, CA, United States). qRT-PCR reactions were performed with 96.96 BioMark TM Dynamic Arrays (Fluidigm, South San Francisco, CA, United States), enabling simultaneous quantitative measurement of up to 96 transcripts in 96 samples each (Spurgeon et al., 2008;Orina et al., 2009). Runs were 40 cycles (15 s at 95 • C; 5 s at 70 • C; 60 s at 60 • C). C T values (Supplementary Tables S2, S3) were calculated using default settings in the Real-Time PCR Analysis Software (Fluidigm) and software-designated failed reactions were discarded from analysis.

Experimental Design and Statistical Analysis
Fresh DVC and CeA samples were collected from rats immediately following rapid decapitation [ZT3 (n = 7), ZT5 (n = 18), and ZT9 (n = 15)]. RNA was extracted, reverse transcribed, and measured with microfluidic qPCR (see above). For each sample, the raw C T values were quantile-normalized ( C T ) in R version 2-11-1 (Bolstad et al., 2003). Quantile normalization has been shown to be viable alternative to normalization via an internal control, or housekeeping gene (Fujita et al., 2006;Pradervand et al., 2009).
C T values were calculated for each gene by subtracting the mean C T at ZT5. − C T values were used for all visualizations of the relative expression changes across time points.
All statistical analyses were performed in R using the C T values. Time point and region effects were calculated with a standard 2-factor ANOVA model. Region-specific tests of time point effects were conducted with one-way ANOVA and both post-hoc two-tailed, Student's t-tests and post hoc Tukey comparisons. Multidimensional scaling analysis (MDS) in three coordinate dimensions was conducted in R using cmdscale. All statistical testing was conducted at a 95% confidence level, and all stated errors are the 95% confidence intervals around the means. By nature of the quantity of genes and statistical comparisons (145 genes per statistical comparison type), we anticipate that these analyses will yield some false positives. However, we accept these errors without correction in order to describe the broad impacts of diurnal time. For each gene measured, the expression data, statistical analysis, and pattern determination (section Diurnal Pattern Classification) are provided in Supplementary

Diurnal Pattern Classification
Genes were determined to have a diurnal expression pattern if: (1) the time point effect in the one-way ANOVA was significant 1 www.universalprobelibrary.com (p ≤ 0.05); (2) t-test comparisons revealed a significant difference in gene expression between any two time points (ZT3 vs. ZT5, ZT5 vs. ZT9, or ZT3 vs. ZT9, p ≤ 0.05), or (3) the mean − C T value at any time point was greater than 0.25 or less than -0.25. These statistics are listed in Supplementary Tables S2, S3. Genes fulfilling at least one of these criteria were further classified into eight patterns according to the − C T magnitudes and the statistically significant differences between time points. Support for these patterns were recorded as either statistics-based (for those genes satisfying either or both of the first two criteria above) or threshold-based (for those genes satisfying only the third criterion). The threshold value of 0.25 was selected to capture the genes with larger average changes (∼20% or more) that are not selected by statistical criteria due to a high level of biological variability. There is no reason to expect that the DVC and the CeA have the same diurnal expression patterns, so each brain region's pattern was determined separately.

Alignment Score Calculation
In order to determine genes whose diurnal expression patterns differ between brain regions, we developed an alignment score based upon the calculated angle between two three-dimensional vectors representing the temporal expression patterns. For each gene, m, and each brain region, i, we calculated the mean C T value for each time point and combined the three values as a three-dimensional vector: We then calculated an alignment score (AS) between the DVC and the CeA: where the final transformation from AS * to AS in Equation 3 expands the scale to accentuate small changes in C T that are meaningful, but do not drastically alter the angle between the vectors, θ. With this transformation, a small (AS) indicates more divergent diurnal expression patterns, whereas a large AS indicates diurnal expression pattern similarity. Note that the denominator in Equation 2 ensures that any changes between brain regions that are not time-dependent do not influence the alignment score. Therefore, the alignment score is a quantitative measure of diurnal expression pattern changes in different brain regions.

Multidimensional Scaling Shows Regional and Time Point Clustering
A multidimensional scaling (MDS) analysis (Figure 2) was performed to reduce the high-dimensional dataset into a FIGURE 2 | 3-dimensional multidimensional scaling plots. Three-dimensional multidimensional scaling (MDS) analysis of gene expression from all genes and all samples reveals that brain region differences contribute the largest differentiating factor, but that diurnal time points within a single region cluster together in these complex dimensional coordinates. Each point represents the gene expression profile from a single sample. Distance between points in 3-dimensional space represents overall differences in gene expression between samples. visualization of the overall gene expression. Figure 2A shows that gene expression across the 145 genes assayed differed substantially between the DVC and CeA brain regions based on tight clustering of same-region samples. Samples from the same time point and region did show some clustering (Figure 2B), and the time point differences were accentuated when considering each region separately with subsequent region-only analyses (Figures 2C,D). These results suggest large and numerous gene expression differences between regions but also complex differences between time points within regions.

Gene Expression Changes in the DVC and/or the CeA Across a 6 h Window
Two-factor ANOVA revealed significant time point effects for 43 of the 145 genes tested, indicating statistically significant changes in gene expression across the 6 h window from ZT3 to ZT9. However, the major source of expression variability in most genes was brain region; 107 of the 145 genes had a significant brain region effect with two-factor ANOVA (111 with one-way ANOVA), suggesting that brain region differences in gene expression may obscure other changes when the two regions are studied together. Subsequent intra-region one-way ANOVAs revealed 33 genes with significant time point effects in the DVC and 20 genes with significant time point effects in the CeA. Eight genes (Adrbk1, Agtrap, Creb1, Ptpn2, Mbd1, Nr4a2, Slc35b2, and Terf1) were found to have significant time point effects in both brain regions. T-tests revealed an additional 18 genes in the DVC and 8 genes in the CeA that had significant differences between time points (the full description for all genes can be found in Supplementary Table S4). In summary, 64 genes in the DVC and 54 genes in the CeA, for a total of 75 of the 145 genes tested had statistically significant gene expression changes in the DVC and/or the CeA due to diurnal variations.

Diurnal Pattern Classification Reveals Asynchronous Dynamic Expression Patterns in the DVC and the CeA
In addition to the 64 genes found in the DVC with statistics-based patterns, 26 genes were determined to have only threshold-based diurnal expression patterns (3 genes demonstrated statisticsbased patterns only). We normalized expression to the ZT5 time point for each gene resulting in eight possible non-zero gene expression patterns to characterize the gene expression dynamics across the three time points in this study (3 2 -1 = 8). There are at least 5 genes with each classified pattern in the DVC (Figure 3, left), indicating that diurnal expression patterns are asynchronous and complex. Similarly, in the CeA (Figure 3, right), 38 genes demonstrated threshold-based expression patterns only (6 genes had statisticsbased patterns only). The increase in threshold-based patterns in the CeA may indicate a greater degree of biological variability in CeA samples or a broader range of interactions between emotional state and basal gene expression rhythms.

Functionality and Cellular Localization of Transcripts With Diurnal Patterns
Diurnal patterns of gene expression do not appear to be associated with cellular localization (Figure 4). Rather, transcripts in all parts of the cell showed large changes in expression across the 6 h period. Interestingly, some functionally related groups of transcripts have similar patterns, suggesting coordinated regulation of expression. For example, transcripts of regulators of G-protein signaling (Rgs proteins), such as those located in the top right corner in each panel of Figure 4, have similar diurnal regulation in the DVC, particularly at ZT9. Glutamate receptors (iGluR's) have a variety of diurnal patterns that are relatively similar in the two brain regions, whereas the expression patterns of gammaaminobutyric acid (GABA) receptors (GABAR) seem to be inverted between brain regions.

Pattern Similarity and Differences Between the DVC and the CeA
Of the 145 transcripts measured, 114 were determined to have a diurnal expression pattern in the DVC and/or in the CeA. Of those, 67 had diurnal patterns in both the DVC and the CeA, and only 22 of that subset were determined to have the same pattern in both brain regions ( Figure 5A). In order to quantify pattern similarity across the two brain regions, an alignment score was calculated for each gene based on the average expression at each time point (see section Materials and Methods). Genes with higher alignment scores have similar diurnal gene expression patterns in the DVC and the CeA, whereas genes with lower alignment scores have divergent patterns in the DVC and CeA. As shown in Figure 5B, 12 genes were identified to have a high alignment score (AS > 5.5, which represents approximately the top 10% of genes measured) and 12 different genes were identified to have a low alignment score (AS < 3.7, representing the lowest 10% of genes measured). An additional 23 genes with a moderately low alignment score (3.7 < AS < 4.1) and 22 genes with a moderately high alignment score (5.0 < AS < 5.5) were also identified. Three families of genes were observed to be of particular interest in this list of 67 genes: the regulators of G-protein signaling (Rgs) family, glutamate receptors, and transcription factors. As shown in Figure 5C, in general, Rgs family member transcripts had different diurnal expression patterns in the DVC and the CeA, indicating differential roles in these two brain regions in diurnal regulation, whereas glutamate receptors had very similar diurnal expression patterns in the two brain regions. For transcription factors, the relationships among gene expression patterns were complex: some had high alignment scores and some had low alignment scores.

DISCUSSION
Temporal gene expression patterns are well-studied in circadian research. Microarray experiments have led to the identification of genes with rhythmic circadian expression patterns in the SCN (Panda et al., 2002;Reppert and Weaver, 2002) and peripheral tissues (Panda et al., 2002;Storch et al., 2002). Other groups have studied the circadian rhythms of clock-related genes such as Per1, Per2, and Bmal1 in the CeA (Lamont et al., 2005;Amir and Robinson, 2006;Perrin et al., 2006;Segall et al., 2006;Ramanathan et al., 2008;Harbour et al., 2014) and in the DVC (Mieda et al., 2006;Kaneko et al., 2009) in rodents. However, we are not aware of any medium-or high-throughput studies of diurnal gene expression patterns in either the CeA or the DVC. The integral viscerosensory role of these nuclei implies that such diurnal gene expression patterns could be critical to first-person emotional experience in response to, for example, withdrawal from addictive substances. Thus, this phenomenon may play an important role in drugseeking.
Discernable circadian rhythms in gene expression were found in less than 10% of all genes in the SCN, liver, and heart (Panda et al., 2002;Storch et al., 2002;Vollmers et al., 2009). Therefore, although we expected to find some genes with diurnal expression patterns, we were surprised to find that more than half of the genes tested in each brain region had at least one statistically significant difference between two individual time points. We note that these high percentages depend on the subset of genes assayed, and may not be representative of global expression in the CeA and DVC. We speculate this could result from overrepresentation of inter-and intra-cellular signaling-involved genes. Furthermore, the patterns of gene expression we describe do not include phasic patterns that are characteristic of circadian expression patterns -three time points across a 9 h period provides a resolution too low to resolve phasic patterning. These experiments were not designed to capture circadian expression FIGURE 3 | Diurnal patterns in gene expression in the DVC and CeA. Heatmaps of − C T values in DVC samples (left) and CeA samples (right) organized according to diurnal pattern (center). In the heatmaps, each column shows data from a single sample. For each brain region, samples are separated according to diurnal time (ZT3, ZT5, and ZT9). Each row shows the expression level of a single gene relative to the mean ZT5 expression in that brain region. Data for the 114 genes with a diurnal pattern in either the DVC or the CeA are shown. Two genes had patterns in the DVC, but expression levels were not detectable in the CeA. Therefore, only data from the remaining 112 genes are shown for the CeA.
patterns and we cannot deduce such patterning from these data. However, the groups of animals at different time points were not different in any respect, except the time of day at which they were sacrificed. Genes assayed were not chosen as circadian rhythm regulators or as those known to play a role in generating diurnal rhythms, yet the expression of these genes showed distinct time-of-day patterns at the three time points investigated. Because a high percentage of the genes assayed demonstrated this disparate expression behavior in brain regions critical to physical and emotional homeostatic regulation, we conjecture that similar expression patterns likely exist throughout the brain, and this kind of "hidden variability" could have considerable impact on the design and execution of gene expression experiments. Additionally, it is prudent to keep in mind that the scope of this study was limited to all male animals on a calorie-restricted diet. Differential expression with respect to diurnal rhythms may be profoundly influenced by hormonal signaling or other physiologic processes that differ between the sexes. Further, calorie restriction is known to influence the circadian clock itself (Froy, 2007). While these limitations may reduce generalizability of the results, they should not lessen the impact of the general suggestions of these findings.
The results from our present study suggest that the DVC and CeA, which are critical in autonomic and emotional regulation, implement time-dependent gene expression patterns that underlie their baseline regulatory activity. We conjecture that these pattern changes would affect corresponding metabolic or protein enrichment changes, but in this study, we are unable to assess downstream effects. Nonetheless, these results serve FIGURE 4 | Cellular and signaling pathway view of diurnal gene expression changes. The center panel shows the genes measured in a layout matching that of Figure 1C in order to depict the gene expression according to the cellular and signaling pathways of interest. Average − C T values for each gene are used to colorize the subpanels using the color scale shown. The three columns show the three diurnal time points (ZT3, ZT5, and ZT9). The top row shows data from the DVC while the bottom row shows data from the CeA. By definition, the average − C T value of each gene at ZT5 is zero, and therefore the two brain regions share a common center panel.
FIGURE 5 | Pattern similarity and differences in the DVC and CeA. (A) Each cell shows the number of genes with the pattern combination found by matching the pattern in the DVC, rows, with the pattern in the CeA, columns (DVC patterns shown at left, CeA patterns shown at top). The 18 genes with the same pattern in the DVC and the CeA are included in the diagonal (yellow). Dotted box outlines highlight several cells with genes that have more diurnal dynamics in the CeA than in the DVC. Double box outlines highlight the reverse: cells with genes that have more diurnal dynamics in the DVC than in the CeA. The cells shaded in gray indicate transcripts with no diurnal pattern in the alternate brain region. (B) Histogram of the alignment scores of all 145 genes. Genes with low alignment scores (<3.7) like that illustrated at the top left and genes with high alignment scores (>5.5) like that illustrated in the top right are shaded in dark gray. The genes in these categories are listed in order of increasing alignment score within the corresponding region. Genes with moderately low (<4.1) and moderately high (>5.0) alignment scores are shaded in light gray. (C) The average − C T values of Rgs family members (column 1), glutamate receptors (column 2), and transcription factors (column 3) with low, moderately low, moderately high, or high alignment scores are shown. All genes within the group with high or moderately high alignment scores are shown in the top row, while all genes within the group with low or moderately low alignment scores are shown in the bottom row. Each gene is shown in a different color with solid lines representing the gene expression data from DVC samples and dotted lines representing the gene expression data from CeA samples. as cautionary advice for subsequent gene expression studies in the brain: diurnal time should be considered explicitly in experimental design regardless of brain region for meaningful comparisons. Furthermore, these data lay a foundation for investigating xenobiotic-induced homeostatic or emotional disturbances due to diurnal changes in gene expression patterns in the CeA and DVC beyond any global changes induced by these perturbations. Diurnal gene expression variations in these nuclei may influence the physical dependence or addictive properties of certain drugs. Conversely, xenobiotics may influence the diurnal rhythms in these nuclei such that dysregulation occurs leading to negative emotional states or behavior change.

DATA AVAILABILITY STATEMENT
The datasets generated in this study are publicly available on the Gene Expression Omnibus (https://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc=GSE147569), under the accession no. GSE147569.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of Thomas Jefferson University.

AUTHOR CONTRIBUTIONS
MS designed the study and conducted experiments, performed analysis, generated figures, and was a main contributor to the writing of the manuscript. SO'S performed analysis, generated figures, and was a main contributor to the writing of the manuscript. RV, GG, and JS designed the study and were involved in analysis, figure design, and editing. KK and BO were involved in analysis and editing. All authors discussed the results and commented on the manuscript.

FUNDING
This work was supported by NIH T32AA007463, R01AA015601, and U01HL133360.