Differential effects of steroid hormones on levels of genetic variance in a wild bird: possible mechanism of maternal-effects x genetic variance interaction?

Genetic variation is one of the key concepts in evolutionary biology and an important prerequisite of evolutionary change. Still, we know very little about processes that modulate its levels in wild populations. In particular – we still are to understand why genetic variances often depend on environmental conditions. One of possible environment-sensitive modulators of observed levels of genetic variance are maternal effects. In this study we attempt to experimentally test the hypothesis that maternally-transmitted agents (e.g. hormones) may influence the expression of genetic variance in quantitative traits in the offspring. We manipulated the levels of steroid hormones (testosterone and corticosterone) in eggs laid by blue tits in a wild population. Our experimental setup allowed for full crossing of genetic and rearing effects with the experimental manipulation. We observed, that birds treated with corticosterone exhibited a significant decrease in genetic variance of tarsus length. We also observed less pronounced, marginally significant effects of hormonal administration on the patterns of genetic correlations between traits expressed under varying pre-hatching hormonal conditions. Our study indicates, that maternally transmitted substances such as hormones may have measurable impact on the levels of genetic variance – and hence, on the evolutionary potential of quantitative traits.


Introduction
Evolutionary change relies on the existence of genetic variance in phenotypic traits (Fisher, 1930;Lande and Arnold, 1983;Lande and Shannon, 1996). According to the general theorem of selection, evolutionary change in a phenotypic trait is equal to the genetic covariance between the trait and fitness (Price, 1970;Lynch and Walsh, 1998;Teplitsky et al., 2014;Walsh and Lynch, 2018). It implies that genetic variance in a trait is a prerequisite of its evolutionary potential. Nevertheless, most of the available evidence for the role of genetic variance in traits' evolution comes from laboratory populations and planned breeding systems (including agricultural and artificial selection for specific, desirable properties of organisms) (Drobniak and Cichoń, 2016;Walsh and Lynch, 2018) -which may bias genetic parameters by exposing organisms to conditions unlike those experienced in nature. Far less estimates come from natural populations -even though recently one can observe a new wave of wild-population heritability estimates, in most cases owing to long-term study systems reaching data capacity required for genetic parameters estimation (Postma, 2014). Moreover, we not only still have limited understanding of evolutionary forces operating in wild populations -but also likely underestimate the extent of variation in evolutionary parameters (such as selection gradients and genetic variances), and factors driving this variation. Variation in quantitative genetic parameters is one of the factors proposed to contribute to maintaining genetic variance in the wild (Houle, 1992;Hoffmann and Merilä, 1999;Teplitsky et al., 2014;Gienapp and Brommer, 2015), hence studying it is one of the most important avenues in evolutionary research.
Available evidence from wild populations and long-term projects following individuals across multiple generations suggests, that levels of genetic variance observed in nature vary greatly with respect to a number factors (Merilä and Fry, 1998;Oltman et al., 2005;Brommer et al., 2008;Galloway et al., 2009;Gunay et al., 2011;Schroeder et al., 2012). Firstly -different classes of traits are characterised by varying levels of genetic variance, with life-history and fitness-related traits conventionally demonstrating lower levels of genetic variance, compared to morphological or physiological traits (Postma, 2014). This pattern is usually attributed to a much closer link of lifehistory traits to actual fitness, and thus stronger directional selection acting on such traits -as compared to morphological or physiological traits (Walsh and Lynch, 2018). However, heritabilities of traits can also greatly vary between contrasting environmental conditions, both in terms of external environments and "internal" environments constituted by specific traits and physiological states of organisms. For example, multiple studies demonstrated that stressful environments may induce decreases in the observed levels of genetic variance, sometimes even to the point of completely removing estimable genetic variance in them (Hoffmann and Merilä, 1999;Teplitsky et al., 2014). Sexes also can differ in heritabilities of sex-specific variants of certain traits, often resulting in lower trait heritabilities in males compared to females (Jensen et al., 2003;Foerster et al., 2007;Poissant et al., 2010;Hallsson and Björklund, 2012;Wyman and Rowe, 2014). This pattern was in turn invoked when explaining existing evidence showing that males are often more sensitive to environmental conditions than females (but see Jones et al., 2009). Finally, recent evidence suggests, that local environment generated by parents can also influence observed heritabilities. Such "environment" may not always reflect typical notions of external habitats and can for example encompass sets of conditions linked to certain characteristics of parents (e.g. parental age, a characteristic demonstrated to affect genetic variance in the offspring (Kim et al., 2011;Drobniak et al., 2015)).
Mechanisms underlying the abovementioned modifications of genetic variances are largely unknown. Irrespectively of the underlying factor influencing genetic variance, several hypotheses were brought up to explain observed patterns of heritability variation. It is possible that external/local/internal environments experienced by individuals modulate expression levels of genes at the very basic, molecular level, resulting in the observed patterns in quantitative genetic parameters (Jensen et al., 2003;Fox and Wolf, 2006). Mutation accumulation, specific to certain environments or individual characteristics, could also be responsible for such patterns (such explanation was so far considered mainly in the context of senescence and genetic variance increase with age) (Wilson et al., 2007;Charmantier et al., 2015). Surprisingly, no study was performed so far to experimentally probe the possible mechanisms behind condition-dependent expression of genetic variance.
In our study we decided to directly follow one of the possible mechanisms. A large body of evidence suggests that on the quantitative genetic level additive genetic effects can interact with genetic maternal effects (Wilson et al., 2005;Galloway et al., 2009;Hadfield, 2012;McAdam et al., 2014). Estimates of a significant covariance between maternal and individual genetic effects were so far discussed mostly in the context of how they influence evolutionary trajectories of phenotypic traits (Wilson et al., 2005). However, maternal effects have also been extensively studied in the context of their proximate mechanisms. Females can influence phenotypes of their offspring via an array of processes, some of which involve maternal transfer of resources and biologically active compounds to the offspring at the stage of eggs or developing embryos (Groothuis and Schwabl, 2008;Wolf and Wade, 2009;Coslovsky et al., 2012). Some of those mechanisms can directly affect the expression of genes in the offspring. Steroid hormones are a good candidate: they have profound effects on offspring development (Hayward and Wingfield, 2004;Tobler and Sandell, 2007;Tschirren et al., 2009;Coslovsky et al., 2012;Ruuskanen et al., 2012;Schweitzer et al., 2013;Lutyk et al., 2017); they directly impact the expression of genes by acting as transcription modulators in the nuclei of cells after binding to their specific receptors (Kawata, 1995;Baker, 1997;Podmokła et al., 2018); they are often involved in modulation and transduction of external cues such as stressful external conditions (Coslovsky et al., 2012;Podmokła et al., 2018;Öst et al., 2019). As such, steroid hormones are well documented as mediators of maternal effects, and one reason for this is the ease of manipulating their levels in eggs of wild birds (Tschirren, Saladin, et al., 2005;Ruuskanen et al., 2009).
Here we employed a direct manipulation of levels of two steroid hormones (testosterone and corticosterone) in eggs of the blue tit (Cyanistes caeruleus) -a model species in evolutionary ecology. Our experiment involved cross-fostering at the egg stage and a fully crossed, factorial design ensuring that genetic and environmental sources of trait variation in these birds would be fully interacted with the hormone level manipulation. The choice of hormones also was not accidental: testosterone and corticosterone are well studied, both in a laboratory and wild population contexts . They differ in physiological impact and mediate different kinds of information (corticosterone being a well-established mediator of stress responses, whereas testosterone being involved in primary sexual characters development and reproductive investment regulation (Groothuis and Schwabl, 2008)). We predicted, that the levels of genetic variation in certain phenotypic traits would be affected by this manipulation: assuming that by supplementing hormones we would simulate a stressful or male-like sex-specific reaction, we expected a decrease in genetic variance in hormone treated birds, compared to control birds receiving a sham manipulation (Hoffmann and Merilä, 1999;Jensen et al., 2003).

General field methods
The experiment was performed in a wild population of blue tits, studied since 2002 on the Baltic island of Gotland, Sweden (57°01' N; 18°16' E) in three breeding seasons (2014)(2015)(2016). In this population blue tits breed in wooden nest-boxes distributed uniformly across 23 study plots of varying size; density of breeding pairs is uniform across plots of different size (unpublished data). Most plots are covered by oak (Quercus robur), ash (Fraxinus excelsior) and poplar (Populus sp.) forests, with dense common hazel undergrowth (Corylus avellana). Some plots lack the undergrowth and are covered by bright, loose oak forests with wet, rich meadows abundant in orchids. In the studied population, tits lay almost exclusively one clutch per year. Females lay on average 11 eggs (range: 5 -17) and incubate them for 13 days; chicks fledge at the age of 17 -20 days.
All breeding attempts were regularly inspected by visiting all available nest boxes every 4-5 days, recording nest construction/egg laying stage, and determining species occupying each nest-box (except for blue tits, the population is also home to great tits (Parus major) and collared flycatchers (Ficedula albicollis)). Selected nests were assigned to experimental triplets (see Fig. 1 and the next section for a more detailed description of the experiment). Figure 1 provides a summary of procedures and measurements performed in each nest. Parents in each nest were ringed with aluminium bands (if not having one already), measured for tarsus length, wing & tail length and body weight. Age of adults was determined based on the presence of moult limits in the tail and between primary and secondary wing coverts (Demongin, 2016). Sex was determined by the presence of a brood patch in females. Experimental nestlings were also injected (in two of the three seasons) with a small dose of phytohaemagglutinin (PHA, Sigma-Aldrich, Germany) to determine their cell-mediated immune responsivity (Sarv and Horak, 2009;Demas et al., 2011). Briefly, 0.2 mg of PHA suspended in 400 ul of buffered saline was injected into the right wing-web of each nestling on the 11 th day post-hatching. The thickness of the web prior to the injection, and 24 h afterwards was determined with three measurements using a pressure-sensitive spessimeter (Mitutoyo, Japan model 7313) to the nearest 0.01 mm. The difference between averaged triplets of "after" and "before" measurements quantifies the amount of swelling resulting from PHA hypersensitivity reaction and is treated as a proxy of cellmediated immune response. Wing web thickness is positively correlated with body size; to account for this, all assayed nestlings were also weighed on the 12 th day. The PHA assay was performed only in years 2015 and 2016. Nonetheless, our analyses still are robust and valid: PHA treatment (or lack thereof) was always applied to all three hormonal groups (i.e. it could not generate observed differences between hormone-treated groups) and all nests not treated with PHA are grouped in one year (i.e. possible effects of not receiving PHA injections on other measured variables, however small, are linked to the year effect and fully explained by the year factor).
Blood samples retrieved from nestlings were used to determine the sex of each chick, using a well-established protocol described (Griffiths et al., 1998). Briefly, after isolating bird DNA a PCR was used to amplify a fragment of the chromohelicase (CHD) gene located on sex chromosomes and exhibiting a sex-specific length dimorphism, scored after separating the PCR products on an agarose gel. In some nestlings (Table S3) the sex could not be assigned due to technical reasons (not enough genetic material for reliable PCR, failure of the PCR reaction or ambiguous result with the gel bands markedly differing in intensity). Field procedures conformed with the legal requirements of Sweden (permit from Jordbruksverket to LG; Swedish ringing licence RC712 to SMD).

Steroid-injection experiment
Nest inspected during egg-laying were grouped into triplets based on their equal laying dates. In each triplet, at the stage of 9 laid eggs, a hormone injection manipulation was performed. The eggs were taken out of their nests and safely transported to the field laboratory. For the time of manipulation, females were left with an equal number of dummy plastic eggs. After transporting to the lab, the eggs were weighed, photographed and labelled. Each egg was assigned by random to one of three experimental groups: testosterone group (T4), corticosterone group (CORT), and control group (C). Each egg was individually marked with a non-toxic marker and then injected with 3 ul of experimental solution. In group C this was pure sesame oil. In group T4 each dose contained 3.8 ng of testosterone (17β-hydroxy-4-androsten-3-on; Sigma-Aldrich, Germany); in the CORT group each dose contained 1.2 ng of corticosterone (11β,21-dihydroxyprogesterone; Sigma-Aldrich, Germany). Testosterone was dissolved directly in the sesame oil, whereas corticosterone (due to its poorer solubility in oil) was first dissolved in absolute ethyl alcohol (Gam et al., 2011), and then 10 μl of such stock dissolved in the sesame oil. To make the groups fully comparable, oil in the C and T4 groups was also spiked with 10 μl of 100% ethanol. Nonetheless, resulting solutions quickly evaporate the residual ethyl alcohol, which anyway would be present in a concentration of 0.5% and less.
The doses of hormones were determined following a hormone assays on randomly chosen unmanipulated eggs from the studied population, sampled in preceding seasons (T4 concentration ± SD: 2.13±0.81 ng/yolk; CORT concentration: 0.61±0.26 ng/yolk; N = 10). These values are close to published estimates from the blue tit and the closely related great tit (Tschirren et al., 2004;Vedder et al., 2007;Kingma et al., 2009;Lessells et al., 2016). Final doses were calculated as the mean concentration plus 2 SDs rounded up to the nearest 0.1 ng, therefore ensuring that hormone concentrations in manipulated eggs would exceed natural levels by at least 2 SD of their natural values.
Injections were performed using a 25 ul 702RN Hamilton micro-syringe with type-4 26s removable needle. Each experimental group had its own syringe, we also used a number of replacement needles kept in 96% ethanol in order to keep them clean and sterile. Prior to each injection the egg was gently swabbed with a small amount of ethanol to disinfect a portion of its shell. Then, a disposable sterile needle was used to make a small hole in the shell, and the Hamilton syringe was inserted through it, under the visual control thanks to a flashlight illuminating the egg from beneath it. In order to make sure that the content of the syringe was injected into the yolk, we performed several trial injections on eggs from deserted nests, using a food dye as the injected liquid. After freezing, these eggs were cut open to verify that the injection procedure delivered the liquid into the yolk and only there. After injection the hole in egg's shell was closed with a drop of Vetbond (3M, Minnesota, USA), a tissue adhesive used in surgical procedures. All egg manipulations were performed on a clean table frequently swiped with ethyl alcohol to minimize risk of egg contamination.
After injection eggs were cross-fostered by randomly assigning each egg to one of the triplet nests (Kruuk and Hadfield, 2007). Hence, this cross-fostering protocol ensured that all combinations of the experimental treatment and nest-of-origin were present in each nest-of-rearing. After crossfostering the eggs were returned to their assigned nests and left there for incubation. On the following days any additional eggs were treated similarly (transported to the laboratory, injected with a randomly chosen C/CORT/T4 solution, and returned to a nest); this protocol was stopped once on a given day incubation commencement was noted (i.e. eggs were not covered with nest material and warm). 1-2 days before the expected hatching date (11-12 days after the incubation start) all experimental nests were visited again. After verifying the development stage of eggs (by egg candling), all eggs were again gently collected and transported in a warmed box to the lab (leaving females with dummy eggs to ensure they would not desert their nests). There, they were placed in individual paper containers and put in an incubator set to 38 degrees and 70% relative humidity for hatching. From that moment the eggs were checked every hour. All chicks hatched between hours 0500 and 2000 were weighed, marked individually by nail clipping and taken back to their foster nests within 1 hour of hatching. Chicks hatched after 8 pm were left in the incubator until 5 am and brought to their nests the following day.
In spite of keeping all procedures as precise and as aseptic as possible, our manipulation had a measurable impact on the chicks' hatchability (likely resulting from water loss resulting from incompletely closed eggs, introduction of microorganisms interfering with embryos' development, or bursting of particularly small egg yolks after delivering additional 3 ul of liquid), which is usual in similar studies (Ruuskanen et al., 2009;Hsu et al., 2019;Sarraude et al., 2020). While natural hatchability in the studied population reaches 98.0%±0.5% (mean±SE, based on a random sample of 57 non-manipulated nests in 2014), the hatchability in nests manipulated by egg injections dropped to 51.0%±23%. Experimental groups differed slightly in hatchability (control group: 43.0%, CORT: 52%, T4: 58%; Fisher's exact test: p = 0.01), however this difference between groups injected with hormones (CORT vs. T4) was not significant (Fisher's exact test: p = 0.07). Average hatchability of ~50% corresponds especially well to similar experiments in birds with similarly sized eggs (Winter et al., 2013;Marri and Richner, 2014); small-egged species may be more sensitive to similar manipulations as their yolks are smaller (hence -more prone to irreversible damage), and their eggs contain smaller amounts of water (making them more prone to dehydration if shell puncture is not sealed completely).

Statistical analyses
To determine patterns of genetic variances in the offspring traits, we applied linear mixed models, fitting them to measured phenotypic traits: body weight measurements at days zero (hatching weight), 2, 8 and 14, tarsus length (on the 14 th day post-hatching) and PHA hypersensitivity response. In all models, response variables were assumed to be normally distributed (assumption checked by visually inspecting model residuals plotted against fitted values).
Each model contained fixed categorical effects of sex (males, females and unknown sex; females as intercept reference group), study year (2014-2016, 2014 as intercept reference group) and experimental treatment (control, CORT, T4, control as intercept reference group). Random effects included: nest-of-origin (genetic family) effect, nest-of-rearing (foster family) effect and residual error effect. Preliminary analyses included also a random term of the nest triad. It was subsequently removed from final models as it consistently explained marginally low amounts of variance (parameter consistently restricted to the parameter space boundary of zero) and its omission had no impact on other estimates.
All remaining random effects were structured to allow for hormone-specific estimates of relevant effects. In all cases the (co)variance structures were set to allow for heterogenous variances among the three experimental groups (3×3 covariance matrices). In the end we have fitted, for each response variable, a set of decreasingly complex models, testing various aspects of model (co)variance structures, starting with the most complex model (that included all possible variances and correlations) and simplifying it to remove redundant (co)variances. The order of tests made sure that factors possibly confounding the genetic variance (i.e. residual and rearing variances) were tested first. The only exception to this pattern were residual correlations between treatment groups which are not identifiable in our setup and thus were not estimated. In each successive step the more complex model (i.e. the alternative hypothesis, Ha) was tested against the simplified one (null-hypothesis; H0) using a likelihood-ratio test. Logged ratios of model likelihoods d = 2log(ℓHa/ℓHo) were assumed to be distributed as a mixture of χ 2 variates with k ∊ {{s, 1, 2, …, s+q} degrees of freedom (Self and Liang, 1987;Stram et al., 1994). The asymptotic distribution against which each likelihood-ratio statistic is tested is where s -the number of tested (co)variance parameters that lie inside the parameter space (e.g. correlations/covariances, where H0 = 0), q -number of tested (co)variance parameters restricted by H0 at the boundary of their parameter space (e.g. H0 = 0 for variances, H0 = 1 or -1 for correlations) (Self and Liang, 1987). For a simple df = 1 test of one variance component (H0: σ 2 = 0) this simplifies to rescaling the p-value of the resulting test by 0.5: . For each phenotypic trait we arrived at the simplest model supported by the likelihood ratio test. However, even if such models did not support heterogenous variances/unconstrained correlations, we estimated them from more complex models for illustrative purposes and plotting. Table 1 provides an overview of all types of fitted models, and (co)variance constraints involved.
In addition to single-trait models that focused on differences between experimental groups we have also fitted a multivariate model that included tarsus length and body weights at days 2, 8 and 14, to estimate between-trait genetic correlations across three experimental groups. In this model all (co)variance matrices were assumed to have a block-diagonal structure (i.e. did not allow for correlations between different traits measured in different experimental groups), and estimated all cross-trait correlations at all random effects' levels.
Estimated (co)variances were used to calculate heritabilities within experimental groups/traits (ratios of genetic variance to the sum of other variance components) of and genetic correlations between experimental groups/measured traits. Standard errors of heritabilities and genetic correlations were calculated using the delta method. Since the model we employed assumes the chicks in one nest of origin are full-siblings, the nest-of-origin effect estimates ½ of the total genetic variance in traits (which includes additive genetic variance, dominance variance if present, as well as variance generated by early maternal effects), and so heritability estimates we derive here are sensu lato. Dominance is assumed to be negligible in similar studies (Tolvanen et al., 2020;Class and Brommer, 2020). Maternal effects are also likely to be small (we cross-fostered eggs before the start of incubation, and so the only source of maternally derived variation could result from differential allocation of resources to eggs or from interactions between maternal and offspring genes; additionally, our hormonal manipulation likely offset any initial maternal contributions resulting from deposited steroids). The full-siblings assumption is also partly violated by small but significant proportion of extra-pair young detected in the studied population in selected breeding seasons (Arct et al., 2013), resulting in some of the offspring being actually maternal half-siblings. However, halfsibships would bias genetic variance estimates downwardly (i.e. conservatively), and in general result in small, negligible biases in genetic parameter estimates (Charmantier and Reale, 2005;Bérénos et al., 2014;Firth et al., 2015).
When reporting heritabilities and genetic correlations, the reported values are obtained from the best possible model (according to sequential LR tests), but -for illustrative purposes -with the reported component (e.g. genetic variances when reporting heritabilities, and genetic variances and correlations, when reporting genetic correlations) left unconstrained. This way of reporting ensures providing meaningful numbers instead of e.g. three identical values in models where no heterogeneity in genetic variances was detected.
All models were run in AsremlR v. 4.1 (Butler, 2019) in the R computational environment v. 3.0.1 (R Core Team and R Core Team, 2014). Before analysis we have removed from the data all individuals where the initial assignments to experimental groups were lost for some reason, e.g. because of premature hatching (and consequent failure to assign chicks to their respective experimental groups). Data to support all analyses is deposited in a repository using the Open Science Framework website (link TBA). Code to reproduce analyses is available through a GitHub repository (link TBA). Some variances or correlations constrained to be equal (here -variances for the C and T4 groups or correlations between C-CORT and CORT-T4) E.g. Ghr=1 + R + E implies heterogenous genetic variances (different VG in treatment groups), correlations constrained to be identical and equal to unity, neat-of-rearing and residual variances homogenous.

Fixed effects
Analyses included between 646 (2-days old chicks) and 603 (14-days old chicks) individuals (butsubstantially less in case of PHA responsivity). Means and standard deviations of raw data, together with sample sizes, are provided in Table S3. Body mass and tarsus length were sexually dimorphic (with males being heavier and larger, Table 2). Sex effect on body mass was not observed in 2 nd day nestlings, although males still tended to be larger in this age group. Hormonal manipulations exerted no statistically significant effect on the measured variables. Nevertheless, day 14 body mass of nestlings tended to be lowest in the CORT-manipulated group (Table 2). Other fixed effects had no statistically significant impact on the variables; in particular, we detected no confounding effect of the person measuring the tarsus length.

Random effects
In all variables we observed statistically significant levels of genetic variance (Table 3), resulting in heritabilities consistent with those reported elsewhere in the literature (Merilä and Fry, 1998;Hadfield et al., 2007;Drobniak et al., 2015;Perrier et al., 2018). For tarsus length the best supported model was the one showing a significant contribution of both nest-of-rearing and nest-of-origin effects, and a significant drop in genetic variance in the CORT-treated group, compared to the other two treatments (T4 and C; Table 3 and Table S1). When estimated separately for each experimental group, the genetic variance was highest in T4 and C, resulting in highest sensu lato heritabilities (h 2 ±SE: 0.39±0.13 and 0.33±0.12, respectively; Figure 2). Heritability in the CORT group was significantly lower (0.10±0.07; Figure 2).

Figure 2. Heritabilities of the five studied traits (symbol-coded) split by the treatment group (colour-coded). Heritability estimates come from models with fully heterogenous genetic variance structures. Whiskers represent SE (estimated via the delta method).
In all body mass variables, both the nest-of-rearing and nest-of-origin effects appeared significant, however there was no sign of any treatment-specific effect on the estimated heritabilities ( Figure 2). When calculated separately for each experimental group, heritabilities for day 2 body mass where the largest (C: 0.55±0.14, CORT: 0.86±0.14, T4: 0.69±0.13). They were of similar magnitude at day 8 (C: 0.46±0.12, CORT: 0.60±0.14, T4: 0.65±0.12) and day 14 (C: 0.38±0.14, CORT: 0.42±0.14, T4: 0.35±0.13). Similarly, there was no statistically significant treatment-specific trend in genetic variance in the PHA hypersensitivity response (C: 0.02±0.05, CORT: 0.37±0.22, T4: 0.53±0.24), with an exception of markedly smaller and statistically indistinguishable from zero genetic variance in the C group. Nevertheless, large standard errors (due to reduced sample size in this variable) prevented this pattern to be corroborated as significant heterogeneity in genetic variances in PHA response. In addition, two traits (body mass on day 2 and PHA response) exhibited significant, albeit relatively small, heterogeneity in residual variances (in both cases accounted for in the most supported models; Table S1). Detailed estimates of all variance components from models selected as best supported are presented in Table S1.
In all cases, genetic correlations between experimental groups were in all cases high and did not differ significantly from unity (which is the proper null hypothesis when testing such correlations within a trait at an inter-individual level) -hence, in all optimal models they are fixed at rG = 1. However, in few cases the differences were close to significance, with strong trends towards lower correlations. In body mass at day 14, although all cross-treatment correlations remained strongly positive, there was a drop (p = 0.07, Table 3) in the genetic correlation between the CORT and C groups (rG = 0.72±0.20) and between the CORT and T4 groups (rG = 0.52±0.19). More in-depth model comparisons indicated similar support when comparing a fully homogenous model (model G; equal variances, all rG = 1) to a fully unconstrained model (Ghr; p = 0.07). Comparing partly constrained model (Gr1,2≡r1,3) with a fully constrained one (G) indicated this trend was mostly driven by low rG between CORT and T4 groups (p = 0.1). Similar non-significant trend was visible in tarsus length (C-CORT rG = 0.39±0.13; T4-CORT rG = 0.52±0.27). In PHA response the correlation between C and CORT was negative not distinguishable from zero (rG = -0.33±0.73) -but it also was accompanied by a very wide SE and should be treated with reserve taken the small sample size of PHA-measured individuals. In none of the models did the rearing-nest associated correlations significantly differ from unity (Tables 3 and S1).

Covariance between traits
Our approach is underpowered to perform proper comparison of full multi-trait and reduced multi-trait models and so we decided to only illustrate the arising patterns using a full model allowing for all cross-trait correlations. Estimates from this model are summarised in Table S2. Figure 3 provides an illustration of the magnitudes of observed nest-of-origin correlations. An interesting pattern was visible for the genetic (nest-of-origin level) correlations: in case of body mass -in all experimental groups they tended to weaken for more distant age classes (Figure 3, middle row). Also, in hormonetreated groups cross-trait correlations tended to be relaxed compared to the control group, even to the point of reversing their sign and correlations becoming negative in the T4-treated group (middle row, last plot, Figure 3).  Rhr=1 + Ehr=0 8 448.15 2 0 * -(0) means cases where log-likelihood of the more complex model is slightly larger than of the simpler model, which effectively means they do not differ statistically; † -lower index indicates which group of parameters the df refers to (vvariances; c -covariances). Other shortcuts: VE, VG, VR -residual, genetic, rearing variance, respectively; rG, rRgenetic/rearing correlations.

Discussion
Differences in heritabilities measured in different biological contexts are not uncommon in natural populations. Apart from population-specificity of heritabilities (which contributes to marked variation in heritabilities even within one species (Lynch and Walsh, 1998)), a great deal of attention has been paid to changes in heritabilities (or more specifically -trait genetic variances) observed under varying biological conditions. Since genetic variance in a trait is one of the most fundamental ingredients of phenotypic evolution (Lande and Shannon, 1996;Walsh and Lynch, 2018), there is a long lasting consensus that changes to genetic variance "expressed" -or visible to natural selection -induced e.g. by individual characteristics or environmental variability should play important role in modulating the course of evolutionary change and conserving the levels of genetic variance in the wild (Kruuk et al., 2008). The latter is especially interesting: traits undergoing strong directional selection are expected to lose genetic variance -and with it their evolutionary potential (Kruuk et al., 2008;Walsh and Lynch, 2018); interactions of genetic variance with external cues may be one of mechanisms maintaining genetic variance in nature.
In this study, we have observed that one of the cues that may modulate observed levels of genetic variance (and consequently -heritability) are in ovo maternal hormones. By using a powerful approach of altering the hormonal environment of developing embryos, in conjunction with a crossfostering experiment that allows for effective separation of phenotypic variance contributors (Kruuk and Hadfield, 2007), we were able to show that steroid hormones acting early in individual development have the potential to alter the observed levels of expressed genetic variance, and to some extent the covariance patterns within and between traits. Of all studied traits, we have observed a strong effect of hormonal manipulation on genetic variance in the tarsus length, a trait typically exhibiting moderate to high levels of heritability in wild populations. Tarsus length heritability in our study (taking the control group as reference) agreed with other published estimates, including in the blue tit (Bonneaud et al., 2009;Teplitsky et al., 2009;Nilsson et al., 2009;Delahaie et al., 2017;Perrier et al., 2018), and with a more general set of published estimates on body size from wild populations (Postma, 2014) -and yet, in corticosterone treated nestlings heritability dropped to statistically indistinguishable from zero.
One insight into the reasons for observing this effect may be the physiological function of corticosterone. This steroid is usually considered to signal stress response and mediate organism's mobilisation after experiencing stress (Schoech et al., 2011;Haussmann et al., 2012;Mora et al., 2012). In birds CORT has been shown to mimic stress induced body weight variations and supressed immune response (Roberts et al., 2007), provoke development of stress-like phenotypes (Roulin et al., 2008), exacerbate behavioural differences along the shy-bold personality axis (Baugh et al., 2012), and impair parental care (Angelier et al., 2009) or learning behaviour (Kitaysky et al., 2003). If corticosterone exposure is regarded as mimicking stress exposure, CORT effects on genetic variance in offspring traits may mimic those observed in genotype-by-environment interactions, when individuals are exposed to stressful or unfavourable conditions. Although the impact of such conditions may differ, the often observed pattern is reduction in observed levels of genetic variance (Hoffmann and Merilä, 1999), supported by a meta-analysis of such results (Charmantier and Garant, 2005). Hoffmann and Merilä (1999) argued that one possible mechanism of such reduction is stopping of offspring growth, caused by stress, before inter-individual variance in achieved body size can fully develop. In this study this explanation doesn't seem to be valid: we haven't observed any systematic differences in body size between the three treatment groups. Other mechanism that can be invoked in cases of condition-varying genetic variances involves changes in gene expression at a molecular level (Hodgins-Davis and Townsend, 2009). Sparse evidence suggests, that in ovo CORT can induce gene expression changes in birds (Ahmed et al., 2016), but such evidence is not unambiguous (Lutyk et al., 2017). Function of corticosteroids' receptors as transcription factors directly modulating expression of certain genes is known (Kawata, 1995;Baker, 1997). Nonetheless, more work is needed -especially at the very basic, molecular level and in early developing embryos -to verify whether transcription modulation could contribute to population-level magnitudes of quantitative genetic parameters such as genetic variance.
In contrast to corticosterone manipulation, testosterone-exposed nestlings did not exhibit any significant changes in genetic variances of their traits. Physiologically, testosterone is traditionally associated with sex-specific effects and is assumed to mediate trade-offs between body maintenance and resource use in production of secondary sexual characters (Peters, 2007;Kingma et al., 2009). Exposure to testosterone early in development was shown to stimulate development of sexually selected traits, bias sex-ratios towards males and increase dominance and competitive behaviours . Unfortunately, due to sample size we couldn't robustly test for sex-bytreatment interaction in genetic parameters. Nevertheless, due to many reported cases of sex-specific genetic variances in a number of traits (Wyman and Rowe, 2014) we expected to see significant impact of hormonal manipulation in our study. One reason for not seeing such effect may be the use of only testosterone (T4), in contrast to many similar studies using both testosterone and androstenedione as two major sex-linked hormones . Our observation is similar to other studies where no significant in ovo testosterone impact was noted (Tschirren, Saladin, et al., 2005) and suggests that future studies should look more closely on complexes of similarly acting hormones, applying them together to better mimic biological reality.
Other studied traits did not exhibit significant patterns of genetic variance with relation to hormonal manipulations. Estimated heritabilities of all body weight variables were higher than in the case of tarsus length -but still remained within ranges known from published estimates, including estimates produced using an animal model rather than full-sib analysis (Jensen et al., 2003;Postma, 2014;Garant, 2020). One explanation for higher heritabilities may lie in reduced hatchability: hatching failure of some eggs may have reduced intra-clutch competition for resources, effectively equalizing prospects of competing offspring and reducing environmental sources of trait variation in traits most closely related to body condition. Body mass would be particularly sensitive to such effects as it is more condition-sensitive than structural body size (Dubiec et al., 2006). Especially interesting is the high heritability of body weight on the 2 nd day of nestling life. Our experiment involved transferring developing eggs to an incubator, and hence removing possible sources of early posthatching body size differences, resulting from changes in pre-hatching incubation behaviour and asynchronous hatching (tits hatchlings rarely hatch at once, usually the process takes several hours to 2 days). Nevertheless, taking into account 2 nd day body mass' standard errors, this estimate also was in line with many published values (including examples of a similar decrease in nestling body weight heritability with age (Bonneaud et al., 2009;Postma, 2014)). Heritability of PHA response was in line with estimates known previously from the same population, generated using an animal model (Drobniak et al., 2015). Lack of the effect of the treatment on the body mass and PHA response's genetic parameters is interesting on its own: it emphasizes that different traits (even partially linked, like body weight and tarsus length) can exhibit very contrasting responses to early developmental conditions. Also, body weight in birds is often interpreted as describing more than just individual body size: usually large variation in residual body weights (when regressed against individual body sizes) emphasizes its complex nature and the fact that body weight should be regarded more as a body condition measure rather than simply body size proxy (Labocha and Hayes, 2012;Janas et al., 2018).
What actually prompted the idea for our experiment in the first place was the observation of significant, positive or negative, correlations between direct and maternal genetic effects from wild populations (McAdam et al., 2002;Wilson et al., 2005;Galloway et al., 2009). It is likely that maternal -and other indirect -genetic effects often are involved in complex covariances with direct genetic effects, contributing ultimately to the potential of traits to evolve (McAdam et al., 2014). It is still too early to conclude that patterns seen in our study may be one of the driving forces of such maternal-genetic covariances. What we know for sure is that maternal deposition of steroid hormones into bird eggs has heritable underpinnings , and so steroids could serve as a possible mechanisms generating maternal-genetic effects correlations. Coupled with a clear role of maternally transmitted steroids in signalling to the offspring the state of external environment (Tschirren et al., 2004;Groothuis and Schwabl, 2008), and with our results, available evidence offers a new perspective on the origin of genotype-by-environment interactions, one that may depend on maternal effects.
On the methodological side, our study may suffer from several shortcomings. Firstly, we based our analysis on full-sibling analyses, which poses a risk of upwardly biasing our genetic variance estimates with potential dominance and early maternal effects. For reason outlined in the Materials and methods section we do not expect this bias to be significant; also our estimates agree with published values based on more robust approaches such as the animal model (Postma, 2014). One of the confounding effects that could theoretically contribute to the observed pattern was hatching success and the fact it was significantly linked to the experimental group, with control group having the lowest hatchability compared to the hormone-treated groups. Smaller embryos generally hatch less successfully (Krist, 2011), which could negatively select for hatchlings' size and eventually reduce body size variation. However, if this would be the case, we would expect drops in variance components in the control group and not in hormone-treated groups -and similar reduction would be expected in all body size-linked variables. Lowered hatchability likely also reduced intra-clutch competition for food compared to broods of natural size. This in principle could also reduce variability in body mass -but should not mask genetic variation and would manifest itself mostly in environmental/residual components of inter-individual variation. Finally, our study proved to be underpowered to detect cross-treatment differences in genetic correlations between different traits, as we could not formally confirm the observed patterns using likelihood-ratio tests. Nevertheless, observed patterns are interesting: according to analyses, testosterone (T4 group) appears to decouple genetic effects influencing tarsus length and body weight in the blue tit, by lowering the magnitude of cross-trait correlations linking these traits. Although we are aware of the limitations of this estimation, it emphasizes the need for more multivariate patterns being studied in similar physiological or hormonal contexts. Univariate analyses may often fail to detect changes in traits' evolutionary potential in cases when only between-trait covariance structures are affected.
In summary, our study provides the first experimental attempt to identifying mechanisms that may be responsible for the modulation of the expressed genetic variance in nature. Complementary studies in wild populations replicating our results, as well us more in-depth analyses looking at the molecular pathways involved, are needed in order to better understand and explain the mechanism of how steroid hormones may mediate the observed effects. Taken the widespread occurrence of environmental effects on the levels of genetic variance in nature, and the commonly accepted role of steroid hormones in mediating environment-induced influences during development, we believe that our study provides a refreshing and novel perspective on the issue.