Fine-tuning of seasonal timing of breeding is regulated downstream in the underlying neuro-endocrine system in a small songbird

Variation in avian seasonal timing of breeding is ultimately rooted in the underlying physiology, as, after transduction and integration of cues, reproductive timing is the outcome of a neuro-endocrine cascade along the so-called hypothalamic–pituitary–gonadal–liver axis (HPGL axis). The hypothalamus, considered to be the final integration point of environmental cues, the pituitary gland and the neural centres are widely assumed to primarily guide top-down hormonal regulation and in this way direct ovarian function to time breeding (Dawson, 2008; Tsutsui et al., 2012). Many studies have therefore focused on these upstream levels of the HPGL axis (Nakane and Yoshimura, 2014, and references therein). Although photoperiod, perceived by three types of photoreceptors (Underwood et al., 2001), is a proximate cue for birds to time breeding (Sharp, 1996; Silverin et al., 1993; Wingfield, 1993), it cannot solely explain individual year to year variation in the timing of breeding, as the change in day length over the season is invariable across years (Bradshaw and Holzapfel, 2007; Visser et al., 2004). A potential explanation for the variation in the timing of breeding is an ‘alternative, female-specific hypothesis’ where females use (changes in) supplementary cues to fine-tune downstream mechanisms at the level of the ovary and/or liver and so may regulate vitellogenesis, follicle development and the timing of egg laying (Caro et al., 2009; Lambrechts and Visser, 1999; Williams, 2012). In general, little work has integrated downstream levels in females, let alone multiple levels of the neuro-endocrine cascade in relation to cues and/or reproductive traits (Cánovas et al., 2014; Laine et al., 2019; MacManes et al., 2017; Maruska and Fernald, 2011; Maruska et al., 2011; Perfito et al., 2015).

Evidence of possible downstream mechanisms regulating the timing of breeding has been found on a few occasions. A study in two wild populations of blue tits (Cyanistes caeruleus) breeding at different times suggested that females have similar photoperiodic sensitivities but that the population differences in seasonal timing could be explained by differences in the response of the ovary to gonadotropins, or the liver to oestrogens (Caro et al., 2009). Work on great tits (Parus major) (Schaper et al., 2012a) and European blackbirds (Turdus merula) (Partecke et al., 2005) showed significant differences in egg-laying dates between females from different temperature treatments and populations, but similar plasma luteinizing hormone (LH) levels. Individual variation in luteinizing hormone receptor (LHR) transcript in the testes and developing follicles was found in dark-eyed juncos (Junco hyemalis) but, again, no differences in LH levels (Bergeon Burns et al., 2014; Needham et al., 2019). A study in male European starlings (Sturnus vulgaris) found that the inhibition of gonadal sex steroid secretion is seasonally regulated within the testes by mechanisms involving melatonin receptors and the gonadotropin-inhibiting hormone (GnIH) system present in the gonads (McGuire et al., 2011). Direct evidence for downstream regulation of the timing of breeding was, however, found in female European starlings housed with or without males (Perfito et al., 2015). Female starlings housed with males showed elevated levels of LHR, follicle-stimulating hormone receptor (FSHR) and vitellogenin (VTG) mRNA only immediately before, or coincident with, rapid yolk development (RYD), together with increased plasma yolk precursor levels (Perfito et al., 2015). This is consistent with a ‘lack of ovarian competence’ to respond to elevated circulating gonadotropins until just before egg laying. In addition, when female starlings housed without males were briefly exposed to males, mRNA levels and yolk precursor levels elevated, indicating that the ovary depends on the ‘supplemental cue’ of male presence (Perfito et al., 2015). Multiple, if not all, levels of the HPGL axis need close and simultaneous examination to enable us to identify where species differ in executing physiological mechanisms resulting in variation in their timing of breeding. This would set the stage for understanding where selection could act and how animals could respond via genetic adaptation to changing environments.

A wealth of studies measuring hormone concentrations in circulation, using endocrine and receptor agonists and antagonists to study physiological and behavioural effects, and assessment of protein levels by immunochemistry, have resulted in the extensive knowledge on HPGL axis functioning so far. However, despite this knowledge and the understanding of which cues (i.e. photoperiod, temperature, food, social cues) influence the timing of breeding, understanding of the mechanisms regulating a female's ‘decision’ to initiate egg laying is far behind. Recent and current developments in genomic technologies have started to provide new options to explore and identify the links between genetic and phenotypic variation (Cheviron et al., 2008; Fidler et al., 2007). For the great tit, a model species in ecology and evolution, such tools, including a well-annotated reference genome, have recently become available (Derks et al., 2016; Kim et al., 2018; Laine et al., 2016).

Here, we used female great tits from selection lines where birds were genomically selected for either early or late timing of breeding (Gienapp et al., 2019; Verhagen et al., 2019). Birds were subjected to two contrasting temperature environments in climate-controlled aviaries. A recent study in these great tits reported differential expression of genes in the hypothalamus under the influence of temperature, and, when females were expected to initiate egg laying, genes were highly differentially expressed in the liver but especially the ovary (Laine et al., 2019). However, because pooled samples (three females per sample) were used in that study, gene expression levels could not be related to individual egg-laying dates. Using samples from the same great tits as in Laine et al. (2019), we assessed (1) whether individual differences in gene expression levels could explain differences in individual egg-laying dates and, if so, (2) where (upstream or downstream in the HPGL axis) these differences in gene expression occur. By making use of the great tit genome (Laine et al., 2016), we took a candidate gene approach and measured individual expression levels using qPCR. Key genes known to be important mediators in reproductive endocrine pathways upstream (i.e. the hypothalamus) and downstream (i.e. the ovary and liver) in the HPGL axis in female great tits were targeted. In addition, we selected genes that are potentially important in reproductive biology from the abovementioned genome-wide study (Laine et al., 2019).

Selection lines for early and late timing of breeding in great tits were created using bi-directional genomic selection (see Gienapp et al., 2019; Verhagen et al., 2019, for details). To summarize, from wild broods of our long-term study population in the Hoge Veluwe National Park, The Netherlands, nestlings (F₁ generation) of which the mother had initiated egg laying either extremely early (‘early line’) or extremely late (‘late line’) in the wild, were brought into the aviary facilities at the Netherlands Institute of Ecology (NIOO-KNAW, Wageningen, The Netherlands), 10 days post-hatching for further hand raising. Subsequently, chicks were genotyped using a 650 SNP chip (Kim et al., 2018) to predict their genomic breeding values (GEBVs, i.e. the value estimating the relationship between genotype and phenotype based on genetic markers). Based on their GEBVs, F₁ generation individuals were selected for early- and late-line breeding pairs to produce the F₂ generation in captivity. The F₂ generation eggs were transferred to wild ‘foster nests’ for incubation and hatching. F₂­ generation chicks were also collected and hand-raised in the laboratory. In turn, the F₂ offspring were genotyped and selected to produce the F₃ generation, which was then genotyped and selected.

The selection line study results are reported elsewhere (Verhagen et al., 2019). Briefly, we found that on average early-line birds laid their eggs about 6 days earlier than late-line birds. Further, the difference in average laying date increased (from about 2 to 10 days) from the F₁ to the F₃ generation, with non-significant line effects for the F₁ and F₂ generation, but highly significant line differences for the F₃ generation (Verhagen et al., 2019).

We would like to point out here that these results were from birds housed in outdoor aviaries. For the present study, we housed the F₁ and F₂ generation birds, in their first year of age, in climate-controlled aviaries for two consecutive breeding seasons (see ‘Experimental setup’, below). As opposed to outdoor aviaries (Verhagen et al., 2019), neither selection line nor temperature environment, nor their interaction, explained the females' reproductive phenotypes (i.e. laying dates and follicle widths) in these climate-controlled aviaries (see Table S1). Those variables were thus left out of further analyses in the present study, meaning that birds originating from both generations of selection line birds and exposed to both temperature treatment were indiscriminately used to increase the sample size.

F₁ generation (n=36) and F₂ generation (n=36) selection line pairs of great tits (Gienapp et al., 2019; Verhagen et al., 2019) were housed in 36 climate-controlled aviaries (2 m×2 m×2.25 m) at the NIOO-KNAW in 2015 and 2016, respectively. Birds were subjected to an artificial photoperiod mimicking the change in natural photoperiod. In each aviary, light was provided by one full-spectrum daylight fluorescent lamp (58 W, 5500 K, True-light) and two fluorescent lamps (58 W, Philips). A roof shaft (SolaTube) provided additional natural light (total average daily light intensity ∼500 lx per aviary; Table S2). A light bulb (7 W, Philips) mimicked dawn and dusk, turning on half an hour before the lights went on and staying on half an hour after the lights went off, respectively (Caro and Visser, 2009). In addition, the pairs were subjected to two contrasting temperature environments mimicking an extreme cold spring (2013) and an extreme warm spring (2014) in The Netherlands (Fig. S1): average (mean±s.e.m.) laying dates were May 5±5.18 days (n=112) and April 11.8±5.46 days (n=124) for 2013 and 2014, respectively, in the wild long-term study population at the Hoge Veluwe. Temperature was changed every hour to follow as closely as possible the observed hourly temperatures in these years (note that the minimum temperature in the aviaries was 2°C so any temperature below 2°C in the temperature time series from outside was set to 2°C). The combination of selection line and temperature environment resulted in four groups of nine pairs: ‘early-warm’, ‘early-cold’, ‘late-warm’ and ‘late-cold’. Note, the variables selection line and temperature environment were left out of further analyses (see above), but are mentioned here to explain the experiment. Birds were fed ad libitum as reported elsewhere (Visser et al., 2011) and had water available for drinking and bathing. All F₁ and F₂ generation pairs went through two experimental breeding cycles: a ‘first breeding season’ and a ‘second breeding season’ (see below and Fig. S2). This study was performed under the approval by the Animal Experimentation Committee (DEC), Amsterdam, The Netherlands, protocol NIOO 14.10 addendum 1.

Pairs of all four groups were put in the climate-controlled aviaries at the beginning of January 2015 and 2016, under the natural photoperiod. We provided nesting material (moss and hair) from the second week of March onwards. Birds went through their breeding season during which reproductive behaviours (e.g. nest building and date of the first egg, i.e. laying date) were recorded. Laying dates were recorded as April dates (i.e. 31 March=0; 1 April=1, etc.). Birds were blood sampled bi-weekly as part of another study (Mäkinen et al., 2019). Females could choose between three nest boxes of which two were accessible to the researcher from the outside to minimize disturbance of the birds.

After this first breeding season, when birds were photorefractory and well on their way into moult (∼mid-July), days were shortened to 9 h light:15 h dark and temperature was decreased to 10°C for 7 weeks to make the birds photosensitive and temperature sensitive again (Dawson, 2015). From September onwards, birds were again subjected to the same contrasting environments as in spring, to bring them into a second breeding season within the same calendar year. Because of this, and two subsequent years (2015 and 2016) with two breeding seasons to fit in 1 year, the second breeding season (of both 2015 and 2016) started with the photoperiod and temperatures corresponding with those of 1 February instead of 1 January. As such, 1 month of photoperiodic and temperature input is missing, but it is likely that the most important period for temperature to affect the timing of breeding is from March onwards (Visser et al., 2006). SolaTubes that bring natural light from outside to the aviaries (see above) were closed, because of the mismatching photoperiods.

Females that did not initiate egg laying (n=4 in 2015, n=4 in 2016) in the first breeding season were replaced with their sisters in the second breeding season. However, the latter were not further used in this study (see ‘Statistical analysis’ section, ‘Explaining variation in mRNA expression’, below). Pairs were divided into three groups and sampled at three time points (see ‘Tissue collection and preparation’, below).

For both generations, pairs were divided into three groups (n=12 pairs per group) based on their laying dates from the first breeding season, resulting in groups with a roughly similar average laying date and distribution (Fig. S3). Three time points were chosen based on the laying dates of 2015 (F₁ generation); (1) 7 October (which corresponds to 7 March of the first breeding season) when gonadal maturation is initiated, i.e. photoperiod exceeded 11 h (Silverin et al., 1993); (2) 28 October (corresponding to 30 March) when nest building occurred in the first breeding season, but prior to laying; and (3) 18 November (corresponding to April 20) when about 25% of the females had initiated egg laying in the first breeding season. The same time points were used in 2016 (F₂ generation) to enable us to compare the experiments of 2015 and 2016, and increase sample size. For each time point, one group was sampled (both males and females, but we focus on the females in this study). Pairs of birds were caught from the aviaries, deeply anaesthetized with isoflurane (IsoFlo, Zoetis, Kalamazoo, MI, USA) and a blood sample of 300 µl was taken from the jugular vein for possible future use. Brain, ovary and liver were dissected out. Brains were flash-frozen on dry ice and stored in 5 ml RNA-free tubes at −80°C (Qiagen), whereas the other dissected tissues were placed in Eppendorf tubes and temporarily stored in liquid nitrogen. The width of the largest follicle was measured to an accuracy of 0.1 mm before freezing. All tissues were stored at −80°C until further processing. From the frozen brains, sagittal cryo-sections (40 µm) were cut (Leica CM3050 S). The hypothalamus and hippocampus were located by the use of online zebra finch brain atlases (Karten et al., 2013), such as ZEBrA (Oregon Health & Science University, Portland, OR, USA; http://www.zebrafinchatlas.org) and directly isolated from the frozen brain sections using surgical punches (Harris Uni-Core, 2.0 mm). Isolated tissue was collected into 1 ml TRIzol (Invitrogen, ThermoFisher Scientific) immediately, homogenized by vigorous vortexing and stored at −80°C until RNA isolation.

For RNA extraction from the hypothalamus, samples were defrosted and 0.2 ml chloroform added to the 1 ml TRIzol. From the liver and ovary samples, a small piece was taken, and RNA extracted using 1 ml of TRIzol. Note that for the ovary samples, we avoided using the largest follicles in order to compare between time points. RNA yield was measured on a Nanodrop 2000 (ThermoFisher Scientific) and used to adjust the concentration for cDNA synthesis.

For cDNA synthesis from the isolated RNA samples, we used the QuantiTect Reverse Transcription Kit (Qiagen). A fixed amount of total RNA (for liver and ovary samples 150 ng in 6 µl RNase-free water, for hypothalamus 50 ng RNA) was incubated in gDNA Wipeout Buffer (1 µl) for the removal of genomic DNA. cDNA was generated (final volume 10 µl) following the manufacturer's instructions (Quantitect-Qiagen). A dilution of 1:5 for hypothalamus and 1:20 for liver and ovary was used for qPCR analysis. Until analysis, all cDNA samples were stored at −20°C.

We made a list of genes (1) known to be important or potentially important mediators of reproductive biology from the literature and (2) based on RNAseq data from the same F₂ generation females used in this study (Laine et al., 2019) (Table S3). In addition, we made a list of reference genes to allow normalization of the gene expression levels (see ‘Reference genes and normalization of candidate gene expression’, below; Table S3). Primers were then built based on the great tit reference genome build 1.1 (https://www.ncbi.nlm.nih.gov/assembly/GCF_001522545.2) (Laine et al., 2016) and annotation release 101 (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/Parus_major/101/) with Geneious version 10.0.2 (Kearse et al., 2012) and tested (see ‘qPCR amplification efficiency’, below). Primers were checked against the great tit reference genome using a BLAST search to confirm that primers were specific for the intended target genes.

Amplification efficiency of each primer pair was determined through qPCR by a 5-point standard curve based on a 5-dilution series (1:10, 1:20, 1:40, 1:80 and 1:160) of cDNA samples. Most assays for the candidate genes studied showed an efficiency (E) within the desired optimal range of 90–110%. Some fell outside this range, but were nevertheless included in the analysis based on a linear relation between the inverse log­₁₀ dilution value and the cycle threshold (C_t) (R²>0.90) and a melt curve showing a single amplicon being formed. Selected primer pairs for the final candidate gene list are listed in Table S4. Relative transcript levels were measured by qPCR using the SYBR Green method (PowerUp SYBR Green Master Mix, ThermoFisher Scientific). Fluorescence was measured with the CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories) and fluorescent data were analysed with the CFX Manager Software (Bio-Rad Laboratories) from which C_t was obtained for subsequent analyses. Amplifications were always run in duplicate (in a different analysis and a different random sample order).

Although cDNA was generated from identical amounts of RNA, variation between samples may arise as a result of different reverse transcription efficiencies and RNA quality. Such variation was corrected for by normalizing the expression level of the target gene to a normalization factor (NF) based on the expression level of a set of reference genes determined for each cDNA sample (Vandesompele et al., 2002). We started out by selecting three candidate reference genes per tissue. Reference gene expression stability was calculated using the application geNorm (Vandesompele et al., 2002) based on which it was decided whether or not to include additional candidate reference genes for accurate normalization of the mRNA expression levels (see  Appendix 2). This resulted in the selection of the following reference genes: protein kinase C alpha (PRKCA), ribosomal protein L19 (RPL19) and succinate dehydrogenase complex flavoprotein subunit A (SDHA) for hypothalamus, beta-2-microglobulin (B2M), PRKCA, RPL19 and SDHA for liver, and hypoxanthine phosphoribosyltransferase 1 (HPRT), PRKCA, ribosomal protein L13 (RPL13), RPL19 and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation (YWHAZ) for ovary. Absolute amounts of cDNA were calculated by conversion of the C_t values (C×E^−Ct, with C=10¹⁰ and E=2) (Dijk et al., 2004). The absolute amounts of the candidate genes were normalized against the normalization factor (NF) calculated by taking the geometric mean from the absolute amounts of the reference genes, resulting in relative mRNA expression levels of the candidate genes (arbitrary amounts).

We used the laying dates in the first breeding season as a measure for whether females were early or late breeders in the second season. Follicle widths were log₁₀ transformed before performing simple linear regression to investigate the relationship between laying date in the first breeding season and follicle width of the largest follicle in the second breeding season. This relationship was subsequently tested for each time point.

Removing females from the data as a result of death or not initiating egg laying in the first breeding season or having unreliable mRNA level measurements resulted in n=59, n=58 and n=59 individual females for hypothalamus, ovary and liver, respectively. Individual mRNA expression data were subjected to both principal component analysis (PCA) and univariate statistical analyses, which were all performed in R (version 3.3.1). Prior to subjecting the data to PCA, we log₁₀ transformed the individual gene expression data. Using the function prcomp, PCA was performed, which consolidates the individual mRNA expression level data into new variables known as principal components (PCs), so reducing the number of dimensions of the data. These PCs allowed simultaneous assessment of expression values of the genes measured for the hypothalamus, liver and ovary and give an indication of the variables that best explain the variation in gene expression levels for each tissue. Horn's analysis was performed to determine which PCs to retain (i.e. eigenvalue>1). Assessment of the association between PCs and explanatory variables (i.e. time point, laying date and generation) was determined by performing ANOVA, with the following model: PC_x∼time point×laying date+generation. P-values were adjusted for multiple comparisons using Benjamini and Hochberg's false discovery rate (FDR) (Benjamini and Hochberg, 1995), accepting an FDR of 0.05. Females did not differ in either laying date or largest follicle width between selection lines, temperature treatments or their interaction within a generation (see above and  Appendix 1). To exclude these variables from the study, we performed an initial analysis to test whether these variables would influence individual gene expression levels. Unlike the study of Laine et al. (2019), in which genome-wide gene expression patterns were tested compared with our limited number of candidate genes, selection line, temperature environment or their interaction did not influence gene expression levels (Table S4) and were therefore left out for further analyses. Subsequently, the same procedure as applied to the PCs was used to analyse the expression level of an individual candidate gene: expression_gene∼time point×laying date (from the first breeding season)+generation.

Pearson's correlation coefficients between every gene pair possible were calculated and visualized with the rcorr and corrplot functions in R, respectively, in order to determine which gene pairs tended to change significantly (accepting a P<0.05) together, within and across the tissues examined.

There was a weak but significant negative linear relationship (r=−0.32, F_1,59=6.88, P=0.01) between laying date and width of the largest follicle (Fig. S4A). When analysing per time point (Fig. S4B), the relationship between laying date and largest follicle width went from no relationship at time point 1 (r=−0.14, F_1,17=0.33, P=0.57) to a moderate negative relationship at time point 2 (r=−0.57, F_1,18=8.81, P=0.01) and a strong negative relationship at time point 3 (r=−0.66, F_1,20=15.18, P<0.001). Also, given the significant difference in follicle widths ( Appendix 1), we are confident that the mRNA expression levels from the second breeding season are representative of the laying dates recorded ( Appendix 3).

PC1 and PC2, the dimensions with eigenvalues >1 according to Horn's analysis, explained together 86.8%, 48.1% and 73.7% of the variance in gene expression among females in the hypothalamus (n=59), ovary (n=58) and liver (n=59), respectively (Tables S5–S7).

Based on the loadings, mRNA expression of iodothyronine deiodinase type 2 (DIO2), opsin 5 (OPN5), thyrotropin releasing hormone (TRH) and nuclear factor interleukin-3-regulated protein (NFIL3) accounted for the variance in PC1, whereas mRNA expression of vasoactive intestinal peptide (VIP) in the hypothalamus explained a large part of the variance in PC2 (Table S5). In addition, the similar loadings (Table S5) and the small angle between the vectors of OPN5 and DIO2 (Fig. 1A) suggest a correlation between these genes. Females showed different candidate gene expression profiles between generations (F_2,54=143, FDR corrected P<0.0001; Table S8), as shown by two distinct but overlapping clusters along PC1 in hypothalamus (Fig. 1A). No distinction in expression profile was found when clustering females per time point for PC1 or PC2 (Fig. 1B), nor did time point explain variance in any of the PCs (Table S8); also, no association between expression of these genes and the interaction between laying date and time point was found.

In the ovary, the variance in PC1 was mainly explained by mRNA expression of the androgen receptor (AR), luteinizing hormone receptor (LHR), matrix metallopeptidase 15 (MMP15) and interferon-related developmental regulator 1 (IFRD1) (Fig. 1C; Table S6). In contrast, the variance in PC2 was mainly explained by mRNA expression of heat-shock protein family B member 1 (HSPB1), cytochrome P450 17A1 (CYP17A1) and very low-density lipoprotein receptor (VLDLR) (Fig. 1C; Table S6). Although not shown in Fig. 1, based on similar loadings in PC1 and PC2 (Table S6), expression levels of CYP17A1, ER and VLDLR were correlated. Females showed distinct differences in candidate gene expression profile between generations (F_1,54=269.57, FDR corrected P<0.0001) along PC1 (Fig. 1C; Table S9) and a gradual change in expression profiles when clustering for time point (F_2,55=22.01, FDR corrected P<0.0001) along PC2 (Fig. 1D; Table S9). PC1 and PC2, together accounting for ∼48% of the total variance, were highly significantly associated with both generation and time point (Table S9).

PC1, accounting for ∼46% of the variance among females, was associated with a laying date×time point interaction (F_{2, 53}=11.019, FDR corrected P<0.001; Table S10) and was mainly explained by mRNA expression of apovitellenin 1 (APOV1; LOC107200088), bestrophin 3 (BEST3), cathepsinE-A-like protein (CTSEAL; LOC10720510) and vitellogenin 2 (VTG2) (Table S7). Generation explained the variation in gene expression (∼28%) along PC2 (F_{2, 53}=38.09, FDR corrected P<0.0001; Table S10). Although not shown in Fig. 1, based on similar loadings in PC1 and 2 (Table S7), BEST3, CTSEAL and VTG2 were correlated in terms of expression among these females, as were mineralocorticoid receptor (MR) and HSPB1. As in the hypothalamus, but along PC2 instead of PC1, females showed overlapping but different candidate gene expression profiles between generations in the liver (Fig. 1E). Similar to the ovary, but again along opposite PCs, females showed a gradual change in expression profile over time points (Fig. 1F).

We found no differences in candidate gene expression in the hypothalamus between time points, laying dates or their interaction (Table 1). The F₁ generation females had significantly higher expression levels at each time point for DIO2, NFIL3, OPN5 and TRH compared with F₂ generation females (DIO2: F_1,57=82.52, FDR corrected P<0.0001; NFIL3: F_1,57=58.03, FDR corrected P<0.0001; OPN5: F_1,57=77.15, FDR corrected P<0.0001; TRH: F_1,57=160.51, FDR corrected P<0.0001; Fig. 2, Table 1). We found no difference in expression levels of VIP between generations (Table 1).

With the exception of FSHR, gonadotropin-inhibitory hormone receptor (GnIHR), prolaction receptor (PRLR) and steroidogenic acute regulatory protein (StAR), all candidate genes showed significant differences between generations and time points in the ovary (Table 2; Fig. S5), but only variation in mRNA expression of IFRD1 (Fig. 3A) and VLDLR (Fig. 3B) was explained by timing of breeding (IFRD1: F_1,55=6.86, FDR corrected P=0.03, VLDLR: F_1,55=13.25, FDR corrected P<0.001; Table 2).

Early-breeding females showed increased mRNA expression for both insulin-like growth factor 1 (IGF1; Fig. 4A) and VTG2 (Fig. 4B) in the liver (IGF1: F_1,55=6.53, FDR corrected P=0.03, VTG2: F_1,58=6.62, FDR corrected P=0.03; Table 3) compared with late-breeding females. Only in the liver did we find differences in mRNA expression levels explained by the laying date×time point interaction (Fig. 5, Table 3). Females showed an increase in gene expression over time points for apolipoprotein B (APOB) (F_1,53=5.10, FDR corrected P=0.03), APOV1 (F_2,53=11.58, FDR corrected P<0.0001), BEST3 (F_2,53=6.53, FDR corrected P=0.010) and CTSEAL (F_2,53=7.21, FDR corrected P=0.01), with higher expression for early-laying females compared with late-laying females at time points 2 and 3. The genes glucocorticoid receptor (GR), HSPB1 and MR only showed a generation effect (Table 3; Fig. S6).

Within and among the tissues examined, candidate genes, whether they reflect differences in timing or not, tended to change in a strong and/or significantly similar way (Fig. 6). For example, CYP17A1 expression in the ovary tended to change in a strong and similar way to expression of APOV1, CTSEAL and VTG2 in the liver. In addition, expression of HSPB1 in the ovary resembled that of APOB and APOV1 in the liver. The mRNA expression of GNIHR in the ovary showed a weak positive, but significant, correlation with VTG2 expression in the liver. Interestingly, the genes examined in the hypothalamus showed a high and significant correlation among each other, but less so when correlated to genes in the ovary and liver. Between the ovary and liver, more genes tended to change in a similar way, both positively and negatively.

Gene expression dynamics within the HPGL axis have not been well studied in avian seasonally breeding females. Using a candidate gene approach, we set out to determine whether individual differences in egg-laying dates (obtained from the first breeding season) are reflected in differences in candidate gene expression levels, and, if so, where these differences occur in the HPGL axis (upstream and/or downstream) and when these differences can be picked up towards the expected laying dates. We found significant differences in mRNA expression of candidate genes between generations in all three tissues examined. However, a correlation of candidate gene expression and egg-laying date (at the three sampling time points) was found exclusively in the ovary and liver, independent of generation. In particular, individual differences in the timing of breeding in females were significantly reflected in mRNA expression for IFRD1 and VLDLR in the ovary and IGF-1 in the liver, and earlier breeding females showed increased expression of APOB, APOV1, BEST3 and CTSEAL over time in the liver. These findings, together with other patterns found, suggest that fine-tuning of avian timing of breeding is regulated downstream in the HPGL axis. This is in concurrence with the ‘alternative, female-specific hypothesis’ (Caro et al., 2009; Williams, 2012), which awards a more prominent role for the ovary and/or liver in fine-tuning the timing of breeding (see ‘Downstream regulation of the timing of breeding’, below).

We compared gene expression levels at different time points approaching laying of the first egg, but with different individuals for each time point. The limitation here is that the same female could not be measured at each time point. There could be individual differences in responses to cues and (reproductive) physiology, which potentially decreased our power to detect patterns over time. In addition, a 3 week interval between time points is quite long, and with a wide range in laying dates, properly determining the last time point, at which most females are assumed to have initiated vitellogenesis or egg laying, posed a challenge. Further, for practical reasons indicated in the Materials and Methods, we had to leave out the January photoperiods and temperatures for the second breeding season. However, because of increased expression levels for genes involved in, for example, vitellogenesis in both this study and the genome-wide study (Laine et al., 2019), and the fact that several females had entered RYD or initiated laying ( Appendix 2), we are positive that, given the narrow time window in which this occurs, the third time point was estimated correctly. Further, we avoided using the largest follicles, which prevented inflated expression levels for (certain) candidate genes and a possible misinterpretation of the results. We used two generations of selection lines in this study, which generated significant differences in gene expression levels in the three organs examined. It is possible that the timing of the experiments and processing of the samples, for example, might have caused these differences. An alternative explanation is that year differences were the cause, but we do not have data for enough years to test this.

Interestingly, temperature treatment affected genome-wide gene expression profiles early in the breeding season (time point 1) in the hypothalamus, but not in the ovary and liver in the same samples as used here of the F₂ generation females (Laine et al., 2019). In addition, we did not find an effect of temperature treatment on gene expression levels or on the onset of egg laying or follicular growth. The latter is contrary to previous studies in great tits housed in climate-controlled aviaries, showing that the pattern of increase in ambient temperature has a direct effect on the onset of egg laying (Schaper et al., 2012a,b; Visser et al., 2009), but agrees with other studies showing that gonadal size is not affected by ambient temperature (Schaper et al., 2012a,b; Visser et al., 2011). It seems that in these females, at the beginning of the breeding season, the brain is able to perceive ambient temperature to ‘switch on’ the reproductive axis at an upstream level (perhaps in a similar way to photoperiod). However, even though temperature could possibly affect other tissues, it does not seem to directly affect gene expression in the ovary and liver to fine-tune egg laying.

The F₂ generation females showed significantly lower expression levels of DIO2, NFIL3, OPN5 and TRH at all three time points compared with the F₁ females, whereas this was not the case for VIP. These genes are involved in circadian rhythms (DIO2, NFIL3) (Cowell, 2002; Yoshimura et al., 2003), photoperiodic perception (OPN5) (Nakane et al., 2014) and regulation of the hypothalamic–pituitary–thyroid axis (TRH) (McNabb, 2007). A possible explanation for this generation difference could be that F₂ females were, on average, ∼7.5 days later in onset of egg laying. However, the F₁ and F₂ generation females followed the same photoperiod. Also, generation differences were found in the ovary and liver, but again not for all genes. We are hesitant to attribute these generation differences to different biological functioning (see ‘Limitations’, below).

The expression of IFRD1, a gene proposed to be involved in regulation of cell proliferation and differentiation (Vadivelu et al., 2004; Vietor and Huber, 2007), decreased at time point 3 compared with time point 1 for all females, as in Laine et al. (2019), but significantly for the early-laying females (Fig. 3). This is in contrast to a study in female Sprague–Dawley rats, where increased expression of IFRD1 was found in granulosa cells and cumulus oocyte complexes after administration of human chorionic gonadotropin (to mimic the LH-surge and induce ovulation), indicating potential involvement of IFRD1 in oocyte maturation (Li et al., 2016). However, this study was performed in a different time frame (hours) and on single cells compared with weeks and ovary homogenates, respectively, in our study.

The mRNA expression of VLDLR increased from time point 1 (early March) to time point 2 (late March) and decreased again at time point 3 (mid-April) in F₁ females. When taking into account that females in climate-controlled aviaries lay ∼3 weeks later (Visser et al., 2009) than wild females, this finding is consistent with expression in ovaries of European starlings (Perfito et al., 2015). However, we expected VLDLR expression to be lowest in non-breeding females (i.e. time point 1) (George et al., 1987) and highest in pre-laying females (i.e. time point 3) (Han et al., 2009).

We find clear variation in expression, though not significantly explaining variation in laying dates, between early- and late-laying females over time for genes regulating, among others, processes involved in steroidogenesis (CYP17A1, LHR) (Johnson, 2015) and follicular development through gonadotropin binding (GnIHR) (Maddineni et al., 2008). In the ovary, both CYP17A1 and LHR expression was higher for females laying early and peaked at time point 3. This increase over time, and nearing egg laying, is consistent with findings in ovary homogenates of European starlings (Perfito et al., 2015) and in follicles of the dark-eyed junco (Needham et al., 2019). These studies support the idea that LHR plays a key role in the ‘competence’ of the ovary to respond to circulating gonadotropins (see ‘Downstream regulation of the timing of breeding’, below). Further, in chicken (Gallus gallus), CYP17A1 and LHR show increased expression when follicle selection takes place and is initiated by the signalling of several receptors via cAMP (Johnson, 2015).

Earlier-breeding females showed increased mRNA expression levels over time in liver for genes involved in vitellogenesis and oocyte growth, which is consistent with differential expression levels found for these genes at time points 1 and 3 (Laine et al., 2019). APOV1 (alias apoVLDL-II) is a protein component of yolk-targeted very-low density lipoprotein (VLDLy), a lipoprotein synthesized by the liver under the influence of oestradiol (E2) and, together with VTG, the primary source of yolk protein and lipid for the developing embryo (Walzem, 1996). APOB, a protein associated with VTG and VLDLy (Walzem, 1996), and VTG2 (one of the three forms of VTG and the most abundant), showed increased expression over time compared with APOV1, BEST3 and CTSEAL. These expression patterns agree with concentrations of VTG and VLDL found in other seasonal breeders (Caro et al., 2009; Challenger et al., 2001). Like VTG and VLDL, synthesis of CTSEAL by the liver is oestrogen dependent (Zheng et al., 2018). Further, CTSEAL is allegedly involved in the sexual maturation of female chickens (Bourin et al., 2012b) and may play a role in processing egg yolk macromolecules (Bourin et al., 2012a), as it is found in egg yolk (Farinazzo et al., 2009). The function of BEST3 in this study is unclear. BEST3 is positioned close to CTSEAL in the genome, and therefore its lower expression might be caused by co-regulation with CTSEAL (Laine et al. 2019; Zheng et al., 2018).

We found expression of IGF1 to reflect individual differences in egg laying, with early-laying females showing higher IGF1 expression compared with late-laying females. There is little knowledge regarding the connection between IGF-1 and reproductive traits in birds. Few studies (mainly poultry) exist; ovaries have IGF-1 receptors and IGF-1 plays a regulatory role in ovarian functions, such as follicular growth and differentiation (Onagbesan et al., 1999), and stimulates ovarian progesterone production (Williams, 1994). Growth and reproduction are closely related and there is cross-talk between the endocrine systems controlling these fundamental processes in vertebrates (Hull and Harvey, 2014, and references therein). Studies in female chicken and rabbit suggest that IGF-1 is also produced by the ovary, together with and under the influence of growth hormone, where they act as paracrine/autocrine regulators during follicular development (Ahumada-Solórzano et al., 2016; Yoshimura et al., 1994, 1996). In addition, different variants of IGF-1 genes, as well as variation in IGF-1 levels in poultry, resulted mainly in variation in productivity, i.e. different numbers of eggs produced or variation in egg quality (Hocking et al., 1994; Nagaraja et al., 2000; Wu et al., 2016).

The limited number of candidate genes, which were not assessed in all the tissues examined, hampers the construction of a gene network and a subsequent co-expression network analysis in order to associate genes (of unknown function in relation to the timing of breeding) with biological processes. Even so, these preliminary results on correlated expression between gene pairs within and across tissues highlight the importance of looking not only within but also across tissues in the HPGL axis. Further, co-expression of these genes might indicate the same transcriptional regulatory programme (e.g. transcription factors, DNA methylation). In addition, these preliminary results emphasize the importance of the communication between the ovary and liver as a potential mechanism in the timing of breeding. For example, CYP17A1 showed significantly correlated expression with genes expressed in the liver (CTSEAL, VTG2 and APOV1; Fig. S5) that are involved in lipid metabolism and yolk formation (Walzem, 1996; Zheng et al., 2018). Of course, E2, for which CYP17A1 is a key enzyme in the steroidogenic pathway underlying its production, stimulates vitellogenesis (Mullinix et al., 1976). However, whether the ‘decision’ to lay is mechanistically linked to follicle selection and development, ovulation and ultimately egg laying remains to be investigated.

Currently, one can only speculate on where the ‘switch’ that initiates egg laying resides within the ovary and/or liver. A potential candidate is the ‘competence’ of the ovary to respond to gonadotropins via their receptors (Johnson, 2015; Caro et al., 2009; Ball, 2007; Williams, 2012; Schaper et al., 2012a; Partecke et al., 2005). Further, in starlings, it has also been shown that sex steroid secretion can be regulated by local GnIH in the gonads (McGuire et al., 2011; Kriegsfeld et al., 2015). As such, the gonadal GnIH system could be a potential mechanism in the timing of breeding in females (McGuire et al., 2011; Needham et al., 2019). Another potential mechanism is the communication between the ovary and liver, where the E2-dependent shift in lipid metabolism or the upregulation of VTG/VLDL receptors could be candidates. These potential mechanisms, however, need to be regulated, and imply a more autonomous role for the ovary together with receiving signals that bypass the classic neuroendocrine pathway. As such, the ovary and brain might act more as ‘partners’ (Ball, 2007). For example, a study in Japanese quail (Coturnix japonica) suggests that the ovary regulates its own functioning through its circadian clock, because the largest follicle, through production of circadian clock gene proteins, controls the LH surge that is essential for ovulation (Nakao et al., 2007).

The exact downstream mechanisms that precede avian timing of breeding and how they are regulated remain to be determined. Though gene expression is not the only mechanism regulating the timing of breeding, we have shown that variation in mRNA expression levels of several candidate genes in the ovary and liver, associated with reproductive functioning, explains variation in the timing of breeding in these females. Our study confirms that shifting the focus more towards females rather than males (Caro, 2012; Williams, 2012) in future experimental studies investigating the timing of breeding is highly important. Also, simultaneous examination of multiple, and preferably all, HPGL axis levels is of the essence in understanding the mechanisms underlying the timing of breeding. This way, we can gain knowledge on the variation in the physiology underlying the timing of avian breeding and what part of this variation is genetically determined. Timing of breeding is currently under selection in wild populations as a consequence of climate change (Both and Visser, 2001; Visser et al., 1998). A better understanding of the variation in the physiological processes underlying seasonal timing will ultimately lead us to a better understanding of a species' adaptive potential to their warming world.

Differences in laying date and largest follicle width between selection lines and treatments were tested by performing ANOVA, with laying date or log₁₀-transformed follicle width as the dependent variable and selection line (i.e. genomic phenotype), treatment (i.e. warm and cold temperature environment, see ‘Experimental setup’ in Materials and Methods) and generation as explanatory variables, with time point added as an extra variable in testing differences in follicle width.

Laying dates (first breeding season) ranged from April 9 to April 69 (8 June). Treatment (F_1,95=0.91, P=0.34), selection line (F_1,95=0.08, P=0.78), or their interaction (F_1,94=0.02, P=0.90) did not explain variation in laying date. We did, however, find an effect of generation (F_2,96=3.38, P=0.04): a difference in mean laying date of ∼7.5 days between the F₁ and F₂ generation (t_61.99=−2.50, P=0.02) (Fig. A1).

Also, we found no effect of selection line (F_1,55=0.27, P=0.36) or treatment (F_1,56=0.27, P=0.61), or their interaction (F_2,50=0.03, P=0.86) on the width of the largest follicle (second breeding season). Follicles were larger for the F₁ generation (F_1,58=7.24, P=0.01) and increased over time (F_2,59=32.68, P<0.0001), with the largest follicles measured at time point 3 (P<0.0001 for comparisons with both time points 2 and 3; Fig. A2).

For proper gene expression analysis, the data must be normalized against the expression level of a set of stable reference genes. This approach controls for factors such as the amount of cDNA load in a sample, variation in the efficiency of the RT reaction, and RNA quality (Vandesompele et al., 2002). Optimal reference genes exhibit stable expression levels that are not influenced by the (experimental) condition. Ideally, a set of multiple reference genes is compiled, for which the NF is then calculated (Vandesompele et al., 2002).

Absolute amounts of reference gene cDNA were calculated by converting the C_t values (C×E^−Ct, with C=10¹⁰ and E=2) (Dijk et al., 2004). Then, the measure of reference gene expression stability (M) was calculated in the application geNorm. M is defined as the average pairwise variation (V_n/n+1) between the normalization factors (NF_n and NF_n+1) of a particular reference gene with all other reference genes (Vandesompele et al., 2002).

There is a cut-off point of V_n/n+1=0.15, below which it is not necessary to include an additional reference gene for normalization (Vandesompele et al., 2002). However, we did not take this cut-off point too strictly. When all V_n/n+1 values were slightly greater than 0.15, but the NF_n and NF_n+1 showed a high correlation, we decided to stop adding reference genes. Using at least three reference genes with highly correlated expression levels is already a significant improvement on the common practice of using a single gene.

The genes PRKCA, RPL19 and SDHA were selected as potential reference genes for hypothalamus samples. All three showed M<1.5 according to geNorm, with M_RPL19 and M_PRKCA=0.632 and M_SDHA=0.691. Further analysis resulted in VPRKCA, VRPL19 and VSDHA=0.217, which is above the recommended V=0.15 (see above). Close inspection of the data revealed three individual samples having strongly deviating amounts, in both the reference gene and candidate gene dataset, indicating decreased cDNA quality in these samples. These individuals were therefore removed from the dataset and the remaining data rerun in geNorm. Stability slightly increased (M_RPL19 and M_SDHA=0.583, M_PRKCA=0.603) and V dropped closer to the cut-off point (V_{PRKCA-RPL19-SDHA}=0.179). We found very low variation (R²=0.980; Fig. A3) between NF₂ and NF₃, meaning that addition of the third reference gene did not add much to the overall normalization. Nevertheless, we used these three reference genes for normalization of mRNA expression data in the hypothalamus.

We started with HPRT, PRKCA and YWHAZ as potential reference genes for ovary samples. All three reference genes showed an M<1.5, with M_HPRT and M_YWHAZ=0.451 and M_PRKCA=1.0262. Further analysis of V in geNorm resulted in V_HPRT, V_PRKCA and V_YWHAZ=0.431, which is far above the recommended V=0.15 (see above). Therefore, we ran two extra potential reference genes, RPL19 and PRL13, in order to decrease V. Analysis in geNorm resulted in M_PRKCA and M_RPL19=0.583, M_RPL13=0.741, M_HPRT=1.010 and M_YWHAZ=1.038. Using these five reference genes enabled us to reduce V (V_HPRT, V_PRKCA, V_RPL13, V_RPL19 and V_YWHAZ=0.189). Here, V was still >0.15, but we found very low variation (R²=0.989; Fig. A3) between NF₄ and NF₅ and decided to not add a sixth reference gene for normalization of mRNA expression data in the ovary.

We started with PRKCA, RPL19 and SDHA as potential reference genes for liver samples. All three reference genes showed an M<1.5, with M_RPL19 and M_SDHA=0.606 and M_PRKCA=0.830. Further analysis of V in geNorm resulted in V_RPL19, V_SDHA and V_PRKCA=0.298, which is far above the recommended V=0.15 (see above). Therefore, we ran an extra potential reference gene, B2M, in order to decrease V. In addition, one individual was removed from the dataset because of strongly deviating amounts. Analysis of the four reference genes in geNorm resulted in all genes showing M<1.5 (M_RPL19 and M_SDHA=0.606, M_PRKCA=0.650 and M_B2M=0.677) and a decreased V (V_B2M, V_PRKCA, V_RPL19 and V_SDHA=0.149). Here, V_B2M, V_PRKCA, V_RPL19 and V_SDHA<0.15 and, together with the high R² (Fig. A3) found between NF₃ and NF₄, adding a fifth reference gene is not necessary for accurate normalization.

The largest follicle width during the second breeding season was 7.3 mm, measured on the day that this female should have laid her third egg in the first breeding season. Therefore, this width was taken as a measure for a fully developed follicle (f1 follicle; note: in this Appendix, f1, f2, etc., are used in relation to follicle order/size to avoid confusion with generation, used elsewhere). We back-calculated the approximate follicle sizes (Table A1), as we did not measure individual size differences of the f5, f4, f3 and f2 follicles, by using the traditional hierarchical model of follicle development (Astheimer and Grau, 1990). This model predicts that the first follicle to enter RYD is the first to ovulate and first to be laid. Hepatic production of VTG and VLDL, both yolk-targeted lipoproteins, is essential for vitellogenesis (i.e. yolk formation through nutrient deposition in the oocyte) and oocyte growth (Bacon et al., 1974; Walzem, 1996). We found increasing follicle width over time (Fig. S4), with three females carrying follicles similar to f3–f1 approximate sizes (Table A1) and five females likely to have entered RYD and others close. We found VTG2 mRNA expression reflecting individual differences in egg laying (F_1,58=6.625, P=0.032) and increasing over time (F_2,57=56, P<0.0001). In addition, we found a significant relationship between laying date and follicle width, especially at time point 3 (Fig. S5).We are therefore confident that the mRNA expression levels from the second breeding season are representative of the phenotypes (i.e. laying dates) recorded and assume these breeding seasons to be similar.

We thank Marylou Aaldering, Coretta Jongeling, Franca Kropman and Anouk de Plaa for taking care of the birds. We also thank Davide Dominoni and Barbara Helm for help in selecting candidate genes, Renske Jongen for assistance during tissue collection, and Jeroen Laurens and Gilles Wijlhuizen for technical assistance prior to and during the experiments. We are grateful to the Netherlands Institute of Neuroscience, Amsterdam, The Netherlands, for making a laboratory available for cryo-sectioning the brain material.