Transcriptome Analyses of Heart and Liver Reveal Novel Pathways for Regulating Songbird Migration

Many birds undertake long biannual voyages during the night. During these times of the year birds drastically reduce their amount of sleep, yet curiously perform as well on tests of physical and cognitive performance than during non-migrating times of the year. This inherent physiological protection disappears when birds are forced to stay awake at other times of the year; thus these protective changes are only associated with the nocturnal migratory state. The goal of the current study was to identify the physiological mechanisms that confer protection against the consequences of sleep loss while simultaneously allowing for the increased physical performance required for migration. We performed RNA-seq analyses of heart and liver collected from birds at different times of day under different migratory states and analyzed these data using differential expression, pathway analysis and WGCNA. We identified changes in gene expression networks implicating multiple systems and pathways. These pathways regulate many aspects of metabolism, immune function, wound repair, and protection of multiple organ systems. Consequently, the circannual program controlling the appearance of the migratory phenotype involves the complex regulation of diverse gene networks associated with the physical demands of migration.

analysis (WGCNA) and pathway enrichment analyses with Ingenuity Pathway Analysis software (IPA) to analyze changes in individual genes and gene regulatory network pathways altered while maintaining the migratory state. Our results demonstrate that the phenotype associated with migratory behavior is extremely complex, involving changes in multiple genetic regulatory events and thus likely not regulated by a single switch. We present multiple regulatory pathways that may provide protective effects on the heart and liver and that can be further mined to understand mechanisms that are engaged for optimal performance during resiliency towards sleep deprivation.

Methods
Animals. 72 White-throated sparrows (Zonotrichia albicollis) were captured during the fall of 2014 using mist nets in Centre County, PA (40.808749, -77.858566) and transported less than 5 km to housing at the Pennsylvania State University animal facility. All birds were housed together in a large flight cage (approximately 2.2 meters on each side) for at least six months to standardize recent life history and avoid confounds such as variations in diet and movement. Birds were given ad-libitum food (Mazuri small bird feed, Mazuri, Richmond, IN) and water. Initially, the light cycle was set to match the environmental photoperiod at the time of capture, but once collection was complete, the light cycle was adjusted by 30 minutes/week in a single step per week until 12.5L:11.5D was achieved. Following acclimation, animals were individually housed in small wire cages (approximately 20 × 28 × 38 cm) on racks within the larger enclosure so that migratory status could be evaluated. While flight was moderately restricted in these smaller cages, birds remained in visual and auditory contact with the remaining birds in the enclosure. In order to accelerate the transition between migratory states for tissue collection, the photoperiod was manipulated to either short (6L:18D) or long (18L:6D) days, which has previously been shown to be sufficient to shift the migratory state of birds in the laboratory. Once animals transitioned to the desired migratory state, they were returned to 12.5L:11.5D for at least one week prior to dissections, while migratory state was continually monitored (see below). This photoperiod stabilizes both the migratory and non-migratory states, and performing dissections under these conditions allowed comparisons between migratory and non-migratory conditions without confounds such as total light exposure or dissections between migratory statuses under different light conditions. All animals successfully transitioned between migratory and non-migratory states within approximately 1-2 month of altered light cycle, and were stable in the migratory condition for another 1-2 months, roughly approximating seasonal timing in the wild. There are two morphs of white-throated sparrows; to reduce variability from sex and morph differences we restricted our analyses to males of the white morph. Sex was verified by gonad inspection during dissections. All procedures followed the guidelines outlined in the NIH Care and Use of Laboratory Animals guide, and were reviewed and approved by the Pennsylvania State University IACUC Committee; Protocol #4558-1. Animals were collected under Commonwealth of Pennsylvania Wildlife Collecting Permit #32318 and USFWS Permit # MB170276-1.
Determining migratory status. Animals were continually monitored by video and behavior was scored for a minimum of 4 nights immediately prior to tissue harvest to ascertain migratory status. Zugunruhe includes a constellation of active behaviors that occur during the night that we observed with video scoring, including increased nighttime vigilance, grooming and perch hopping. However, we have found the most robust marker that separates migratory versus non-migratory birds is a wing-whirring behavior 6,8 . Animals are unable to fly in the laboratory housing conditions; however, they are able to orient to a desired direction for migration and perform a rapid flutter of wings, which has been interpreted as a proxy for flight. This behavior is seen frequently during the night in birds in the migratory state, and never in birds that are non-migratory 6,8 . Therefore, recordings were scored in 10 minute bins for presence of or absence of wing-whirring. If the animal exhibited Zugunruhe (i.e. wing whirring) in a 10 minute bin it was given a score of one, and if it did not it was given a score of zero. In addition, fat deposit scores have previously been used to determine migratory state and were recorded during tissue collection as verification. tissue collection. All tissue collection was performed under the same housing conditions (ad libitum food and water availability, 12.5 L:11.5D light cycle). We selected two time points across the day for tissue collection. The primary behavioral change between birds that are migrating and those that are not is seen at night when migrating birds are flying while non-migrating birds sleep. Based on preliminary studies from our lab and others 6,8 , migrating birds may sleep for a few hours at the beginning of lights out. Therefore, we chose a point in the middle of the night -Zeitgeiber time 18 (ZT18) or 18 hours after lights on -when all migratory birds were awake and displayed Zugunruhe compared to sleeping non-migratory birds. This time point captures the greatest phenotypic difference between birds in the migratory and non-migratory dispositions, but adds the potential confound of sleep status. To address this, we also collected tissue at ZT6, when both migratory and non-migratory birds are awake and feeding. Animals were decapitated and the heart and liver were quickly removed, cut into small ~5 mm diameter chunks and placed into a volume of RNAlater ™ solution (ThermoFisher Scientific, Waltham, MA) 5-10 times the tissue volume.
Tissue was allowed to stabilize for 1-2 days at 4 °C and then moved to storage at −20 °C until all samples were collected.
RNA extraction, library prep and sequencing. RNA was extracted from ~50 mg of tissue using the RNeasy mini kit (Qiagen, Hilden, Germany) according to manufacturer's instruction. DNA-digestion was performed on column, and the elution was analyzed on a Nanodrop (ThermoFisher Scientific, Waltham, MA) for concentration and Bioanalyer (Agilent Technologies, Santa Clara, CA) for RNA quality. Library preparation and sequencing was performed by the Penn State Genomics Core Facility at University Park. A barcoded library was prepared from each sample using the Illumina TruSeq Total RNA Library Prep Kit with RiboZero Gold (Illumina Inc., San Diego, CA) according to the manufacturer's protocol. Initially, 8 samples (1 sample for each experimental group, for each tissue) were sequenced on the Illumina HiSeq2500 (Illumina Inc., San Diego, CA) in Rapid Run mode using 150 nt paired-end sequencing. These data were used for a de novo assembly of the transcriptome www.nature.com/scientificreports www.nature.com/scientificreports/ in each tissue. The remaining 5 samples per group were sequenced using the HiSeq2500 in Rapid Run mode using 100 nt single-end sequencing for transcript quantification. transcriptome assembly. Raw reads were trimmed for quality and adapter content with Trimmomatic 9 using the default settings. After trimming, the sequence quality was assessed using FastQC (Babraham Bioinformatics, Babraham, England) and compared to untrimmed data. This trimmed data was used as input for Trinity 10 to assemble a de novo transcriptome for each of our tissues individually. Finally, we used TransPS 11 to post-process and annotate the de novo assemblies based on the closely related, but much more complete and well annotated zebra finch genome (Taeniopygia guttata, release 103). Differential expression and weighted gene co-expression network analysis. After transcriptome assembly and annotation, we quantified expression at the gene-level by mapping the remaining 5 biological replicates to the de novo transcriptome using Bowtie2 12 followed by eXpress 13 for transcript quantification. The resulting counts were cleaned by excluding any genes with low expression (fewer than 2 counts) and then assessed for differential expression using edgeR 14 with a 2 × 2 experimental design (time of day, migratory status, and the time x status interaction as factors) using the generalized linear model option.
As a secondary method of analysis, normalized FPKM values from eXpress were used in weighted gene co-expression network analyses 15 (WGCNA) for each tissue. Similar to the above analysis, genes with low expression were filtered out during cleaning (any gene with FPKM < 1 in any sample). Signed networks were constructed with the following parameters: deep split of 4, minimum module size of 30, and merge cut height of 0.1. www.nature.com/scientificreports www.nature.com/scientificreports/ Soft thresholding power was set to 10 for the heart networks and 9 for the liver networks based on calculated scale-free topology determined for each tissue. All other parameters remained at default values. pathway analysis. Output from differentially expressed genes or differentially expressed WGCNA modules were analyzed using Ingenuity Pathway Analysis software (IPA; Qiagen, Hilden, Germany) to determine over-represented canonical pathways and potential upstream regulators. The criteria for entry into analyses were FDR corrected p-value less than 0.1, and a 1.5-fold change in expression. For WGCNA modules, only hub genes were subjected to pathway analysis. Hub genes were defined as those that were highly correlated (r = 0.8 or higher) with module eigengene values.
statistics. All statistical tests were performed as 2-way ANOVA in R Statistics (Vienna, Austria). Factors were migratory status, zeitgeber time (ZT), and migratory × time interaction. All p-values were corrected for multiple testing using the Benjamini-Hochberg method. We considered significant effects to be those with corrected p-values of less than 0.05, whereas trends were those with p-values between 0.1 and 0.05.

Results
Determining migratory status. Previous studies from our lab have shown that the most robust and consistent measure of migratory status in individually housed birds is expression of wing-whirring behavior during the dark phase 6,8 . We quantified this behavior in 10 minutes bins across the entire dark phase, and found that birds we classified as migratory had significantly higher scores on this behavior compared to those who we classified as non-migratory (t = 7.93, df = 11.2, p < 0.00001; see Supplementary Fig. S1a). For additional verification, we also scored fat deposits, which have previously been used to determine migratory status 16  www.nature.com/scientificreports www.nature.com/scientificreports/ transcriptome assembly. Each sample was sequenced to approximately 40-million reads (raw counts 40,347,295 +/− 986,413; mean +/− SE) before being trimmed for quality and adapter content, resulting in a mean count of 33,416,401 +/− 959,794 (mean +/− SE) reads that survived trimming. Trinity assembly of heart reads resulted in 921,784 transcripts, and assembly from liver reads resulted in 930,732 transcripts. These very large transcript numbers are typical for a Trinity assembly, which include all potential isoforms of each transcript. These raw transcripts were put into TransPS to condense to the gene level and annotate against the zebra finch genome, resulting in final transcriptomes of 15,129 and 14,874 annotated gene-level transcripts for heart and liver respectively. Each assembly was evaluated by BUSCO (Benchmarking Universal Single-Copy Orthologs; a measure of transcriptome completeness) 17,18 scores. Both heart and liver achieved very high (89.7% and 90.4% respectively) completeness scores, suggesting good quality transcriptomes. Mapping reads from biological replicates was carried out with Bowtie resulting in approximately 60% mapping efficiency. Raw and processed data in the form of sequence reads and assembled transcriptomes have been deposited in GenBank and are available under the BioProject PRJNA478852 (transcriptome assembly) and GEO Series accession number GSE116989 (raw reads and quantification). Differential expression. After removing genes with low expression (FPKM < 1), edgeR was used to examine differential expression, defined as a 1.5-fold (or greater) difference between groups, and a p-value < 0.1 after FDR correction. These relatively lenient thresholds were chosen to allow the greatest amount of information to be included in subsequent pathway analysis steps. Heatmaps of differentially expressed genes (DEGs) across the four groups were constructed and are shown in Fig. 1a (heart) and Fig. 1b (liver). In the heart, 379 genes were  www.nature.com/scientificreports www.nature.com/scientificreports/ differentially expressed by time of day, 684 genes by migratory status, and 69 genes showed a time × migratory status interaction. The 10 genes with the lowest corrected p-value from each factor are presented in Table 1, and the complete list is provided in Supplementary Materials S2. In the liver, 1,757 genes met criteria for differential expression in the time of day compairson, 576 genes were differentially expressed according to migratory status and six genes showed a significant interaction between the main effects. The 10 genes with the lowest corrected p-value for each factor are presented in Table 2, and the complete list is available in Supplementary Material Table S3.
An abbreviated summary of the most significantly enriched pathways is shown in Table 3 (heart) and Table 4 (liver), and is fully enumerated in Supplemental Table SI4. For the main effect of time, pathways including circadian rhythm signaling and the unfolded protein response are significant. A large number of pathways are enriched for the main effect of migratory status, but broadly, there appear to be three main categories that overlap between both tissues: (1) metabolism, (2) immune function, and (3) cytoskeleton & cell growth.

Weighted Gene Co-Expression Network Analysis (WGCNA).
To reduce the dimensionality of the data and find closely related groups of genes that are differentially regulated by migration, time of day or the interaction between them, we used WGCNA 15 . Network construction in the heart resulted in genes grouped into 38 modules (See Supplemental Fig. S5(a) for dendrogram of grouping) and 28 modules in the liver (Dendrogram in Fig. S5(b)). Each of these modules were tested with a 2-way ANOVA to determine module significance. Overall, there were 18 modules with significant main effects of time of day, 14 with main effects of migratory status, and 10 modules were significant for the interaction in heart data. In the liver, 12 and 16 modules were significant for main effects of time and migratory status respectively. We did not find any modules with significant time × migratory status interactions in liver tissue. Genes with expression values that correlated highly with the composite measure of genes in the significant modules -the so-called hub genes -were analyzed with IPA to determine common pathways for each module. Complete summaries of these analyses can be found in Supplementary Tables S6 for heart and S7 for liver. However, we highlight several interesting modules here. Module 10 from the liver (Fig. 2b) showed significant main effects of both ZT time (F (1,16) = 18.23, q < 0.005) and migratory status (F (1,16) = 28.19, q < 0.001), but not the interaction between the two. IPA analysis revealed significant enrichment for metabolism-regulating and signaling pathways, including stearate and cholesterol biosynthesis and tryptophan degradation (Tables 5, SI7). In addition, module 35 (Fig. 2c) (Fig. 2d) from the heart was significant for migratory status (F (1,16) = 9.59, q < 0.05), but not ZT time or the interaction. This module was enriched for numerous pathways as seen in Table 5 (and SI6), but most prominent among them are several cytoskeletal and growth pathways, including Rho family signaling, actin cytoskeletal signaling, and reelin signaling.  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Behavioral changes that characterize the migratory state in birds have fascinated biologists for years. However, the physiological changes that underlie this behavioral phenotype remain poorly understood. In our study, we chose to investigate the heart and the liver because these tissues are modified during migration to meet increased cardiovascular and metabolic demands, respectively, from flight and because lack of sufficient sleep can promote cardiovascular and metabolic disorders. Recent studies have used transcriptome analysis to compare gene expression within specific tissues of migrating and non-migrating or searched for polymorphisms that might be associated with migratory subpopulations 19 , however some of the findings are inconsistent across species, as described in further detail below. Because migratory birds are a non-traditional model organism, we first constructed the tools for examining wide scale gene expression. The de novo transcriptomes that were generated for both the liver and heart appeared to be relatively complete (as measured by BUSCO scores), and the novel TransPS tool allowed us to match our white-throated sparrow sequences to well-characterized and annotated zebra finch genes. While this step precluded detection of novel or greatly divergent genes within the sparrow, it provided a higher quality annotation than would otherwise have been possible. Because the zebra finch is a non-migratory species, it is possible that by annotating via the zebra finch genome we may have missed genes completely novel for migration. However, evolution tends to be conservative, and thus it is more likely that existing genes evolved the additional function of migration, which would still be captured with the zebra finch annotation.
Unsurprisingly, we found that pathways regulating many aspects of metabolic activity were differentially enriched between migratory states, in particular for fatty acid and carbohydrate metabolism. These findings were consistent with Singh et al. 20 and Sharma et al. 21 who also found wide scale changes in gene transcripts regulating metabolic pathways for lipids, carbohydrates and protein in the black headed bunting. Previous behavioral data show that prior to migration, birds become hyperphagic and increase fat stores, and that failure to do so restricts the ability to enter the migratory state 22,23 . This accreted fat is then used as an energy source during migration (Reviewed in Guglielmo, 2018) 24 , and our data are consistent with increased fatty acid utilization. However, in addition to differences in fat metabolism, we found enzymes responsible for amino acid metabolism, such as tyrosine aminotransferase, and nucleic acid salvage, such as uridine phosphorylase 1 & 2, were greatly upregulated during migration. These pathways present a mechanism by which tissues within the body can be catabolized for energy to fuel the oversized energetic demands of long-distance flight. Module 10 in the liver was heavily enriched for pathways controlling metabolism of carbohydrates, nucleic acids, lipids, and amino acids. Fudickar et al. similarly found differences in genes responsible for fatty acid and carbohydrate metabolism between resident and migratory dark eyed juncos 25 . These processes have implications not only for fueling flight, but also for influencing other uncovered pathways. For example, degradation of the amino acid tryptophan can induce the kyneurin signaling pathway, a process leading to immune tolerance and further changes in metabolism such as increased NAD+ synthesis 26 . NAD+ is an important co-factor for energy generation via glycolysis and lipolysis 27 . Interestingly, WGCNA module 5 in the liver shows that pathways for ketogenesis are enriched during migration, suggesting that birds may be fueling their journey, in part, by producing and burning ketones, most likely derived from fatty acid stores. Conversion of fatty acids to ketones is important as the neurons in the brain cannot  www.nature.com/scientificreports www.nature.com/scientificreports/ metabolize fatty acids directly 28 and the inflammasome of the immune system can be downregulated by the main ketone metabolite β-hydroxybutyrate 29 .
In the general category of cell signaling, we see several interesting pathways associated with the migratory phenotype. Although the literature on signaling changes that occur during migration is limited, our data are consistent with alterations in glucocorticoid signaling 30 and adiponectin signaling 8 . Interestingly, we find previously unreported alterations in downstream pathways involved in mediating the effects of thyroid hormone, retinoic acid, mineralcorticoids and prolactin. Migratory behavior in Swainson's thrushes was also strongly correlated with differential expression of gene transcripts involved in retinoic acid and thyroid hormone signaling in the brain 31 . These findings affirm our transcriptome results, in the case of glucocorticoids and adiponectin, and also provide novel hormonal pathways to explore in regards to the regulation of migratory behavior. Another system in which changes are known to occur during migration is the immune system. For example, innate immune response to a challenge with phytohemagglutinin is reduced during migration 32 . Similarly, our data showed differential gene expression and immune-related pathway enrichment during migration. In particular, module 35 in the heart demonstrates broad changes in immune signaling during migration (Fig. 2c). The immune system can also modulate the efficiency of whole body metabolism by regulating insulin sensitivity at the site of target tissues 33 . The immune system requires significant amounts of energy to be maintained in an active state, so it may be advantageous to reduce immune function during this period of intense energy demand for flight. In general, this hypothesis is borne out in our data, where we observed widespread downregulation of proinflammatory cytokine pathways. However, there was one notable except to this overall downregulation of the immune system: those portions of the proinflammatory pathways geared towards increasing levels of scavenging and remodeling of damaged tissues appear to be upregulated. These processes are similar to those found during and after successful recovery from tissue damage due to insult, such as myocardial infarct 34 . For example, some athletes involved in extreme exercise programs may experience either "athlete's hepatitis" 35 or cardiac myopathies as a result of either overuse or repeated transient ischemia 36 . We propose that the immune system is modified in a novel way during migration to help remove and repair the physical damage to these tissues caused by overexertion during long distance flight. Additionally, it is intriguing that many gene pathways associated with the formation of fibrosis are significantly downregulated during migration, as fibrosis formation often occurs during tissue damage as part of  NEDD4, PSMA6, UBE4B, CUL1, DNAJC3, HSPA5, DNAJA1, SKP2,  HSPA8, PSMD11, HSP90B1, PSMC1, PSMC6, USP47, PSMB1, HSP90AA1 www.nature.com/scientificreports www.nature.com/scientificreports/ the wound repair process 37 . These findings further suggest that part of the migratory phenotype is the utilization of the immune system to efficiently repair and remodel tissue as it is catabolized to fuel flight. Franchini et al. (2017) found that the receptor for the cytokine TGF-beta is differentially expressed in the blood of migrant and resident European blackbirds, further suggesting that investigating the role of the immune system in migration may be useful in uncovering systemic regulation of these seasonal events 38 .
One unexpected finding was the general downregulation, in the migratory group, of transcripts for cellular matrix and adhesion proteins. This group of transcripts comprised roughly 6% and 6.5% of all of the differentially expressed genes in the liver and heart, respectively. When analyzed with WGCNA, we found that module 14 in the heart was particularly enriched with pathways controlling fibrosis, cytoskeleton rearrangement, and integrin regulation. These proteins are largely responsible for maintaining the extracellular spaces and providing scaffolding for the attachment of transmembrane proteins on the cell surface 39 ; our findings include changes in transcripts for molecules such as cadherins, integrins, glycoproteins and multiple types of collagens. Similarly, Johnston et al. 31 and Sharma et al. 21 found general downregulation of genes involved in cell adhesion and cell motility, consistent with our results 31 . We hypothesize that this downregulation of extracellular matrix proteins may help modulate cell to cell communication, cell adhesion and cellular differentiation. Regulation of extracellular matrix composition may be another mechanism that birds use to facilitate repair of tissues damaged from overuse during long distance flight in a way that does not trigger the traditional wound repair mechanisms that could lead to fibrosis or scarring. This coordinated mechanism for organ rejuvenation programed within the migratory phenotype would be necessary to avoid fibrosis or scarring since birds engage in migration twice per year and have a lifespan of 10 years or more, at least in captivity. Although the life expectancy of songbirds in the wild is only a few years, due to predation and other risks, our lab currently houses several white-throated sparrows over 11 years old.
The main focus of this study was to examine changes in RNA levels that occur during migration. One interesting aspect that we attempted to capture in our data is that while both migratory and non-migratory birds share a similar phenotype during the day (i.e. awake and foraging for food), dramatic behavioral differences are seen at night (i.e. sleep when birds are non-migratory vs. awake, navigating and flying when migrating) 6 . To examine this behavioral difference, we included two time points (one during the day and one at night) for transcriptome analysis. Consistent with previous literature, transcripts and pathways known to be a part of the core molecular circadian clock or known to be directly regulated by the molecular circadian clock are significantly different between our two time points (Reviewed in Cassone et al. 40 ). Circadian clock and clock target genes were also found to be differentially expressed with time of day in the hypothalamus of migratory Swainson's thrushes and the liver of black headed buntings 20 . Some of our most interesting findings are in the interaction between time of day and migratory status, which implies some aspects of clock control differ with migratory phenotype, unlike what was observed in Swainson's thrushes. Because, the number of samples required for sufficient power to detect significant gene-level interactions is high, we likely underestimated the number of genes affected by this interaction. Nevertheless, we found 10 WGCNA modules within the heart that showed significant interactive effects. These data suggest that there are suites of genes that are differentially expressed while migrating birds are awake at night, and that these are distinct from those expressed during the day in either the migratory or non-migratory states. These suites of genes, in particular, may prove the most useful in dissecting candidates for investigating clock control of migration in future work.  www.nature.com/scientificreports www.nature.com/scientificreports/ In summary, the work presented here provides new insights into, and key understanding of, the endogenously generated seasonal changes in the transcriptome that occur during nocturnal migration in birds. These physiological changes include increased endurance and the ability to overcome the typically deleterious effects of sleep loss. Our data are consistent with previous findings such as increased rates and efficiency of lipid metabolism. However, we were able to extend those findings to elucidate novel molecular mechanisms, such as the transcriptional control of enzymes responsible for regulation of different metabolic pathways. Also, key amongst the morphological and physiological changes during migration are those regulated by molecular pathways involved in direct cell-to-cell signaling and adhesion, such as actin and integrins. These pathways, coupled with the extensive downregulation of many immune-related pathways, are unique. They simultaneously maintain the tissue-repair and remodeling machinery at a high efficiency for both catabolism of tissue for fueling flight, fixing damage to tissue through prolonged use and ultimately providing a scaffolding network for tissue regrowth after migration, while reducing the energetic requirements of the full immune system. Our work suggests that these processes are undertaken via previously unreported mechanisms. Taken together, our data will allow for the generation of multiple novel hypotheses on the physiological regulation of migration and will lead to a more complete understanding of the unique ability of birds to dynamically reduce sleep requirements while completing ultra-endurance events without concomitant negative consequences.