Developmental stress does not induce genome‐wide DNA methylation changes in wild great tit (Parus major) nestlings

The environment experienced during early life is a crucial factor in the life of many organisms. This early life environment has been shown to have profound effects on morphology, physiology and fitness. However, the molecular mechanisms that mediate these effects are largely unknown, even though they are essential for our understanding of the processes that induce phenotypic variation in natural populations. DNA methylation is an epigenetic mechanism that has been suggested to explain such environmentally induced phenotypic changes early in life. To investigate whether DNA methylation changes are associated with experimentally induced early developmental effects, we cross‐fostered great tit (Parus major) nestlings and manipulated their brood sizes in a natural study population. We assessed experimental brood size effects on pre‐fledging biometry and behaviour. We linked this to genome‐wide DNA methylation levels of CpG sites in erythrocyte DNA, using 122 individuals and an improved epiGBS2 laboratory protocol. Brood enlargement caused developmental stress and negatively affected nestling condition, predominantly during the second half of the breeding season, when conditions are harsher. Brood enlargement, however, affected nestling DNA methylation in only one CpG site and only if the hatch date was taken into account. In conclusion, this study shows that nutritional stress in enlarged broods does not associate with direct effects on genome‐wide DNA methylation. Future studies should assess whether genome‐wide DNA methylation variation may arise later in life as a consequence of phenotypic changes during early development.

mechanisms that mediate these effects are largely unknown (Sepers et al., 2019). Yet, knowledge on these mechanisms is essential for understanding the role of environmental variation on phenotypic variation and associated fitness consequences in natural populations (Beldade et al., 2011;Bossdorf et al., 2008).
Early ontogenetic effects have been studied extensively in mammals, including humans, and altricial birds (Lindstrom, 1999;Lummaa & Clutton-Brock, 2002;Wells, 2011), since their offspring rely on parental care for a significant part of their development.
One of the conditions that offspring experience during their rearing period that is known to have profound phenotypic effects is the number of siblings the early development needs to be shared with (Spencer & Tilbrook, 2009). The effect of the number of siblings (i.e. brood size) has been extensively studied in altricial birds as brood size varies naturally (Perrins & McCleery, 1989) and is relatively easy to manipulate experimentally. Brood size is known to affect stress levels during the nestling period (Naguib et al., 2011), due to for example increased parasite loads, high thermoregulation costs and stronger sibling competition (Carere et al., 2005;Fargallo & Merino, 2004;Mertens, 1969). Indeed, experimental brood enlargement is associated with increased testosterone (Naguib et al., 2004) and corticosterone levels (Saino et al., 2003), leading to several negative pre-fledging effects, such as comprised immunocompetence (Naguib et al., 2004;Saino et al., 2003), accelerated telomere loss (Boonekamp et al., 2014), impaired growth, shorter wing length (Martyka et al., 2021;Nettle et al., 2013), delayed fledging (Nettle et al., 2013) and ultimately higher nestling mortality (De Kogel, 1997). Apart from morphological consequences, experimental brood enlargement and the accompanying increased sibling competition might also affect nestling behaviour, as the handling stress response has been found to negatively associate with the received amount of caterpillar biomass in wild great tit (Parus major) nestlings (van Oers et al., 2015).
Epigenetic mechanisms have repeatedly been suggested to mediate early environmental effects on the phenotype. Epigenetic mechanisms are biochemical mechanisms that alter gene expression without a change in the nucleotide sequence of the genome (Richards, 2006). The most-studied epigenetic mechanism is DNA methylation (Korochkin, 2006), which is the addition of a methyl group (-CH 3 ) to a DNA nucleotide, usually cytosine (C). DNA methylation interferes with the binding of proteins required for transcription initiation and usually represses gene expression, especially if the promoter is methylated (Weaver et al., 2004). Epigenetic changes can be caused by genetic variation and can be induced as a response to the environment (Richards, 2006). Methylation variation, therefore, encompasses a promising mechanism explaining early developmental effects in birds (Sepers et al., 2019).
Genome-wide DNA methylation changes occur during the nestling period in wild great tits at least up until Day 15 after hatching (Watson et al., 2019), suggesting the existence of a window in which DNA methylation can mediate the environmental role in shaping the phenotype. Early environmental effects on DNA methylation in blood have been found in three previous studies of experimental brood size (Jimeno et al., 2019;Sepers et al., 2021;Sheldon et al., 2018). In a previous study (Sepers et al., 2021), we manipulated brood size and found differentially methylated CpG sites (CpGs) in genes involved in development, growth, metabolism, behaviour and cognition, when comparing pools of 14-day-old nestlings reared in reduced or enlarged broods. This suggested that the effects of early rearing conditions on physiology and behaviour are mediated by changes in DNA methylation. Although this study and others underline the potential role for DNA methylation in the regulation of developmental phenotypic plasticity in altricial birds, these studies focused on genome-wide DNA methylation levels, single genes, anonymous loci and most importantly, were based on small sample sizes or even pools of individuals. In our former study, we used pooled samples to assess treatment effects, but we also found differences between pools of control broods, indicating the likelihood of picking up false positives is large in such underpowered studies. Furthermore, the previous studies could only speculate on the relationship between DNA methylation and phenotypic traits, as none of them linked individual variation in DNA methylation to early life effects on phenotypic traits. Thus, even though studies have now predicted that environmental effects on phenotypic traits may be mediated by epigenetic mechanisms (Ledón-Rettig et al., 2013;Romano et al., 2017), we have limited insight in the influence of the early environment on links between epigenetic mechanisms, functional genes and phenotype traits in natural populations.
To assess whether experimentally induced phenotypic changes are induced by DNA methylation changes during early development, we cross-fostered great tit (P. major) nestlings and manipulated the brood size in which they were reared, in a natural study population. This enabled us to assess brood size effects on biometric measures, developmental (handling) stress and DNA methylation. We used epiGBS2 to measure DNA methylation of individuals at the CpG level . EpiGBS2 is a cost-effective and timeefficient reduced representation DNA methylation analysis tool, based on a bisulfite sequencing version of Genotype by Sequencing (GBS), allowing for larger sample sizes. The previous version of this tool, epiGBS (van Gurp et al., 2016), lacked a user-friendly pipeline and has mainly been applied to plant species and has hardly been used in vertebrates (McNew et al., 2021). However, the performance of epiGBS2 has been successfully validated using both plant and great tit samples .
We expected that an experimental increase in brood size induced nutritional stress and impaired nestling conditions (De Kogel, 1997;Dubiec et al., 2006;Martyka et al., 2021;Naguib et al., 2004;Nettle et al., 2013;Saino et al., 2003;Smith et al., 1989;Tinbergen & Boerlijst, 1990) and increased the handling stress response (van Oers et al., 2015), caused by decreased parental feeding rates per nestling (Gow & Wiebe, 2014;Mathot et al., 2017). If the effects of brood size dependent-stress on biometry and behaviour are mediated by DNA methylation, we hypothesise to detect differential methylation in functional regions of genes involved in metabolism, morphology and growth that are proportional to the phenotypic effects (Jimeno et al., 2019;Sepers et al., 2021). Furthermore, since we expect a functional relationship between the handling stress response and gene expression in the brain and because significant correlations between DNA methylation in blood and brain have been found in the great tit (Derks et al., 2016), we hypothesise to detect differential methylation in functional regions of genes involved in neuron development and brain development as well (Sepers et al., 2021).

| Subjects, study site and general procedures
This study was conducted from April to June 2018 in a long-term nest box population of great tits (P. major) in the Westerheide estate near Arnhem, the Netherlands (52°01′00N, 05°50′30E).
Westerheide is a mixed wood forest with 228 nest boxes distributed over a 1000 × 1200 m area. From the first week of April onwards, we checked nest boxes twice a week to determine the initiation of breeding activity, date of first egg-laying, clutch size and start of incubation, for calculating hatch dates. We visited broods 2 days before the expected hatch date to determine the exact hatch date.

| Cross-fostering and brood size manipulation
A total of 32 broods (with 226 nestlings) were assigned to a crossfoster pair (hereafter: CF pair, N = 16) based on similarity in hatch date and brood size. Four broods in two CF pairs had one day difference in hatching date and were therefore excluded from the analyses. Therefore, the numbers below refer to broods with the same hatch date within a CF pair.
We conducted a partial cross-foster design once the nestlings were one to two days old. Out of the 14 CF pairs, five were crossfostered on Day one after hatching (N = 10 broods, 75 nestlings) and nine were cross-fostered on Day two after hatching (N = 18 broods, 151 nestlings). The method of van Oers et al. (2015) was used for cross-fostering. A schematic overview of the cross-fostering procedure and the brood size manipulation is provided in Figure S1 in Sepers et al. (2021). In brief, on the Day of cross-fostering, nestlings within broods were ranked based on their weight (which was measured using a digital scale, ±0.01 g) and were then assigned to a brood. One brood within the CF pair received all the nestlings with even ranks and the other brood received all nestlings with odd ranks.
With this method, weight differences between nestlings that stayed in the brood of origin and cross-fostered nestlings were minimized (van Oers et al., 2015). Half of the nestlings were unmoved, while the other half were cross-fostered. In addition to the cross-fostering procedure, the original brood sizes were manipulated. One brood within a CF pair received three nestlings more than the original brood size (+3 nestlings, N = 14 enlarged broods, 153 nestlings), while the other brood received three nestlings fewer than the original brood size (−3 nestlings, N = 14 reduced broods, 73 nestlings). This approach was favoured over having a specific number of nestlings in reduced and enlarged broods (e.g. 5 and 11), as a breeding pair's original brood size is its optimal brood size in terms of reproductive success (Tinbergen et al., 1990). Therefore, an enlarged brood size of for example 11 nestlings could be manageable for one breeding pair, because it is close to the original brood size, while it is extremely challenging for a different breeding pair. Despite the brood size manipulation, all broods contained unmoved and cross-fostered nestlings and we aimed for minimal weight differences between unmoved and crossfostered siblings. A brood within a CF pair was randomly assigned to become either an enlarged brood or a reduced brood, unless reducing the number of nestlings would result in three or fewer nestlings in a brood (N = 4 CF pairs). In those cases, the smaller brood was assigned to become an enlarged brood, while the other brood in the CF pair was assigned to become a reduced brood. In all other cases, the manipulation was independent of the original brood size. Before the manipulation, the brood size did not significantly differ between the treatment groups (treatment, LMM: F 1,16 = 1.12, p = .31, mean brood size ± SE = enlarged: 7.81 ± 0.36; reduced: 8.25 ± 0.36). Directly after the manipulation, enlarged broods were significantly larger than reduced broods (treatment, LMM: F 1,16 = 182.25, p = 3.67 × 10 −10 , mean brood size ± SE = enlarged: 10.88 ± 0.37; reduced: 5.19 ± 0.37).
Right before cross-fostering and weighing, the down tufts on the head, wings and back of the nestlings were selectively plucked to be able to identify individuals and their brood of origin (van Oers et al., 2015) up until Day six after hatching. Six days after hatching, the nestlings were weighed and ringed with uniquely numbered aluminium bands (Vogeltrekstation). Since some eggs hatched after implementation of the partial cross-foster design and some nestlings were missing or found dead or complete broods were deserted, we weighed 232 nestlings from 27 broods on Day six after hatching. Of these nestlings, 158 were raised in 14 enlarged brood and 74 in 13 reduced broods (14 complete CF pairs and one incomplete CF pair).

| Feeding rate
On Days nine to ten after hatching, parents were caught using spring traps placed inside the nest box. They were ringed, weighed, aged (second calendar year or older) and the length of their third primary (P3) was measured (ruler, ±0.5 mm). Subsequently, adults were fitted with a passive integrated transponder (PIT) tag. This is a radiofrequency identification (RFID) tag embedded in a polypropylene leg ring (Eccel Technology Ltd, Glenfield, UK). On Day 11 or 12 after hatching, the adult feeding rate was recorded for one full day, making use of RFID readers (Dorset ID), with antennas mounted around the nest box opening. When a parent entered or exited the nest box, its individual PIT tag number was recorded. From the RFID data, the total feeding rate per day per adult was determined. For females, records within 38 s of each other were treated as one feeding visit, and for males, records within 36 s of each other were treated as one feeding visit. These figures were based on comparing video recordings of actual feeding visits with PIT tag recordings. Records within these minimum time intervals were likely caused by an adult sitting near the antennae or exciting the nest box. For each female, the first recording was ignored as females roost in the nest box and the first recording therefore reflects her exiting the nest box. Subsequently, the total feeding rate was divided by the total number of nestlings in the brood to obtain parental feeding rates per nestling per day.
One brood that was measured on Day 13 after hatching and another brood with only one provisioning parent were removed from the analysis. As a result, the provisioning behaviour of 16 parents from nine enlarged broods and 25 parents from 13 reduced broods (nine complete and four incomplete CF pairs) were included in the analysis.

| Handling stress test
The nestlings were subjected to a handling stress test as described in Fucikova et al. (2009) on Day 14 after hatching. The nestlings were taken from the nest box and placed in a bird bag. The nest box' entrance hole was blocked to prevent the parents from seeing an empty nest. Individual nestlings were taken from the bird bag one by one and held on their back in the palm of one hand while counting the number of breast movements (i.e. breath rate) during four subsequent bouts of each 15 s (period 1). After recording the breath rate, each nestling was socially isolated from the other nestlings by placing it in a wooden box with individual compartments for about 15 min. Subsequently, the breath rate was recorded again as described on the previous page and in the same order (period 2).
Since breath rate is affected by ambient temperature, the air temperature at the test location was also measured. Once all nestlings were tested, they were weighed, their tarsus length was measured (calliper, ±0.1 mm), and they were blood sampled from the brachial vein (approximately 10 μL). Each sample was stored in 1 mL cell lysis buffer (Gentra Puregene Kit, Qiagen) at room temperature until analysis. Once all nestlings were measured and bled, they were placed back in the nest box. Broods within a CF pair were always measured and bled on the same day. On Day 14 after hatching, 191 nestlings from 25 broods were measured. Of these nestlings, 118 nestlings were from 12 enlarged broods and 73 nestlings from 13 reduced broods (11 complete CF pairs and three incomplete CF pairs).
We were able to record the breath rate in all nestlings, expect for one nestling from an enlarged brood. Dead nestlings were removed from the nest box but not replaced by other nestlings. Despite this, a significant difference in brood size was still present on Day 14 after hatching (treatment, LMM: F 1,9.52 = 103.22, p = 2.03 × 10 −6 , mean brood size ± SE = enlarged: 10.1 ± 0.51; reduced: 5.30 ± 0.50).
To obtain individual estimates for the handling stress response, a similar method as described in Fucikova et al. (2009) andvan Oers et al. (2015) was used. In short, we ran a linear mixed model (LMM) with ML estimation, with the number of breaths per 15 second bout as dependent variable. The interaction between individual and bout nested in period was included as a fixed effect and individual was included as a random effect. Time, temperature and temperature squared were included as covariates. As we did not want to correct for the effect of the experiment, we did not include weight, tarsus length and date. The model produced two estimates for every individual, which refer to the individual deviation from an average slope in breath rates for each period. The estimates of the second period (i.e. handling stress) were extracted from the model for further analysis.

| Variance partitioning
To determine how much of the variance in weight on the Day of cross-fostering, Day six and Day 14, tarsus length and handling stress could be attributed to either the prehatching or the rearing environment, five LMMs with REML estimation were fitted. All dependent variables were standardized (rescaled into a mean of zero and a standard deviation of 1). In each model, CF pair was included as a fixed factor to account for variation between CF pairs of broods due to for example seasonal effects. Brood of origin and brood of rearing, both nested within CF pair, were included as random effects. Furthermore, for the model on weight on the Day of crossfostering, the age of the nestlings at cross-fostering (1 or 2 days after hatching) was included as fixed effect to correct for age effects on weight. The significance of the random effects was determined with a Likelihood Ratio Test (LRT) using the ranova function in R.

| Treatment effects
To analyse the effect of the treatment (experimental brood size) on handling stress, weight and tarsus length on Day 14, separate LMMs with ML estimation were used. Treatment, hatch date and hatch date squared (second-order raw polynomial term) were included as fixed effects. Since environmental conditions such as food availability change over the breeding season (Naef-Daenzer & Keller, 1999;Naef-Daenzer et al., 2000;Wilkin et al., 2009), the interactions between treatment and hatch date and hatch date squared were also included to test whether the treatment effect was consistent over the breeding season. CF pair was not included as a fixed effect to prevent correlation between CF pair and hatch date as differences between CF pairs might be caused by seasonal effects. Brood of origin and brood of rearing both nested within CF pair were included as random effects.
To analyse the effect of the treatment on parental feeding rates, two LMMs with ML estimation were fitted. In the first model, the parental feeding rate per day was the dependent variable and in the second model, the parental feeding rate per nestling per day was the dependent variable. In both models, treatment, age of the nestlings (Day 11 or 12), hatch date and hatch date squared were included as fixed effects. The interactions between treatment and hatch date and hatch date squared were also included to test whether the treatment effect was consistent over the breeding season. Brood ID nested within CF pair was included as a random effect.
To be better able to compare the results of the current study to the ones of our previous study performed in 2016 (Sepers et al., 2021), we performed several analyses to check for any differences between the two years. To check if the original brood size (before cross-fostering and brood size manipulation), rearing brood size on Day 14 after hatching (after cross-fostering and brood size manipulation) and parental traits (tarsus length, weight and P3 length) differed between the years, we fitted a linear model (LM) for each independent variable. Year and sex were included as fixed factors, although the latter was only included in the models involving parental traits to correct for parental sex.

| DNA methylation analysis: Sample selection and processing
For each brood from a complete CF pair, we randomly selected the Day 14 samples of two to three cross-fostered and two to three unmoved nestlings (i.e. 2 × 2 or 2 × 3 siblings). This resulted in a total sample size for DNA methylation sequencing of 67 samples originating from 11 enlarged broods and 55 samples from 11 reduced broods, together belonging to 11 complete CF pairs. 2.6.1 | epiGBS2 library preparation and sequencing Genome-wide DNA methylation levels were assessed using epiGBS2 , which is a reduced-representation DNA methylation laboratory protocol and analysis tool. Library preparations were done at the NIOO-KNAW and as described in Gawehns et al. (2022) with some improvements. These improvements entail increasing the amount of DNA per sample from 400 to 800 ng and conducting size selection before adapter ligation (see laboratory protocol in the Supporting Information). In short, from each sample, 2 to 5 μL blood was used. DNA was isolated using the FavorPrep™ 96-Well Genomic DNA Extraction Kit (Favorgen).
From each sample, 800 ng DNA was extracted and checked on a NanoDrop 2000 spectrophotometer (Thermo Scientific). Next, the DNA was digested with the restriction enzyme MspI (NEB), which cuts the genomic DNA at C^CGG motif sites. Subsequently, large fragments were removed using beads (0.8× AMPure XP beads; Beckman Coulter). The fragments were ligated to customized barcoded adapters. A unique adapter combination was used for each sample within a library. We varied the used BA adapters within a library as much as possible to increase the complexity of the first read, as the first couple of bases of the first read are used the calibrate the PCR signal. Next, the fragments of all samples were pooled (i.e. multiplexed). Small fragments (<60 bp, e.g. adapter(s) (dimers)) were removed with the NucleoSpin Gel & PCR Cleanup Kit (Macherey-Nagel). Subsequently, the nick between adapter 5′ ends and fragment 3′ ends was repaired and the DNA fragments were treated with the chemical sodium bisulfite. This chemical does not affect methylated cytosines, but converts unmethylated cytosines (Cs) into uracils, which are amplified as thymines

| Bioinformatics and statistical analysis DNA methylation
We analysed the epiGBS2 output according to the recommendations outlined in Laine et al. (2022). 2.7.1 | Demultiplexing, quality control and trimming Raw reads were demultiplexed using epiGBS2, checked for quality and adapter content and merged as described for the great tit samples in Gawehns et al. (2022). To trim, the 3′ adapter sequence cutadapt v2.10 (Martin, 2011) was used. First, the standard illumina sequence (specifying AGATCGGAAGAGC) and short reads (< 20 bp) were removed. From illumina sequence trimmed reads, additional 10 bp were removed from the 3′ end, to also discard the custom part of the adapter sequence. Quality improvement of the reads was verified by FastQC v0.11.8 (Andrews, 2010), FastQ screen v0.11.1 (Wingett & Andrews, 2018) in bisulfite mode and MultiQC v1.8 (Ewels et al., 2016 (Laine et al., 2016) and to align the trimmed reads. The reads were aligned in non-directional and paired-end mode (average mappability: 46.61%, range: 39.30%-50.60%).
Methylation calling was also done using Bismark in paired-end mode and overlap between read pairs was removed using the -no_overlap option to prevent double scoring of overlapping methylation calls.
Calling was only done in CpG context as methylation hardly occurs outside CpG context in great tit erythrocyte DNA (Derks et al., 2016;Laine et al., 2016). The results were summarized using MultiQC, and careful observation of the m-bias plot revealed a "fluttering" pattern in DNA methylation in the first 14 bps in both the R1's and the R2's.
This could indicate a bias in DNA methylation due to the introduced cytosine during end-repair (Bock, 2012), sequencing errors or basecalling errors (Hansen et al., 2012;Taub et al., 2010). To prevent a bias in DNA methylation due to read position, methylation was called again while ignoring the first 4 bps in the R1's and R2's.

| Filtering of methylation calls
Filtering of the data was done with Rstudio v1.4.1717. The R package methylKit v1.16.1 (Akalin et al., 2012) was used to read the datafiles into Rstudio and to to merge complementary CpG dinucleotides.
MethylKit was also used to exclude sites that were covered by fewer than 25 individuals in each treatment, as this removed a great portion of redundant sites and to keep enough free memory for the subsequent filtering steps. Next, sites with low coverage (<10×) were removed. To avoid a possible PCR bias in the statistical tests (i.e. higher coverage), percentile filtering (99.9%) was applied for every individual. Again, sites that were not present in enough individuals (minimally 25) in each treatment were removed to prevent significant differences between treatments due to differences in sample size. Subsequently, per CpG site, the level of methylation was calculated by dividing the number of methylated Cs by the coverage of that site. Sites with a mean methylation level across all samples of lower than 5% or higher than 95% were discarded. This way, many sites with low variation in methylation level across all individuals were excluded. After filtering, 268,951 out of the 3,983,861 CpGs were used for further analysis (Table S1).

| Statistical analysis DNA methylation
To analyse the effect of the treatment on DNA methylation level per CpG site, generalized linear mixed models with a binomial error structure and a logit link function (negative binomial) were used.
We accounted for variation in coverage (Lea et al., 2017;Zhang et al., 2018) by modelling the dependent variable as the fraction of the number of Cs (number of methylated Cs) over the number of Cs plus the number of Ts (unmethylated Cs) using the cbind function.
To assess differential methylation between nestlings of different experimental groups, two GLMMs for every CpG site were run. In the first model, treatment and CF pair were included as fixed effects.
To test whether treatment dependent methylation differences occurred as a function of hatch date, we ran a second model in which we included treatment and hatch date and the interaction between the two as fixed effects. Hatch date was standardized (rescaled into a mean of zero and a standard deviation of 1). We did not test for a quadratic relationship between treatment and hatch date as we did not have methylation data for broods that hatched at the beginning or the end of the breeding season.
The two models described above contained brood of origin and brood of rearing as random effects. For every CpG site, the significance of treatment or the interaction of treatment with hatch date was determined with a likelihood ratio test (LRT) by comparing a model with the factor of interest with a model without the factor of interest using the R packages lme4 v1.1.27.1 and lmerTest v3.1.3 and the anova function in R. p-Values were corrected for multiple testing using a false discovery rate (FDR) threshold of 0.05 (Benjamini & Hochberg, 1995), based on 268,951 analysed CpGs. To correct for potential overdispersion, we calculated the dispersion statistic for each GLMM (Zuur et al., 2013). The 95% Highest Density Interval (HDI) for the distribution of the dispersion statistic was calculated using the R package HDInterval v0.2.2. Sites that fell out of the 95% HDI were excluded. An overview of the statistical filtering steps and the number of CpGs is provided in Table S2.
A small number of nestlings were sampled on Day 13 (N = 5 nestlings) or 15 (N = 5 nestlings) after hatching instead of on Day 14 (N = 112 nestlings). We did not exclude those nestlings, since they were present in both treatments (mean age ± SE = enlarged: 14 ± 3.68 × 10 −2 ; reduced: 14 ± 3.67 × 10 −2 ). Including all 122 nestlings created more variation within the treatments group if there would be an effect of age, making this approach conservative. However, we repeated the analyses described above while correcting for age of sampling, which did not change the results. Furthermore, some nestlings were cross-fostered, while others stayed in the brood where they were born. Although this significantly affected DNA methylation levels in 211 CpGs (Table S2), we did not exclude those sites, since cross-fostered and unmoved nestlings were present in all broods and in both treatments (in the epiGBS2 data set: 34 cross-fostered and 33 unmoved from enlarged broods and 25 crossfostered and 30 unmoved from reduced broods). We therefore continued with the original GLMMs (without age of sampling and cross-fostered/unmoved).
In the case of significant interactions of treatment with hatch date, post hoc comparisons and estimates of slopes were performed with the emtrends function. p-Values were corrected for multiple testing using an FDR of 0.05 (Benjamini & Hochberg, 1995).

| Gene annotation
The CpGs were annotated with the P. major reference genome

| Variance partitioning
The variance in the weight of nestlings right before cross-fostering was not explained by the brood of origin or by the brood of rearing ( effect increased to 68% of the variance, while brood of origin explained only 8% of the variance (Table 1).
Variance in tarsus length on Day 14 after hatching was explained significantly by both genetic and rearing effects and 21% and 13% of the variance could be attributed to brood of origin and brood of rearing, respectively (Table 1). Variance in handling stress response on Day 14 after hatching was also significantly explained by origin and rearing effects as 32% and 12% of the variance could be attributed to brood of origin and brood of rearing, respectively (Table 1). where nestlings reared in enlarged broods were lighter (Figure 1a).

| Treatment effects on biometry and behaviour
There was a lack of such an effect for nestlings that hatched during the first half of the breeding season.
Similarly, a non-linear pattern was visible for tarsus length, where smaller tarsi were found for nestlings from enlarged compared to reduced broods, but only for late hatching broods (treatment × hatch date 2 , LMM: F 1,169.45 = 11.81, p < .001; Figure 1b;

| Treatment effects on parental feeding rate
Possibly the parents of enlarged broods were able to compensate for the sudden increase in the number of nestlings during the early part of the breeding season, but not later on. The parental feeding rate per day was significantly higher in the enlarged compared to the reduced broods (treatment, LMM: F 1,41 = 17.02, p < .001; Figure 3a;  Figure 3b; Table S7) and this effect was also independent of the hatch date of the nestlings. These results indicate that parents provisioning enlarged broods were not fully compensating for the increased number of nestlings in their brood.

| Treatment effects on DNA methylation
The effect of the treatment on DNA methylation changed linearly over hatch date in one CpG (treatment × hatch date 2 , GLMM: In all other CpGs, we did not find any significant differential methylation when comparing DNA methylation between nestlings from experimentally enlarged and reduced broods (treatment, GLMM: all FDR corrected p-values ≥ .05) or when assessing a hatch date-dependent treatment effect (treatment × hatch date 2 , GLMM: all FDR corrected p-values ≥ .05; Table S2).

| DISCUSS ION
Studies have suggested DNA methylation to be one of the key mechanisms underlying developmental phenotypic plasticity (Bentz et al., 2016;Sepers et al., 2021;Sheldon et al., 2018). Here, we experimentally tested this by assessing phenotypic and DNA methylation effects of experimental manipulation of brood size in a partial cross-foster experiment in wild great tits. Although brood enlargement negatively affected nestling condition in the second half of the F I G U R E 2 Handling stress response Day 14. Relationship between hatch date (April day) and the handling stress response on Day 14 after hatching for both treatments. Raw data points, regression lines and 95% confidence intervals are plotted for nestlings from the reduced treatment and the enlarged treatment .
F I G U R E 3 Parental feeding rate. Parental feeding rate (a) per day and (b) per nestling per day for both treatments. Raw data points (parents) and estimated mean ± SE (indicated by error bars) are plotted for nestlings from the reduced treatment and the enlarged treatment .

| Origin and rearing effects of early development
The minimal contribution of the common origin environment and the major contribution of the common rearing environment (Brinkhof et al., 1999) we found on variance in nestling weight and an increase of the rearing effect have been demonstrated before (van Oers et al., 2015). Our results show that the cross-foster experiment successfully separated brood of origin effects from rearing effects. The variation in handling stress response and tarsus length was associated with the origin environment and hardly with the rearing environment, indicating heritable variation and not supporting the large scope for early developmental effects in contrast to earlier findings (van Oers et al., 2015). Most of the variance could not be explained by either brood of origin or brood of rearing (i.e. residual), which points towards considerable within-brood variation.

| Provisioning and nestling biometry
Brood enlargement negatively affected nestling condition, as reflected by a lower body mass and shorter tarsus length in nestlings from enlarged broods. However, this effect only appeared in the second half of the breeding season. These results are in agreement with previous studies that have found negative effects of experimental brood size on growth and development (De Kogel, 1997;Dubiec et al., 2006;Naguib et al., 2004;Nettle et al., 2013;Saino et al., 2003;Smith et al., 1989;Tinbergen & Boerlijst, 1990). The effects on condition indicate that nestlings in enlarged broods suffered from nutritional stress, which could have several causes. Despite an increase in the total feeding frequency, parental feeding rates per nestling were lower in enlarged broods, suggesting that parents did not compensate fully for the increased food demand in enlarged broods, which is line with earlier research (Baldan et al., 2019;Gow & Wiebe, 2014;Hinde & Kilner, 2007;Mathot et al., 2017;Wegmann et al., 2015).
The phenotypic effects of our experiment were present only during the second half of the breeding season. Several possible reasons may cause this season-dependent treatment effect. It is likely that the parents were able to (partly) compensate for an increased food demand by delivering bigger or more nutritional prey only early in the season. Due to a decrease of such profitable prey later in the season or a decrease in food availability in general, this compensation was probably not possible later in the season. Indeed, previous studies have found seasonal effects on caterpillar density and mass (Naef-Daenzer & Keller, 1999;Naef-Daenzer et al., 2000;Wilkin et al., 2009).
Another cause of these date-dependent effects could be that late-breeding adults were of lower quality (Verhulst et al., 1995) and therefore less able to sufficiently compensate for increased demand.
Although the relationship between laying date and parental quality is still unclear, laying date and several physiological variables are linked (Norte et al., 2010), indicating that females in good physiological condition were able to breed early. As male quality has been related to the proportion of spiders fed (Pagani-Núñez & Senar, 2014), it might be that the low quality parents did not deliver enough of the right prey species in enlarged broods. Therefore, effects of brood enlargement might have only appeared in broods of low-quality parents. It could also be that the late breeding parents occupied lower quality habitats. For example, the proportion of caterpillars in the nestling diet depends on the distance to the nearest oak tree (Wilkin et al., 2009) and this might hinder prey selectivity, which is already reduced in enlarged broods (García-Navas & Sanz, 2010;Mathot et al., 2017;Wright et al., 1998). However, it is very hard to disentangle seasonal variation from environmental variation and variation in parental quality throughout the season, and they are probably not mutually exclusive.

| Handling stress response
We also found a date-dependent effect of experimental brood size on the handling stress response, Contrary to our expectations, the handling stress response decreased over the breeding season, especially in nestlings from enlarged broods. As our results strongly suggest that the nestlings from enlarged broods suffered from nutritional stress, especially in the second half of the breeding season, and as the handling stress response is negatively associated with the amount of caterpillar biomass (van Oers et al., 2015), F I G U R E 4 DNA methylation level of the significant CpG located downstream of FOXC1 and GMDS. The relationship between hatch date and CpG methylation level for both treatments is displayed for the CpG that showed a significant hatch date-dependent treatment effect. Raw data points, regression lines and 95% confidence intervals are plotted for samples from the reduced treatment and the enlarged treatment .
we expected a stronger stress response in nestlings from enlarged broods. It is possible that the nestlings from the reduced broods showed a stronger stress response due to for example higher thermoregulatory costs (Neuenschwander et al., 2003). Indeed, captive great tit nestlings from small broods showed an increased stress response compared to control broods, although this was not the case in wild great tit nestlings (Naguib et al., 2011). Our somewhat contradictory result as compared to Naguib et al. (2011) andvan Oers et al. (2015) might be explained by the experimental induction of nutritional stress in our study. Nutritional stress might affect physiological pathways related to stress, which could explain a different physiological response towards handling compared to previous studies. Overall, the treatment effect on the stress response was likely mediated via differences in nestling condition (Fucikova et al., 2009).

| DNA methylation
Despite observed effects on nestling condition, we found an effect of experimental brood size on DNA methylation in only one CpG and only if we took the timing of the season into account. This CpG was located downstream of the genes FOXC1 and GMDS. In nestlings from the enlarged broods, methylation of this site decreased with hatch date. High FOXC1 expression is associated with a low percentage of abdominal fat in chickens (Sun et al., 2013) and dynamic expression of GMDS has been found in during the laying period in chicken breast muscle (Zhang et al., 2022). Therefore, methylation of this CpG might regulate or reflect nestling condition, especially since nestling weight also decreased with hatch date in the enlarged broods in our study. However, as this CpG did not occur in the regulatory region of either FOXC1 or GMDS, there is no evidence for a potential change in gene expression (Laine et al., 2016) and related phenotypic traits.
We did not find an effect of experimental brood size on DNA methylation in any other CpG. Thus, although developmental processes are affected by nutritional stress, this does not seem to be mediated by DNA methylation. Effects of brood size on DNA methylation in blood have been found in three previous studies (Jimeno et al., 2019;Sepers et al., 2021;Sheldon et al., 2018). In wild zebra finches, nestlings from naturally larger broods showed higher levels of DNA methylation, although this was not found in experimentally changed brood sizes (Sheldon et al., 2018), indicating that DNA methylation differences are likely caused by quality differences between the parents.
Also in a previous study, we manipulated brood size and found 32 CpGs that were differentially methylated when comparing pools of great tit nestlings reared in reduced and enlarged broods (Sepers et al., 2021). There were some differences between the previous and the current study. First, the cross-fostering took place at a slightly earlier age in the current study. This could have affected our results, as this would have exposed the nestlings to the treatment longer than the nestlings in our previous study. One might expect this to induce a stronger effect on DNA methylation as compared to our previous study. As we observed the opposite, we do not think this explains the (more or less) lack of a treatment effect in the current study. Second, in the year that our former study was conducted, the brood sizes before and after brood size manipulation were significantly smaller than in the current study.
As great tits seem to adjust their clutch size in response to predation and breeding pair density (Julliard et al., 1997;Perrins & McCleery, 1989), Hence, methylation levels in a pool have a high chance of being biased towards a single individual, which might lead to an increased chance of false positives. Since we also found 17 differentially methylated sites between pools from control broods in our former study, we assume that due to the low sample sizes we picked up several false positives. Therefore, we consider the current results from our improved experimental design as much more robust. This clearly shows the importance of a large enough sample size when conducting differential methylation analyses with bisulfite methods (Laine et al., 2022).
CpGs in the putative promoter region of the glucocorticoid receptor gene (NR3C1) were differentially methylated between captive zebra finches reared in large broods and small broods (Jimeno et al., 2019). Although methylation of several CpGs in or near this gene were analysed in our study as well, methylation was not significantly affected by the treatment. Where the study on captive zebra finches focused on one candidate gene using pyrosequencing, our studies focused on genome-wide CpG methylation. The number of tests that are done simultaneously in our study, made our multiple testing burden quite intense, but at the same time diminishes the risk of picking up false positives. In combination with the lower accuracy of genome-wide methods, we may have been unable to pick up slight difference in DNA methylation such as with more precise methods such as pyrosequencing, where lower differences can be detected (Laine et al., 2022).

| Future direction and limitations
The lack of a genome-wide treatment-effect during the nestling period suggests that DNA methylation variation is largely genetically programmed during early development (Dubin et al., 2015). In this study, we assessed CpG methylation in red blood cells.
It is under discussion how well blood DNA methylation levels reflect DNA methylation in other tissues (Husby, 2020). Significant correlations between DNA methylation in blood and brain (Derks et al., 2016) as well as correlations between changes in DNA methylation in blood and other tissues  have been found in the great tit. However, the direction and strength of the relationship between DNA methylation levels in blood and other tissues might strongly differ between CpGs and between individuals . Thus, even though we found an effect of experimental brood size on DNA methylation in only one CpG, we cannot exclude the possibility that other sites were affected in other tissues. Also since we did not measure gene expression in the tissues we expect the molecular processes related to phenotypic changes to take place, and the relationship between DNA methylation and gene expression is tissue-specific . Further research, therefore, has to validate these findings.
In contrast to blood, non-CpG methylation occurs in vertebrate brain cells (de Mendoza et al., 2021;Pinney, 2014), including great tit brain cells (Derks et al., 2016;Laine et al., 2016). As we expect a functional relationship between the brain and one of the traits assessed in this study, the handling stress response, we cannot rule out the possibility that non-CpG methylation contributes to the observed variation in this behavioural trait. Although we recognize the limitations that come with assessing DNA methylation in blood, it would have been necessary to sacrifice wild animals to obtain inaccessible tissues such as brain tissue. Sacrificing individuals would have ruled out the possibility to follow them throughout their lives and to study long-term or carry-over effects of brood size on DNA methylation.
Furthermore, assessing blood methylation levels allowed for comparisons with other epigenetics studies (Jimeno et al., 2019;Sepers et al., 2021;Sheldon et al., 2018). In conclusion, given the largely tissue-generality in DNA methylation (Derks et al., 2016;, we expect the number of CpGs affected by experimental brood size in blood to be representative of the number of affected CpGs in other tissues, but functional validation is needed to verify this.

| CON CLUS ION
This study shows that nutritional stress has consequences for nestling body mass, tarsus length and the handling stress response, but not for genome-wide DNA methylation in blood during the nestling period. Therefore, DNA methylation is unlikely to reflect the current physiological condition during early development. To conclude, this study shows the limit for the direct role of early environmental effects on DNA methylation and emphasizes the need for studies on how differences in early developmental affect the post-fledging phenotype and its plasticity, possibly via changes in patterns of DNA methylation.

AUTH O R CO NTR I B UTI O N S
Kees van Oers and Bernice Sepers conceived the study. Bernice

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest. j1bv.

E TH I C S S TATEM ENT
All animal experiments involved in this study were reviewed and approved by the Institutional Animal Care and Use Committee (NIOO-IvD) and were licenced by the CCD (Central Authority for