Using the canary genome to decipher the evolution of hormone-sensitive gene regulation in seasonal singing birds

M1. Animals, testosterone levels, tissue sampling, and RNA isolation

We used zebra finches and common domesticated canaries bred at the Seewiesen animal facility. For the genome sequencing, a female canary was used; for all other procedures, male canaries and zebra finches were used (Additional file 1). Animal handling was carried out in accordance with the European Communities Council Directive 2010/63 EU and the legislation of the state of Upper Bavaria (record number 55.2-1-54-2532-68-12). Reproductively active adult males (one-year-old canaries, n = 21; zebra finches, n =9) were housed pairwise with a female under long-day photoperiod (LD; 14 hours light:10 hours dark). Outside of the breeding season, male canaries (n = 15) were housed pairwise with a female under short-day photoperiod (SD; 10 hours light:14 hours dark). Songs were recorded before killing to confirm that the males had adult typical songs. Six canaries were treated with testosterone. Testicle weight (mean ± standard deviation) was 158 ± 27 mg for the LD canaries, 1.2 ± 0.5 mg for the SD canaries, 0.9 ± 0.6 mg for the SD canaries treated with testosterone, and 58 ± 14 mg for the zebra finches. These testicle weights are typical for reproductively active adult zebra finches and canaries and canaries during the non-breeding season. For testosterone treatment, the birds were subcutaneously implanted with custom-made silastic implants filled with testosterone, as described previously [64].

Testosterone levels (mean ± standard deviation) of the LD males, SD males, and SD + T males were 2.84 ± 1.11 ng/ml plasma, 0.09 ± 0.05 ng/ml plasma, and 6.29 ± 3.79 ng/ml plasma, respectively, on the day of killing. A testosterone radioimmunoassay was carried out as described before for canaries [16].

For RNA-seq, RNA of the HVC, RA, or ENT of three males was pooled for each biological replicate; that is, we used a total of nine male zebra finches and nine male canaries. The other males were used for microarray procedures. The birds were sacrificed quickly with an overdose of isoflurane. Afterwards the brains were quickly removed and snap-frozen over liquid nitrogen. For RNA preparation, brains were cut into four 40 μm para-sagittal sections, mounted each on one slide, followed by two 20 μm sections on a cryostat, mounted on RNase free Superfrost slides as two different series. One series of 20 μm sections was used for Nissl-stainings. These Nissl-stainings were used to verify the location of the HVC, RA, and ENT and to guide the dissection of these areas. The areas were dissected from the 40 μm sections with micro-scalpels (a different one for each brain area to avoid contamination) under a stereo-microscope, transferred into Qiazol (Qiagen Inc., Valencia, CA, USA), and stored at −80°C until RNA preparation.

For in situ hybridizations, the left half of each brain of 6 male canaries was cut into 10 adjacent series of 15 μm sections and mounted on RNase free Superfrost slides. One series was used for Nissl-staining; the others were stored at −80°C for in situ hybridizations with various probes.

RNA was extracted by the automated Qiacube system using a RNeasy Micro Kit according to the manufacturer’s protocol (Qiagen Inc., Valencia, CA, USA ) using the optional DNA digest step. RNA quality and concentration were assessed using an Agilent Model 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and Nanodrop 1000 spectrometer (Thermo Fisher Scientific, Wilmington, MA, USA), respectively.

M2. Whole genome shotgun sequencing with Illumina and Roche 454

For whole genome shotgun sequencing, we used the Illumina GAIIx [102] to generate a high sequencing coverage (40×) of the genome based on short reads (2 × 100 bp or 2 × 115 bp) and a lower coverage of the genome (5×) based on long reads (400 bp or 750 bp) from Roche/454 pyro-sequencing technology [103]. All sequencing libraries for whole genome shotgun sequencing were constructed from about 200 μg genomic DNA of a single female domesticated canary. Illumina paired end libraries with 200 to 300 bp insert size as well as the 454 single read libraries were prepared manually according to the manufacturers’ protocols (Illumina TruSeq DNA sample preparation kit, v2, and Roche rapid library prep kit, v2). The Illumina 300 to 500 bp paired-end library was prepared semi-automatically by a Beckman-Coulter SPRI-T robot. Illumina paired end libraries were sequenced on an Illumina GAIIx sequencer.

To reconstruct large genomic scaffolds from the assembled contigs, mate-pair libraries with 2 to 3 kb insert size were prepared manually using Illumina’s biotinylation/circulization protocol (Illumina mate pair v2 library preparation protocol) and sequenced on an Illumina Hiseq2000 sequencer, while Roche/454 20 kb mate pair libs were constructed at Eurofins/MWG using the Cre-recombinase circularization approach from Roche (Roche Diagnostics GmbH, Mannheim, Germany). All 454 libraries were sequenced on GS FLX(+) sequencers.

M3. Processing of raw sequencing data

While Roche/454 read data can be used for de novo assembly directly, Illumina data need to be filtered to improve de novo assembly results. Furthermore, filtering Illumina reads is strongly recommended to reduce the need for computing resources (CPU time, RAM and disk space) during genome assembly. During read processing we clipped the reads from 3′ and 5′ ends by choosing the longest region of the read that had no base with quality lower than 11 and trimmed adaptor sequence, if at least 15 bp of the adaptor were found. We used only reads longer than 63 bp (as 64 bp is the minimum read length for assembly with the Celera Assembler; see below). Duplicate fragments generated due to a PCR bias are misleading for genome assemblers and cause false coverage estimates; hence, they were removed and only reads whose first 64 bp had a unique sequence were considered for analysis. Scripts for filtering reads were implemented in Perl and Awk. The average error rate of the trimmed sequence was below 1% according to statistics from Newbler Assembler.

M4. Genome assembly and scaffolding

We performed different whole genome shotgun assemblies during the different phases of the project using SOAPdenovo (v.1.05) [104], Newbler (v.2.6; Roche) or Celera Assembler (CA6.1) [105]. Among those three assemblers CA6.1 gave the best results in terms of N50 contig and scaffold statistics, while taking the most computing time. The reported canary genome was assembled using CA6.1. The screening of the 2 to 3 kb mate pair data by mapping them to scaffolds of prior assemblies to remove unmated/chimeric read pairs and duplicate read pairs was especially helpful to improve running time and assembly quality of CA6.1.

The genome assembly quality is reflected in the contig N50 length and scaffold N50 length of the canary genome assembly (Table 1). A few potential mis-assemblies in CA6.1 scaffolds became evident as inter-chromosomal rearrangements when mapping them against the zebra finch genome. Scaffolds were split at these locations and subsequently improved by the SSPACE scaffolder [106] with 2 to 3 kb and 20 kb mate pair reads. In a final SSPACE scaffolding step, de novo assembled scaffolds were ordered into larger scaffolds, now called superscaffolds, if their ends could be mapped closely (applied windows: 40 ± 20 kb, 200 ± 100 kb, 500 ± 250 kb) and with correct orientation in zebra finch chromosomes. This collinearity (see Materials and methods section M6) to the zebra finch chromosomes resulted in the canary having the qualitatively best bird genome among those produced using second-generation technologies alone (Figures 2 and 3; Table 2; Additional file 2).

M5. Genome alignment with other bird genomes

The LAST tool [107] was used to align available bird genomes from Ensembl or UCSC databases (budgerigar (Melopsittacus undulatus), chicken (Gallus gallus), collared flycatcher (Ficedula albicollis), medium ground finch (Geospiza fortis), turkey (Meleagris gallopavo), zebra finch (Taeniopygia guttata)) as well as genomes from the CLIMB database (Adelie penguin (Pygoscelis adeliae), emperor penguin (Aptenodytes forsteri), rock dove (Columba livia)) and NCBI GENOME database (mallard (Anas platyrhynchos), scarlet macaw (Ara macao), Tibetan ground tit (Pseudopodoces humilis), white-throated sparrow (Zonotrichia albicollis)) to the canary and zebra finch genomes (Figure 3). We found that the speed of whole genome alignments significantly increased when setting the parameter -m from 10 (default) to 1, while maintaining high sensitivity needed to align canary versus chicken and turkey genomes. The raw alignments were clustered according to their reference and query coordinates to remove suboptimal alignments (single_cov2 from TBA/MULTIZ package [108]) and further processed by the BlockDisplaySatsuma script from the Satsuma v.1.17package [109]. Resulting blocks of collinearity were converted to bed format for visualization in the canary genome browser. A complete overview of alignments against the zebra finch chromosomal sequences was done using CIRCOS [46].

M6. Collinearity improves the genome assembly

While sequence similarity on the nucleotide level goes down to the 70 to 80% range when aligning distantly related bird genomes, collinearity is much more conserved (Figures 2 and 3). Most species have 50% of their genomes in collinear regions larger than 1 Mbp, even when draft assemblies are aligned. To demonstrate that collinearity drastically improves the genome assembly, we used the chicken and zebra finch genome, removing the assembled parts that were based on physical maps or bacterial artificial chromosomes (BACs). We downloaded the zebra finch sequences from the UCSC genome browser. The chromosomal sequences were split into their original scaffolds using information from the assembly track of the browser, and we removed chrUN sequences (it seems that chrUN harbours many sequences redundant with the chromosomal sequences of the assembly). The resulting scaffolds comprised 1,928 scaffolds with a total length of 1,057,275,226 bp. N50 scaffold length was 12.3 Mbp (which is comparable to the canary assembly without collinearity scaffolding). These scaffolds were aligned against the G. gallus gal4 assembly (UCSC) using the LAST aligner. The alignments were filtered to represent collinear regions of both genomes. We then used several homemade scripts to rebuild the chicken chromosomes with the zebra finch sequence (meaning zebra finch nucleotide sequence ordered according to the chicken genome). The rebuilt genome was split into random paired-end reads of 40 kb (10 × 10⁶ pairs), 200 kb (2 × 10⁶ pairs) and 500 kb (0.8 × 10⁶ pairs) insert size by applying the ‘simulate_reads’ tool from the clc assembly simulator cell v.3.1.0. These reads were used for scaffolding of the original zebra finch scaffolds with the SSPACE v2.0 program (allowing a ±50% change from the original insert size). This resulted in superscaffolds with an N50 length of 36,752,594 bp (a nearly three-fold increase). We compared our collinearity scaffolding results with the zebra finch reference assembly. We found about 93 to 95% (depending on including or excluding unordered pieces of the zebra finch chromosomes in the analysis) of the scaffolds were correctly linked by our approach. The collinearity scaffolding approach was also able to assign 100 zebra finch sequences comprising about 19 Mbp of sequence to scaffolds that were assigned to a chromosome, but not located in the chromosomes of the reference sequence (called chrXY_random in the reference assembly) [41].

M7. The karyotype of the canary

The karyotype was produced as follows: cell cultures of canaries were performed as described [110] with modifications, using dissociated cells following incubation in collagenase for 1 h. Chromosomes were obtained by standard arrest with colcemid (Gibco, Invitrogen, Carlsbad, CA, USA), hypotonic treatment with 0.075 M KCl, and cell fixation in methanol/acetic acid (3:1). Diploid number definition and karyotype ordering were performed in conventionally stained metaphases (Giemsa 5% in 0.07 M phosphate buffer, pH 6.8).

M8. Genome annotation

To perform homology/transcript-based annotation we downloaded protein datasets from Uniprot (aves, human), GenBank (aves) and Ensembl (zebra finch, chicken, and turkey) and aligned them with SPALN2.04 [111] to the genome, to the de novo assembled transcriptome and to the cufflinks/cuffmerge assembled transcriptome.

Resulting coding sequence (CDS) models from alignments to the genome were clustered and the best scoring model was kept for each cluster. Alignment of proteins against the transcriptome data was used to assign UTR/CDS coordinates to the transcripts’ genomic coordinates and remove putative residual introns as well as remove transcripts that contained stop codons or frame shifts. Assembled transcripts as well as CDS predictions were applied to train the AUGUSTUS gene prediction tool [112], which was used together with hints from transcripts or protein alignments to predict ab initio gene models.

To construct a reference gene set, we clustered transcript-based gene models and protein alignment models by their CDS coordinates and selected the best scoring model (the transcript CDS scores were weighted 20% higher to give preference to those models). Finally, AUGUSTUS gene models whose CDS coordinates were non-redundant with transcript-based or homology-based CDS annotations were added to the final gene set.

Automated annotation of the 18,818 reference genes with BLAST2GO [113] resulted in 62% (11,671) BLAST2GO annotated, 9.4% (1,769) annotated but not meeting BLAST2GO criteria for valid annotation, and 16.2% (3,044) having BLAST hits but no GO assignment; 12.4% (2,334) of the genes had no hits to any database. Thus, combining strong and weak assignments, 87.6% (16,484) of the reference genes have an assigned function. To be stringent we based the analysis of the transcriptional activities of the HVC, RA and ENT on a functional (manual) annotation, which identified 10,573 expressed genes for the HVC, 10,572 for the RA, and 11,427 for the ENT. Our annotations increased the number of functionally annotated songbird genes by approximately 40% (4,900 newly annotated genes) compared with previous work [41].

M9. Repeat analysis and CpG island prediction

To annotate and mask repetitive sequences in the genome, we applied the RepeatModeler package for de novo repeat finding [114]. The resulting canary repeat library was combined with known aves repeats from the RepeatMasker libraries and used to mask and annotate repeat elements in the canary genome. CpG island detection was done by CpGplot from the EMBOSS package.

M10. Transcriptome assembly and differential expression analysis with RNA-seq

We deeply sequenced the transcriptomes of the HVC (n = 3 biological replicates), RA (n = 2 biological replicates), and ENT (n = 3 biological replicates) of canaries and of the HVC (n = 3 biological replicates) of zebra finches using RNA-seq (Additional file 1) and obtained >170 million paired-end reads on average for each sample. All these animals were reproductively active adult males as judged by the size of their testes and their blood testosterone levels (Materials and methods section M1). After quality filtering and adapter removal, we mapped approximately 150 million reads per sample to each reference genome with TopHat using sample-wise insert size parameters, a common splice junction library, and otherwise default parameters. The remaining reads either represented repetitive genomic regions or did not pass the quality threshold. Next, we applied reference-based transcriptome assembly with cufflinks using default parameters [115]. In addition, to include genes and transcripts with very low expression that might otherwise be missed from our transcriptome assemblies, we mapped heterologous protein sequences from public annotations of human, zebra finch and chicken against the canary genome with SPALN [111]. We merged the RNA-seq-derived tissue transcriptomes with cuffmerge and clustered them together with homology predicted transcripts to create organism-specific integrated meta-transcriptomes. The reconstructed transcriptomes are visualized and accessible as combined and tissue-specific tracks in the genome browser [57].

To identify relationships of homologous genes between species, we established two methods. First, we aligned all species meta-transcriptomes against each heterologous species genome and identified co-clustering of potential homologous transcripts within genomic regions with cufflinks. This method turned out to be very robust and, in many cases, the homology could be identified on the level of transcript isoforms. Second, in cases where no homology relationship could be assigned with genomic co-clustering, we used reciprocal blast [116] to identify the likely homologous gene.

To quantify gene expression, we first calculated raw read counts with HTseq and estimated normalized read counts (baseMean) for each of the 12 samples on the basis of the respective organism-specific meta-transcriptomes with DESeq. Differentially expressed genes for each pairwise comparison were identified using the negative-binomial model of read counts as implemented in the DESeq Biocondoctor package. The DESeq R/Bioconductor package [117,118] implements a model based on the negative-binomial distribution and a local regression-based estimation for variance and mean. P-values were adjusted for multiple testing with the Benjamini-Hochberg procedure. Differential expression with adjusted P-values of <0.05 were considered statistically significant.

M11. Microarray gene expression and data analysis

Total RNA (100 ng) was processed for hybridization on the microarray using the Ambion WT Expression Kit and the Affymetrix WT Terminal Labeling and Controls Kit. The resulting cDNA was hybridized to the Custom Affymetrix Gene Chip® MPIO-ZF1s520811 Exon Array. For each experimental group and brain area, we made six biological replicates. The 5.7 million zebra finch-specific probes spotted on this array correspond to approximately 232,000 probe sets, hence 33,000 transcripts published on public databases (NCBI-ENSEMBL) and detected in studies performed at the Department of Behavioural Neurobiology, MPIO. The functional final gene annotation was optimized comparing the gene IDs with orthologous human proteins found in EL DORADO homology groups (Genomatix toolbox) and the annotation performed in this study for RNA-seq data. The hybridization was carried out over 16 hours at 45°C and 60 rpm in a GeneChip Hybridization Oven 640. The arrays were washed, stained and scanned using the Affymetrix GeneChip Fluidics Station 450 and Affymetrix GeneChip scanner 3000 7G. The CEL files were generated by the Affymetrix® GeneChip® Command Console® Software (AGCC) and quality control for evaluating the success of individual hybridizations was assessed by the Affymetrix® Expression Console™ software. Affymetrix CEL files were imported into ChipInspector software, version 21 (El Dorado Database version: E26R1204 Genomatix GmbH [119]). Differential expression was analyzed using the group-wise exhaustive analysis with false discovery rate set to zero and 10-significant probe minimum coverage [120].

Since our microarray was designed subsequent to the publication of the zebra finch genome [41], it contains a more complete representation of the zebra finch genes than the microarray used in previous studies (for example, [62]). This and the exon-based design of the presently used microrarray (see above) make comparison with the seasonality study of Stevenson et al. [62] difficult.

The microarray data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus [121] and are accessible through GEO Series accession number GSE50070 [122].

M12. Validation of area-specific gene expression by in situ hybridization

In order to validate the RNA-seq-based expression of genes that were differentially expressed in the HVC, RA and ENT of canary and zebra finch, we designed primers using Primer3web [123]. cDNAs were prepared from total mRNA isolated from a male adult canary forebrain using Superscript III First Strand Synthesis Kit (Invitrogen, Carlsbad, CA, USA) following the recommended protocol and the amplifications were carried out using Hot FIREPol DNA polymerase (Solis BioDyne, Tartu, Estonia). The amplicons were purified and inserted into the PCRII TOPO vector accurately following the supplied protocol using a TOPO® TA Cloning® Dual Promoter Kit (Invitrogen), host strain TOP10. The inserts were subjected to restriction enzyme analysis and sequencing to confirm the orientation and sequence. The in situ hybridization of sagittal cryostat sections was performed as described [20]. For autoradiography, sections were hybridized with the canary- or zebra finch-specific S³⁵-labelled RNA anti-sense probes and subsequently exposed to X-ray film (Figure 6). Films were exposed for 4 weeks. For each probe, we used sections of three adult canaries.

In total we performed in situ hybridzations for 30 selected genes. For the 60 comparisons (considering the comparison of HVC and ENT and of RA and ENT), the in situ hybridizations confirmed the results found in RNA-seq in 53 cases. Here we depict the in situs hybridizations of six genes to i) confirm differential gene expression in HVC, RA and ENT derived from the RNA-seq approach, ii) indicate a reason for mismatches, and iii) evaluate the extent of area-specific expression of such genes. The reported genes (Figure 6) are ALDH1A2 (aldehyde dehydrogenase 1 family, member A2; GenBank accession number KF571935), CNTFR (ciliary neurotrophic factor receptor (KF571934)), GLI2 (GLI family zinc finger 2 (KF571936)), GPR83 (G protein-coupled receptor 83 (KF571938)), KCNH5 (potassium voltage-gated channel, subfamily H (eag-related) (KF571937), member 5), and RASGRP1 (RAS guanyl releasing protein 1 (calcium and DAG-regulated) (KF571933)). Expression, either up- or down-regulated in HVC and RA compared with ENT, was as expected from the RNA-seq for all six depicted genes. Nevertheless, down-regulated genes in HVC or RA might still have considerable expression.

M13. Sample preparations for the creation of the canary brain protein ion library

For the creation of the canary brain protein ion library required for SWATH data processing and protein quantification, tryptic digests were prepared from a pool of RA, HVC-I, HVC-II, ENT, forebrain and cerebellum tissues as described in Shevchenko et al. [124]. We pooled the HVC and the entopallium, respectively, of three LD canaries. Briefly, 1 ml of lysis buffer (1% SDS, 10 mM Tris/HCl pH 7.4, 0.15 M NaCl, 1 mM EDTA, 1X Roche Complete Protease Inhibitor in phosphate-buffered saline) was added to the brain tissues and the tissues homogenized for 1 minute on ice using a Dounce glass homogenizer. After 1 h incubation at 4°C, lipids were removed by adding 1.4 ml tri-n-butylphosphate/acetone/methanol mixture (1:12:1) and incubating on ice for 90 minutes. Pellets were successively re-suspended in tri-n-butylphosphate (100%), acetone (100%) and methanol (100%). For trypsin digestion, pellets were re-suspended in 200 μl digestion buffer (1% n-octyl-beta-D-glucopyranoside in 50:50 acetonitrile (ACN)/8 M urea) and proteins were denatured by incubating with 5 mM dithiothreitol (DTT) for 15 minutes at 50°C and alkylated using 10 mM iodoacetamide for 15 minutes at room temperature in the dark. Then, samples were applied to 3 kDa cutoff filters (Millipore, Merck KGaA, Darmstadt, Germany) and impurities removed by repeated centrifugations with 2% ACN in 50 mM ammonium bicarbonate (ABC), 50% ACN in 50 mM ABC and 50 mM ABC. Trypsin was added to the protein samples on the filter (trypsin:protein ratio 1:40) and incubated at 37°C overnight. Tryptic peptides were eluted from the filter by centrifugation and lyophilised. In order to increase the number of peptides identified by liquid chromatography (LC)-tandem MS analysis, a first dimension of separation was achieved by high pH reverse phase chromatography. A reverse phase column (Waters, Milford, MA, USA, BEH C18, 2.1 × 150 mm, 1.7 μm) was utilised in combination with a 10 step gradient (33 minute gradient). Eleven fractions were collected per gradient and supplemented with 1× iRT standard peptide kit (Biognosys, Schlieren, Switzerland) after solvent removal and re-suspension in 5% ACN/0.1% formic acid. A second dimension of separation was achieved by using a nanoLC Reverse Phase column (Eksigent technologies, Dublin, CA, C18 ChromXP nanoLC column 75 μm id × 15 cm, ChromXP C18 3 μm 120 Å) in combination with a 200 minute gradient reaching 35% organic phase in 155 minutes and then 80% organic phase in 5 minutes. Data were acquired on an AB Sciex TripleTOF5600 mass spectrometer in a data dependent acquisition mode (AB Sciex, Framingham, MA, USA). Data processing was performed using ProteinPilot Software 4.5 (AB Sciex) with a search database containing the transcribed canary open reading frames and common contaminant peptides. Peptide identifications with a predicted false discovery rate of <5% were regarded as significant. For the generation of a canary SWATH ion library, the recorded retention time (RT) for each peptide were transformed into indexed RT (iRT) using the iRT-calculator (Biognosys, CH). Then, the eight most intense fragment ions of each peptide were extracted and combined.

M14. Protein preparation for SWATH-MS analysis and data processing

Protein samples were prepared according to Vowinckel et al. [77] from ENT and HVC by incubating the tissues at 90°C for 20 minutes in lysis buffer (0.1 M NaOH, 0.05 M EDTA, 2% SDS, 2% β-mercaptoethanol, 0.1 M acetic acid) and then homogenizing in a TissueLyser (Qiagen) operated at 18 m/s for 2 minutes with 5 mm metal beads (Qiagen). Protein precipitation was performed in 10% trichloroacetic acid and incubation at −80°C for 2 h. Protein pellets were washed twice with 80% acetone before trypsin treatment. Proteins were finally re-suspended in 0.2% Rapigest SF Surfactant (Waters) in 50 mM ABC and quantified using the BCA protein assay kit (Pierce, Rockford, IL, USA). For trypsin digestion 50 μg and 1.6 μg protein were used for ENT and HVC, respectively. Protein samples were denatured by incubating at 60°C for 30 minutes with 5 mM DTT and then alkylated with 10 mM iodoacetamide at room temperature for 30 minutes in the dark. A first aliquot of trypsin was added to the protein samples (enzyme:protein ratio 1:40) and incubated for 2 h. Finally, a second aliquot of trypsin was added to the samples (enzyme:protein ratio 1:40) and digestion performed at 37°C overnight. Surfactant removal was achieved by centrifugation after acidification of the samples with 0.5% trifluoroacetic acid. ACN was added to the samples to a final concentration of 5%.

Targeted proteomics was conducted with minor modifications as described in Gillet et al. [76]. Samples were spiked with 1× iRT standard peptide kit and subjected to the nanoLC separation as described above. Data were acquired on an AB Sciex TripleTOF5600 operating in SWATH mode, setting the SWATH m/z acquisition window to 25 Da. Data post-processing was conducted in Skyline [125]. Quantification was carried out by summing the peak areas of the three most intense SWATH transitions per peptide and considering all the peptides identified for each protein of interest. Only unique peptides were considered for quantification of the proteins of interest.

M15. Gene network analysis

We carried out different analyses using the Genomatix software suite combining several mining sources (MatInspector, Region Miner and Genomatix Pathway Systems (GePS) [126-128]). The canonical pathway and GO-term enrichment analysis were performed using the GePS-Gene Ranker option (Genomatix tool). The full normalized data derived from DESeq analysis for each species (Materials and methods section M10) was filtered by setting the measure of read abundance, that is, expression level cutoff baseMean, at 10 reads, which is a conservative criterion [129]. ClueGO, a Cytoscape plug-in [130-133], was used to interpret the differential gene expression analysis. This application fuses terms of GO as well as KEGG (Kyoto Encyclopedia of Genes and Genomes)/BioCarta pathways and creates a functionally organized GO/pathway term network that can be visualized by functionally grouped terms [63].

M16. Extracting putative promoter sequences from the genome assemblies for the detection of AREs and EREs

We created subset gtf files for assembled transcripts only present in canary or zebra finch HVC as well as a subset gtf file for HVC genes common to both species. These files were uploaded to the UCSC genome browser (canary and zebra finch). Using the table browser we extracted 1,000 bp of genomic sequence upstream of the putative transcription start sites. Thirty-five percent of the promoter sequences contained gaps in the canary, and 39% contained gaps in the zebra finch. This suggests that both genomes have similar quality regarding promoter regions. For comparisons of ARE and ERE sites between orthologous genes of the canary and zebra finch we considered only those genes without sequence gaps in the promoters of either species.

In order to determine the potential hormone sensitivity of genes, we searched for AREs (V$AREF; matrix names ARE.01, ARE.02, and ARE.03) and EREs (V$EREF; matrix names ER.01, ER.02, and ER.03) within 1,000 bp upstream of the transcription start sites. The matrix selection constraint was the presence of IR3-inverted repeated sequences (palindromes); half-palindromic sites were not considered [34]. Clearly, this approach would discard AR and ER binding sites that are far away from the transcription start sites or that differ from the canonical binding sites. However, ChIP-seq results of the ER-regulated genes in breast cancer showed that most functional ERα binding sites contain EREs [94] and that ERE sites are required for ERα-dependent cell proliferation and differentiation [95]. The putative promoter extracted sequences were analyzed using Genomatix-MatInspector tool [126].

Major datasets