Transcriptome analysis of a wild bird reveals physiological responses to the urban environment

Identifying the molecular basis of environmentally induced phenotypic variation presents exciting opportunities for furthering our understanding of how ecological processes and the environment can shape the phenotype. Urban and rural environments present free-living organisms with different challenges and opportunities, which have marked consequences for the phenotype, yet little is known about responses at the molecular level. We characterised transcriptomes from an urban and a rural population of great tits Parus major, demonstrating striking differences in gene expression profiles in both blood and liver tissues. Differentially expressed genes had functions related to immune and inflammatory responses, detoxification, protection against oxidative stress, lipid metabolism, and regulation of gene expression. Many genes linked to stress responses were expressed at higher levels in the urban birds, in accordance with our prediction that urban animals are exposed to greater environmental stress. This is one of the first studies to reveal transcriptional differences between urban- and rural-dwelling animals and suggests an important role for epigenetics in mediating environmentally induced physiological variation. The study provides valuable resources for developing further in-depth studies of the mechanisms driving phenotypic variation in the urban context at larger spatial and temporal scales.


Results and Discussion
Urban and rural birds exhibit different gene expression profiles. Differential gene expression analyses, together with a principal components analysis (PCA; Supplementary Fig. S1), revealed that the urban-and rural-dwelling great tits exhibited significantly different expression profiles ( Fig. 1; Supplementary Data S2 and S3). This pattern was evident in both liver and blood tissues: 304 and 372 genes were significantly differentially expressed (q < 0.1) between the two environments in liver and blood, respectively ( Fig. 1; Supplementary Data S3). Supplementary Tables S1 and S2 show the top 20 most significant annotated genes that were differentially expressed in each of liver and blood tissues. Of those genes that were differentially expressed, the majority were expressed at higher levels in the urban (liver: 238 [78%]; blood: 321 [86%]), relative to the rural, environment. Neither the heatmaps (Fig. 1) nor PCA ( Supplementary Fig. S1) revealed high similarities among the two different age classes (see Methods), and so we have no reason to believe that differences in age are driving the observed variation in gene expression between urban and rural individuals. Gene enrichment analysis revealed significant overrepresentation of 318 and 281 GO terms assigned to the differentially expressed genes in liver and blood transcriptomes, respectively (Supplementary Data S4). While only 9% (n = 53) of the differentially expressed genes were common to liver and blood ( Supplementary Fig. S2A), 28% (n = 130) of the significantly overrepresented GO terms were common to both tissues ( Supplementary Fig. S2B). This result indicates that although some of the environmentally induced functional responses operate in both liver and blood tissue, they largely involve different genes. Differences in expression profiles are expected, given the different functions of the two tissues, and emphasise the value of investigating transcriptional responses in other tissues besides blood. Many of the genes that were more highly expressed in urban, compared with rural, birds are involved in physiological and biochemical responses to the biotic and abiotic stressors that are characteristic of urban environments. In the proceeding sections, we evaluate and discuss the functional significance of the observed variation in gene expression. It should be noted that gene expression at the RNA level might not always correlate to actual protein levels 30 , though several studies have shown good correlation between mRNA and protein concentrations 31,32 . Urban upregulation of adaptive immune responses. Genes involved in the development, activation and regulation of adaptive immune responses were significantly overrepresented among the differentially expressed genes from both liver and blood ( Fig. 2A,B; Supplementary Data S4). Indeed, many of the most significant differentially expressed genes are involved in functions related to adaptive or innate immunity (Supplementary Tables S1 and S2). Overrepresented functions primarily relate to activation and differentiation of T-cells (antigen receptors) and B-cells (involved in antibody production; Fig. 2A), antigen presentation and processing (e.g. via dendritic cells and MHC; Fig. 2A), and their relevant signalling pathways ( Fig. 2A,B). T-cells and B-cells are the key mediators of adaptive immunity, and when presented with an antigen or a signal from microbes, proliferation and differentiation of lymphocytes proceeds 33 . Remarkably, nearly all genes associated with the GO terms 'immune system process' (liver: n = 53/53; blood: n = 57/59) and 'immune response' (liver: n = 29/29; blood: 29/30) were expressed at higher levels in the urban, relative to the rural, birds (Supplementary Data S3). Urban birds also demonstrated significantly higher expression of numerous genes whose products are involved in the detection of, and responses to, bacteria and compounds originating from bacteria (Fig. 2B). These differences were more evident from blood (14 genes associated with 13 significant GO terms), compared with liver (6 genes associated with 3 significant GO terms), transcriptomes (Supplementary Data S3 and S4). Since the circulatory system provides a fundamental role in the body's defence against bacterial invasion, environmentally induced changes in bacterial defences would be expected to be more pronounced in the blood, compared with other tissues.
Scientific RepoRts | 7:44180 | DOI: 10.1038/srep44180 Consistent higher expression of genes linked to immune functions in urban birds could be explained by higher pathogen diversity and/or prevalence -and thus increased frequency of exposure to microbes and antigens -in the urban, compared with the rural, environment. It has been shown that urban environments can influence host-pathogen interactions, via changes in prevalence, routes of transmission and susceptibility of wild animals to disease 29,34 . Yet, while some studies report higher rates of disease and parasitism in urban, compared with rural, bird populations 5,35 , others report reduced pathogen prevalence and intensity in urban birds 4 . If pathogen diversity and/or prevalence are indeed higher in urban environments, increased activation of immune responses may increase the ability of urban great tits to cope with urban-associated changes in host-pathogen dynamics. Indeed, a previous comparative analysis of 39 phylogenetically-independent events of urbanisation suggested that species that showed stronger immune responses (based on the size of the bursa of Fabricius) were more successful in invading urban areas than species with weaker immune responses 36 . However, enhanced immune activity may not necessarily mean birds can avoid the expected negative impacts on survival or reproduction associated with higher pathogen diversity and prevalence, since there may also be fitness costs associated with an upregulation of immune defences 37 . Regardless of the mechanisms involved, the transcriptome data from our urban and rural study populations suggest that characteristics of the adaptive immune system could be central to the underlying differences between urban-and rural-dwelling great tits.
Urban birds have enhanced innate immune activity. Among the genes that were significantly differentially expressed between the urban and rural birds, GO terms linked to inflammatory response pathways were significantly overrepresented in both liver and blood transcriptomes ( Fig. 2B-D). More specifically, overrepresented processes primarily relate to the production, secretion and receptor-binding of cytokines, and chemotaxis (liver: n = 41 genes associated with 29 significant GO terms; blood: n = 24 genes associated with 6 significant GO terms; Supplementary Data S3). Cytokines (e.g. tumour necrosis factors, chemokines, and interleukins) play an important role in host response to infection and trauma, by mediating inflammatory responses 33,38 typified by redness, swelling, heat and pain 39 . Similar to our findings concerning the adaptive immune system, all but one . Individuals are denoted U1-U6 and R1-R6 for urban and rural individuals, respectively. Normalised gene expression has been variance-stabilised, and scaled and centred around zero to produce Z-scores. A darker colour indicates higher expression and a lighter colour indicates lower expression. Scatterplots show mean normalised gene expression levels (log 10 -transformed) in (C) liver (n = 12 individuals) and (D) blood (n = 11 individuals). Red points indicate significant differentially expressed genes (q < 0.1), whereas black points indicate genes that are not significantly differentially expressed. The blood transcriptome from individual R2 was removed from the analysis since it was identified as an outlier (see Supplementary Fig. S4).  of the genes associated with overrepresented inflammatory processes were expressed at higher levels in urban, relative to rural, birds (Supplementary Data S2 and S3). Increased expression of cytokine-encoding genes has been widely demonstrated in response to bacterial 40 , viral 41 and parasitic 42 infections in avian macrophage cells. The CP gene, coding for ceruloplasmin, is just one example of an inflammatory-linked gene that is expressed at higher levels in the urban birds; as well as playing a key role in iron detoxification (see below), ceruloplasmin is an acute phase protein that is synthesised in the liver during the acute phase of an inflammatory response and, as such, is an indicator of acute inflammation in birds 43 . Markers of inflammation are widely used in toxicological and biomedical studies, and enhanced inflammatory responses in humans have, in addition to being linked with infection and trauma, been associated with exposure to pollutants, such as heavy metals 44 , particulate matter 45 and nitrogen oxides 46 . The observed differences in gene expression linked to innate immune activity between the urban and rural birds could therefore be explained by differences in exposure to pathogens, traffic-generated particulate matter, noise and/or artificial night light; increased exposure to all of these sources of environmental stress are likely to cause upregulation of inflammatory responses 12,47 .
While inflammatory responses are part of the normal innate immune response, excessive or uncontrolled inflammation contributes to tissue damage and disease 38 . Inflammation triggers an increased production of reactive oxygen species (ROS), which, if not effectively removed by antioxidant defences, leads to oxidative stress and damage to cellular macromolecules such as nucleic acids, lipids and proteins 48 . Thus, inflammation and oxidative stress are tightly linked. Increased oxidative stress is widely implicated in senescent declines in organismal function and the pathology of many cancers and cardiovascular diseases 48 . Although we have some understanding of how oxidative stress is affected by urban pollution, we know little, if anything, about inflammatory responses to urban stressors in wild animals 12 . While enhanced inflammation could be beneficial, providing defences to infection, injury and pollution associated with urban environments, persistent high activity of innate immune responses could also be degenerative for urban birds.
Urban birds have more active detoxification and repair systems. Although the GO term 'metal ion binding' was not significantly overrepresented, genes associated with this term were numerous among the differentially expressed genes in both liver (n = 38) and blood (n = 63; Supplementary Data S3). Furthermore, the majority of genes were expressed at higher levels in the urban, compared with the rural, birds (liver: n = 23/38; blood: n = 54/63). Two genes coding for metallothioneins (B5G2T6_TAEGU and MT4) were not only expressed at significantly higher levels in the livers of urban, relative to rural, great tits, but were among the most significant of all the differentially expressed genes from liver transcriptomes (Supplementary Table S1). Urban great tits also exhibited significantly higher expression, in blood tissue, of a metalloprotein gene (CP) and three genes coding for metallopeptidases (MME, MMP28, ADAM19; Supplementary Data S3); the primary biological functions of all of the encoded proteins concern metal ion binding and metal detoxification.
Although metals are essential micronutrients for organisms, playing key roles in metabolic processes and signalling pathways, they can become toxic at higher concentrations 25 , partly as a consequence of oxidative stress 49 . The environmental toxicological literature abounds with studies demonstrating oxidative stress induced by exposure to pollutants in both terrestrial and aquatic environments 45,50 . While there is evidence for pollutant-induced oxidative stress in urban birds, our knowledge concerning the impacts and capacity for modulating detoxification pathways in birds is limited 28,51 . Our results suggest that the expression of genes involved in metal detoxification pathways could be important in underlying the functional differences between urban and rural birds. Metallothioneins (MT) play an important role in the detoxification of heavy metals 52 and have been shown to be efficient scavengers of hydroxyl-radicals in vitro 53 . Concentrations of MT in liver and kidney have been shown to strongly correlate with levels of cadmium in great tits 24 , while exposure to a number of different metal ions has been shown to result in enhanced activity of various antioxidant enzymes and proteins, including MT, in fish 54 . Similarly, in mice, MT gene expression is upregulated in response to oxidative-stress-inducing agents, including metal ions, xenobiotics, and inflammation, as well as oxidative stress itself 52,55 . The observed high expression of genes involved in metal detoxification pathways in urban birds, in the present study, could represent a response to limit the damaging effects of elevated levels of pro-oxidants and reduce the toxic effects arising from exposure to heavy metals in the urban environment.
Among the genes that were significantly differentially expressed between urban and rural birds, there was significant overrepresentation of genes involved in oxidation-reduction activity (i.e. associated with the GO term 'oxidoreductase activity' and all directly descending terms) in both blood (n = 30; Fig. 2C) and liver (n = 3; Fig. 2C) tissues. For example, the urban birds exhibited higher expression of the TXNRD1 gene, coding for thioredoxin reductase 1, responsible for the catalysis of the reduction of thioredoxin, which is an important cellular antioxidant. This suggests that the urban birds may have higher circulating levels of thioredoxin, which could be a response details. Abbreviated terms: (a) myeloid cell activation involved in immune response; (b) macrophage activation involved in immune response; (c) antigen processing and presentation of peptide or polysaccharide antigen via MHC class II; (d) antigen processing and presentation of peptide antigen via MHC class II; (e) adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains; (f) immune response-regulating cell surface receptor signaling pathway; (g) detection of mechanical stimulus involved in sensory perception; (h) oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular oxygen to two molecules of water; (i) oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen; (j) positive regulation of sequence-specific DNA binding transcription factor activity; (k) regulation of sequence-specific DNA binding transcription factor activity; −ve = negative; +ve = positive.
Scientific RepoRts | 7:44180 | DOI: 10.1038/srep44180 -either adaptive or non-adaptive -to exposure to high levels of oxidative stress. In addition, blood transcriptomes of the urban great tits revealed significantly higher expression of two genes encoding nucleases (EXO1 and FEN1) that play key roles in DNA repair 56 . This result provides a further indication that urban birds could be exposed to higher levels of oxidation-induced DNA damage. A recent study from the same populations confirms this finding, by experimentally demonstrating that being reared in the urban environment leads to shorter telomeres 57 . Short telomeres have been widely associated with reduced lifespan and fitness, and oxidative stress is likely a key mediator of accelerated telomere attrition 58 . Taken together, our results suggest that differences in expression levels of genes related to mechanisms involved in detoxification and repair of oxidative damage correspond with the likely differences in exposure to oxidation-inducing chemicals between urban and rural individuals. Differential regulation of fatty acid physiology. Several genes that were significantly differentially expressed between the urban and rural birds are involved in lipid metabolism, binding and storage in both liver and blood (Supplementary Data S3). In particular, liver transcriptomes revealed significant overrepresentation of GO terms linked to the synthesis and elongation of fatty acids (n = 2 genes associated with 6 significant GO terms; Fig. 2C; Supplementary Data S4). The fatty acid composition of tissues is partly determined by the diet (which is likely to differ between urban and rural habitats), and has important implications for many physiological functions and processes including thermoregulation 59 , inflammation and oxidative stress 39 . These processes themselves are sensitive to extrinsic stressors, such as pollutants, associated with the urban environment 28,45,50 .
The rural great tits exhibited significantly higher expression, in the liver, of enzymes -two elongases (ELOVL2, ELOVL5) and a desaturase (FADS6) -involved in the biosynthesis of long-chain polyunsaturated fatty acids 60,61 (PUFAs; Supplementary Data S3). Long-chain PUFAs are prone to lipid peroxidation 62 , thus having higher levels increases an organism's susceptibility to oxidative damage upon exposure to pro-oxidants. As discussed above, our transcriptome data -in line with other studies -suggest that urban birds already face a greater oxidative threat; maintaining low expression of these elongases in the urban environment could therefore be a response to reduce the extent of oxidative damage to membrane lipids.
PUFAs are converted -via elongation and desaturation pathways -to eicosanoids 62 , which modulate the intensity and duration of inflammatory responses 39 . Typically, eicosanoids derived from the omega(ɷ)-6 PUFA pathway have pro-inflammatory effects (e.g. fever, vasodilation, oedema), while those derived from ɷ-3 PUFAs, have anti-inflammatory effects 39 . Although inflammatory responses act to combat harmful stimuli, excessive and uncontrolled inflammation can incur costs as a result of the generation of pro-oxidants and subsequent damage to macromolecules 48 . Interestingly, it was recently shown in our study populations that urban great tits had higher circulating levels of the ɷ-6 PUFA arachidonic acid, while rural great tits had higher levels of the ɷ-3 PUFA eicosapentaenoic acid 21 , suggesting that urban birds have higher levels of pro-inflammatory eicosanoid precursors. Although the relative preferences of the two elongases (ELOVL2 and ELOVL5) for substrates from the ɷ -6 and ɷ -3 pathways 61 have not been determined in birds (for mammalian studies, see 60,61), reducing expression of the genes in the urban environment may be a response to limit endogenous synthesis of pro-inflammatory PUFAs and a further increase in inflammation and subsequent oxidative damage.

Epigenetic mechanisms as a possible mediator of responses to urban stressors. Expression
profiles from liver and blood tissues revealed significant overrepresentation of several GO terms linked to the regulation of gene expression, via epigenetic modifications such as DNA methylation and histone methylation/ acetylation, and via regulation of transcription factors ( Fig. 2D; Supplementary Data S4). Furthermore, the specific genes concerned (liver: n = 15 genes; blood: n = 11 genes; Supplementary Data S3) are linked to DNA damage repair, tumour suppression, senescence and metal-ion binding, all of which relate to other observed transcriptional differences between the urban and rural great tits. Epigenetic mechanisms act to dynamically regulate the organisation and function of the genome, by regulating transcription and thus gene expression 63 . One of the most significant differentially expressed genes from liver transcriptomes encodes a methyltransferase (MTR; Supplementary Table S1), which catalyses the addition of methyl groups to DNA. DNA methylation is the most well-understood epigenetic mark and plays a crucial role in regulating gene expression. Furthermore, MTR was expressed at higher levels in rural, compared with urban, birds, suggesting that DNA-methylation levels of rural birds could be higher than those of urban birds in the study populations. Since methylation is typically associated with reduced gene expression 63,64 , higher expression of this methyltransferase gene in rural birds could contribute to overall lower levels of gene expression observed among the significant differentially expressed genes in these individuals (Fig. 1).
The results indicate a potential key role for epigenetic mechanisms in mediating the observed variation in gene expression, and likely subsequent phenotypic variation, between urban and rural birds. Indeed, recent sequencing of the great tit methylome indicated an important role for methylation in evolutionary processes within the species 65 , while methylation has also been shown to be important in shaping great tit personality 66 . The epigenome is highly responsive to environmental cues and aberrant patterns of DNA methylation have been linked with exposure to a wide range of environmental contaminants in humans and other animals 64,67,68 , as well as cardiovascular diseases, cancers and ageing 69,70 in humans. Studies on laboratory rodents and humans have shown that patterns of methylation are particularly sensitive to early-life nutrition and may underpin disease susceptibility in later life 70,71 . More recently, diet was also shown to influence whole-genome methylation in wild baboons (Papio cynocephalus) 72 . The inferred differences in levels of DNA methylation between the urban and rural birds could therefore be a result of differences in pollutant exposure and/or diet between the two different environments. Importantly, since epigenetic modifications can be inherited during mitotic cell division, environmentally induced changes in gene expression can persist throughout life, resulting in life-long changes in phenotypic traits. Very little is known about epigenetic mechanisms in birds 73 , how patterns of methylation in birds are influenced by the environment, and whether epigenetic modifications are primarily the result of "developmental programming".
Conclusions. This is one of the first studies to quantify variation in genome-wide gene expression of a wild animal in the context of urbanisation. The existence of consistent and marked differences in gene expression profiles between urban and rural great tits, in the present study, suggests the regulation of responses at the molecular level to local environmental conditions. While the present study is limited by considering only a single pair of populations, many of the upregulated genes and functional differences were in the predicted direction and in accordance with previous physiological studies. For example, many genes linked to stress responses were expressed at higher levels in urban, compared with rural, birds, suggesting that the observed differences in gene expression are a result of differential stress-exposure in the two contrasting environments. The results suggest that differences in diet and exposure to pollutants and pathogens could all contribute to the observed variation in transcriptomes between urban and rural great tits. Future experimental studies should aim to determine the specific role of these different stressors and associated fitness consequences.
By offering novel insight into the biochemical and physiological functions involved in mediating responses to urban environmental stressors, our findings will facilitate identification of key regulatory genes in both liver and blood tissues. Since current sequencing costs restrict the viability of performing replicated transcriptome analyses across multiple pairs of urban/rural populations (while also retaining sufficient replication at the individual and tissue level), our results should be used to inform the application of candidate gene approaches to verify the consistency of the observed responses across broader spatial and temporal scales. It remains to be understood to what extent the observed transcriptomic differences between the urban and rural birds are the result of genetic adaptation, epigenetic modifications, and/or simply direct responses to ambient stressors. Our results, however, indicate that epigenetic mechanisms could play a key role in mediating the environmentally induced variation in gene expression in the context of urbanisation. Importantly, epigenetic mechanisms could generate life-long changes in gene expression and subsequently phenotypic traits.

Methods
Ethics statement. This study was carried out in accordance with Swedish legislation and approved by the Malmö-Lund animal ethics committee (Dnr M454 12:1).
Study species and sites. Twelve wild male great tits -six from an urban environment and six from a rural environment -were euthanised in late winter 2014 in southern Sweden. Urban birds (denoted as U1-U6) were captured in Malmö (55°35′ N 12°59′ E; Supplementary Fig. S3A) -the third largest city in Sweden with c. 300,000 inhabitants. The study site in Malmö is an urban park that is characterised by a mix of coniferous and deciduous trees, managed grassland, ponds and associated urban infrastructure including paved footpaths and buildings. Rural birds (denoted as R1-R6) were captured in Vombs fure (55°39′ N 13°33′ E; Supplementary Fig. S3B), a rural site located 35 km northeast of Malmö. Vombs fure comprises ~12 km 2 of, primarily coniferous, forest interspersed with broadleaved trees; the surrounding area is sparsely inhabited by humans (< 5 inhabitants km −2 ). For further details of study sites, refer to 21 . Rural birds were captured seven days after urban birds. Since breeding occurs up to 7 days earlier in Malmö, compared with Vomb, we do not expect there to be differences between urban and rural birds with respect to physiological changes in anticipation of the breeding season. Birds were captured while roosting in nestboxes at night and transported to the laboratory in Lund where biometrics were recorded, and a blood sample was collected and immediately transferred to storage at − 80 °C. Post-mortems were carried out immediately after sacrifice; liver tissue was collected and transferred to storage at − 80 °C within 5 min of death. Rural birds (mean ± SE: 19.5 ± 0.25 g) were significantly heavier than urban birds (mean ± SE: 17.8 ± 0.22 g; T-test: t = 3.46, df = 6.8, P = 0.011), though there were no differences in body size, as measured by tarsus length (mean ± SE: rural: 22.9 ± 0.20, urban: 22.6 ± 0.79; T-test: t = 0.35, df = 4.5, P = 0.739). All birds were adults: rural birds were all aged 3 K + (i.e. 3 yr or older); the urban birds U1-3 were 2 K (i.e. 2nd calendar year), while U4-6 were 3 K+ .
RNA isolation, library preparation and sequencing. Total RNA was isolated from 25 mg liver tissue and 15 μ l of whole blood for each of the 12 individuals. Samples were first homogenised in 60 μ l (liver)/20 μ l (blood) Autoclaved Millipore water using a TissueTearor (BioSpec, Bartlesville, USA). A single Trizol/chloroform extraction step was performed, followed by purification using Qiagen's RNeasy kit (Hilden, Germany), according to the manufacturer's instructions. RNA yield and purity were checked using a NanoDrop spectrophotometer (Thermo Scientific, USA) and an Agilent 2100 Bioanalyser (Santa Clara, USA). All samples were confirmed to be of acceptable quality as indicated by the 260/280 nm absorption ratio (median ± SD: 2.02 ± 0.05) and RIN (RNA Integrity Number; median ± SD: 7.1 ± 1.1). DNase treatment was performed prior to sequencing. Truseq library construction and 100 bp paired-end RNA-seq of individual samples was carried out using Illumina HiSeq 2000 by BGI Tech Solutions (Hong Kong). All samples were sequenced twice, in different lanes, and a mix of urban and rural samples were run in each lane. For reasons beyond our control, the RNA from the blood of R2 was sequenced separately from all other blood RNA samples; since the resulting transcriptome was highly dissimilar to all other transcriptomes ( Supplementary Fig. 4), this individual was removed from the analysis of blood transcriptomes.
Transcriptome analysis. Raw RNA-seq reads were multiplexed, and subjected to adaptor trimming and quality screening. We obtained 763 million and 772 million high-quality paired-end reads, from liver and blood tissues, respectively. Reads (1.5 billion; 767 million pairs) were mapped using TopHat2 2.0.12 74  genome has recently been published 65 , within the passerines, the zebra finch has undoubtedly the highest-quality genome assembly, with 18,618 annotated genes compared with 13,036 genes in the great tit genome. Furthermore, GO information exists for 14,601 genes in the zebra finch genome, while GO data does not yet exist for the great tit. Avian genomes are highly conserved between species, making the use of a well-annotated genome from another species possible 75 . A maximum mismatch rate of 25% was allowed to account for species divergence. A total of 582 million (76.3%) and 632 million (81.9%) reads, from liver and blood tissues respectively, successfully aligned to the genome. Unique reads that mapped unambiguously to genes were counted using HTSeq 0.5.3p9 76 and SAMtools 0.1.19 77 . All raw reads have been deposited in NCBI's Sequence Read Archive (PRJNA314210).
Differential gene expression analysis. Differential gene expression analyses were performed using DESeq2 1.6.3 78 . As per Love et al. 78 , individual counts were normalised using the geometric mean and modelled with a negative binomial distribution. Shrinkage of dispersion and fold-change estimates are implemented in DESeq2 to improve stability and interpretability, meaning that logarithmic fold-changes will have more shrinkage when there is little information available for a gene (i.e. low read count, high dispersion, or few degrees of freedom). P-values were corrected for false positives using the Benjamini and Hochberg false discovery rate (FDR) correction 79 for multiple testing, within DESeq2. Genes with an FDR (or q-value) < 0.1 were regarded as statistically significantly different. Hierarchical clustering, using Principal Component Analysis (PCA), was performed on expression data from liver and blood tissues ( Supplementary Fig. S1A-C). Expression values were transformed using the recommended variance stabilizing transformation procedure for use in PCA and to generate heatmaps (Fig. 1A,B). One individual (R2) was removed from the analysis of blood transcriptomes, since the PCA and distance heatmaps revealed it to be an extreme outlier ( Supplementary Fig. S4), possibly due to it being sequenced in a separate sequencing run.
Gene enrichment analysis. Functional annotations and overrepresentation of GO terms were analysed using the Cytoscape plug-in BiNGO 3.0.2 80 . A full ontology file for the zebra finch was downloaded from www.geneontology.org (08- . Gene enrichment analyses were performed using a hypergeometric test for overrepresentation of GO terms from the differentially expressed gene set, using the zebra finch GO as the background. Again, P-values were corrected using the Benjamini and Hochberg FDR, within BiNGO, and GO terms with a q-value < 0.1 were regarded as statistically significantly different. We evaluated GO terms from the biological process and molecular function domains. For the production of heatmaps in Fig. 2, GO term sets were condensed using REVIGO 81 . Data accessibility: RNA-seq raw reads are deposited in NCBI's Sequence Read Archive (PRJNA314210).