Carotenoid pigments produce most red, orange, and yellow colours in vertebrates. This coloration can serve as an honest signal of quality that mediates social and mating interactions, but our understanding of the underlying mechanisms that control carotenoid signal production, including how different physiological pathways interact to shape and maintain these signals, remains incomplete. We investigated the role of testosterone in mediating gene expression associated with a red plumage sexual signal in red-backed fairywrens (Malurus melanocephalus). In this species, males within a single population can flexibly produce either red/black nuptial plumage or female-like brown plumage. Combining correlational analyses with a field-based testosterone implant experiment and qPCR, we show that testosterone mediates expression of carotenoid-based plumage in part by regulating expression of CYP2J19, a ketolase gene associated with ketocarotenoid metabolism and pigmentation in birds. This is the first time that hormonal regulation of a specific genetic locus has been linked to carotenoid production in a natural context, revealing how endocrine mechanisms produce sexual signals that shape reproductive success.

(a) Plasma sample collection and quantifying circulating carotenoids

            We collected samples and conducted experiments (below) on free-living red-backed fairywrens captured in mist nets at our long-term study site in Samsonvale, QLD, Australia (27°27’ S, 152°85’ E). We collected blood samples (20–70 µl) from the wing vein using heparinized microcapillary tubes from May - August 2017 and 2018 during the non-breeding season, a period when most males were actively moulting into their breeding season plumage. Blood was centrifuged for 5 min at 10,000 rpm, after which the plasma was separated from the packed cells and stored at -20° C until transport to the US where samples were stored at -80° C until high performance liquid chromatography (HPLC) analysis.

            Plumage score was recorded at time of capture following prior methods (Karubian, 2002), yielding total ornamentation scores that ranged from 0 (brown) to 100 (red/black). Based on this score, males were labelled as either “unornamented” (brown plumage, score <33), “intermediate” (mixed plumage, score between 33 and 66), or “ornamented” (red/black plumage, score >66). Females always have completely brown plumage (plumage score = 0) and are therefore considered unornamented. Timing of moult into ornamented plumage is variable (Welklin et al. unpub. data), similar to other recorded Malurus species (Mulder & Magrath, 1994), meaning that male plumage at time of capture and sample collection may differ from the “final” plumage colour score the male expressed later in the breeding season. Males can breed in unornamented or ornamented  plumage or serve as auxiliaries with unornamented plumage (helpers) at the nest (Webster, Karubian, & Schwabl, 2010). Because we were interested in differences between unornamented and ornamented plumage, we documented the “final” plumage score of colour-banded individuals on November 1^st (the approximate mid-point of the breeding season) and used this “final” score for analyses of unornamented males, ornamented males, and females; we excluded from analysis the relatively small sub-set of birds whose “final” plumage score was intermediate (n = 11). We assigned either minimum or known age (age range 1–7 years, 77%, 123 of 160, were of known age) to all birds at the time of sample collection using nestling banding records or extent of skull ossification (ossification scale modified from (Pyle, Howell, Yunick, & Desante, 1987), and we have validated this scale within this species multiple times). Qualitatively similar results were obtained in analyses run with these age criteria, or using only known-age birds (see electronic supplementary material, Table S1).

            We used high performance liquid chromatography (HPLC) to identify and quantify the concentration of carotenoids in the plasma, following the methods of Rowe and McGraw 2008 (Rowe & McGraw, 2008). We analysed carotenoids in 160 plasma samples (n=42 females, n=29 unornamented males, n=89 ornamented males).  To assess the relationship between circulating ketocarotenoid levels and plumage phenotype, we ran a linear mixed effect model with the lme function in the R package nlme (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2018), R v. 3.6.0 (R Core Team, 2019). The model included total circulating ketocarotenoid concentration (i.e. the sum of alpha-doradexanthin, astaxanthin, adonirubin, and canthaxanthin concentrations) as the response variable and the following predictor variables: (1) phenotype (female vs. unornamented male vs. ornamented male); (2) age (as a continuous variable); (3) year of sample collection; and (4) the interaction between phenotype and age. To control for repeated measures of the same individual across years, we added individual as a random effect (n=13 individuals sampled both years). Residuals were inspected visually for homoscedasticity, and we used the varIdent function to control for heterogeneity of variance between groups. Year of sample collection did not improve model fit (i.e, it did not improve AIC by more than 2 and the p-value of the Year variable was greater than 0.05) and was therefore dropped. We tested the model for significance of phenotype with a Tukey’s posthoc test using the glht function in the R package multcomp (Hothorn, Bretz, & Westfall, 2008).   

(b) Testosterone implantation and liver sample collection

            We collected liver samples from breeding, but not auxiliary helper, red-backed fairywrens in November 2017, to control for potentially confounding underlying differences in endocrine or genetic profiles that may exist between auxiliary non-breeding vs. breeding individuals (Lindsay, Webster, Varian, & Schwabl, 2009). First, three breeding unornamented males were implanted with testosterone. At time of initial capture and implantation, around 10 feathers were plucked from the centre of the back to induce feather replacement at that location. Implants were composed of beeswax (73% by weight; Sigma-Aldrich, St. Louis, MO) and hardened frozen peanut oil (24% by weight; ACROS Organics, NJ, USA) that were mixed in a water bath at 67° C. Once the beeswax/peanut oil mixture was melted, crystalline testosterone (3% by weight; Sigma-Aldrich, St. Louis, MO, USA) was dissolved in 2.5 μl of 200 proof ethanol (Fisher Bioreagants^TM). The implants were formed by feeding partially solidified wax through the tip of a syringe, resulting in implants of 2 x3.2mm weighing between 19.8 and 20.7 mg. Testosterone concentration in the beeswax carrier was scaled to produce high physiological concentrations found in circulation during the breeding season (Lindsay et al., 2009). Implants were inserted subcutaneously using forceps above the thigh into a small (2-3 mm) skin incision that was sealed with veterinary skin adhesive. After confirming that the incision was completely sealed and the bird was in good condition, the bird was released. Three unornamented males were implanted with sham controls (beeswax/peanut oil implant with no testosterone) and also had around 10 back feathers plucked.

            Implanted birds were recaptured 10-12 days post-implantation for liver sample collection, a time period that allowed for growth of pin feathers in the plucked plumage patches (red pins in testosterone-implanted males, and a mix of red and brown pins in the sham-implanted males, consistent with what had been observed in another feather-plucking study in this species (Karubian, Lindsay, Schwabl, & Webster, 2011)). We were unable to recapture one of the sham-implanted birds after implantation. Instead, we captured and obtained samples from one additional breeding unornamented male, who was not implanted, to include in our control group. Following re-capture on territories in mist nets, birds were immediately sacrificed by cervical dislocation. Body dissection was performed in the field, and the right lower lobe of the liver was removed and stored in 1mL of RNAlater storage buffer (ThermoFisher Scientific), and immediately placed on dry ice. Samples were stored at -80° C until RNA extraction. In addition to the three testosterone-implanted unornamented males and the three control unornamented males (two with sham implants, one without an implant), we also collected liver samples for three ornamented breeding males and three breeding females (without implants). All birds sacrificed were seen paired with a male or female within two weeks prior to sample collection, and all samples were collected within a period of 9 days. Circulating androgens were measured using an established radioimmunoassay protocol for this species (full methods in (Barron, Webster, & Schwabl, 2015; Lindsay et al., 2009)); the intra-assay coefficient of variation was 8.68%. Testosterone-implanted birds were confirmed to have high concentrations of circulating androgens at time of collection (mean=3027 pg/ml, range=2198-4065 pg/ml), which is within the natural range of androgens for breeding ornamented males in this species (Lindsay et al., 2009). Ornamented males also had similarly high levels of androgens at time of collection (mean=1834 pg/ml, range=1124-2925 pg/ml). We were unable to obtain samples to assay testosterone concentrations for control unornamented males or females.

(c) Quantifying relative expression of CYP2J19

            To extract messenger RNA (mRNA), we removed liver tissue from the RNAlater buffer and homogenized it in a Qiagen TissueRuptor. We used a Qiagen RNAeasy mini-kit, following manufacturer’s instructions, and reverse transcribed the mRNA to cDNA with a Superscript IV first strand synthesis kit (Invitrogen). All qPCR reactions were run on CFX96 Touch™ Real-Time PCR Detection System (BioRad) with CFX Maestro Software (BioRad), using PowerUp SYBR Green Master Mix (Thermofisher Scientific). For measurements of CYP2J19 expression, we used qPCR primers CYP2J2-2F and CYP2J2-2R (Mundy et al., 2016). We assayed gene expression in triplicate for each sample and normalized the data using the housekeeping gene GAPDH, using primers Gg_GAPDH_qPCR_F and Gg_GAPDH_qPCR_R (Lopes et al., 2016). Reaction conditions for qPCR were tested and optimized using a standard curve produced by creating a serial dilution of a pool of all cDNA samples. Efficiencies ranged from 95%-105%, and we analysed qPCR data using the delta-delta Ct method (Livak & Schmittgen, 2001), further described in electronic supplementary material, Methods.

            We confirmed that there was no effect of presence or absence of the sham implant on gene expression within unornamented males (see electronic supplementary material, Methods), and homoscedasticity was confirmed with a Breusch Pagan test (p>0.05). We tested for statistical differences in liver CYP2J19 expression (log fold change) between phenotypes with an ANOVA, using the aov function in R, followed by a Tukey’s posthoc test using the TukeyHSD function in R.

Supplementary Methods:

Analyzing qPCR data with delta-delta CT method

We analysed qPCR data using the delta-delta Ct method [S1], which reports gene expression (mRNA abundance) for the gene of interest as the fold change in expression (2^-DDCt), normalized to a housekeeping gene and calibrated to a “calibrator sample”. We use GAPDH as our housekeeping gene, as it is often used as a housekeeping gene for qPCR analysis in birds [S2–S5]. We averaged the DCt (Ct CYP2J19 – Ct GAPDH) for all three females and used that average as our calibrator sample. We present our results as log fold change. Using the average female DCt value as our calibrator sets the average log fold change value for females at zero, allowing for easier visualization of the difference in relative gene expression between phenotypes. We also ran this analysis with only one female as our calibrator sample, and our statistical results remain exactly the same since these are all relative expression levels.

Assessing the effect of the presence of absence of a sham implant

In order to assess whether or not having a sham implant affected expression levels for those unornamented males, we combined the two sham-implanted males and the unmanipulated unornamented male into one phenotype category we called “control unornamented male.” We ran a linear model using the lm function in R with CYP2J19 expression (log fold change) as the response variable, and presence of an implant (yes vs. no) nested within phenotype (female vs. control unornamented male vs. testosterone-implanted unornamented male vs. ornamented male) as the response variable. We found no interaction between phenotype and implant, and the implant variable did not improve the fit from the model, suggesting there was no significant effect of not having the sham implant within our control unornamented males. In addition, we evaluated raw data and confirmed that the log fold change for the unmanipulated male (2.20) was similar to that of the two sham implant males (1.95 and 2.19, and can be seen in figure S2), and the variance of the log fold change for control unornamented males (all three males) was very small, as can be seen in the standard error bars in Figure 2. Taken all together, we therefore dropped the implant variable from our model.

The first data set “Khalil_etal_2020_HPLC_data.csv” is the circulating carotenoid data from the High Performance Liquid Chromatography (HPLC) Analysis. This data was used for the circulating kecarotenoid model.

The data set consists of (in this order):

• Year of sample collection
• Sample number
• Day of sample collection
• ABBBS number – metal band number, a unique identifier
• Color bands – color band combination, another unique identifier, and the one used in the model for the random effect
• Sex
• Age - those with a “+” were of minimum age, since they are considered at least that age.       Those without a “+” are of known age.
• Age_min – same information as age, but without the “+” to be used in analysis. When running the model with only known age birds (Table S1), we filtered out individuals that had a “+” in the Age column.
• The next six columns are the circulating carotenoids -ador (alpha-doradexanthin), asta (astaxanthin), lut (lutein), adoni (adonirubin), zea (zeaxanthin), and canth (canthaxanthin).
• Breeding_phenotype – the breeding phenotype (as of November 1^st) used for the analysis. F is female, DM is unornamented male, and BM is ornamented male.

The second data set “Khalil_etal_2020_qpcr_data.xlsx” is the qPCR ct data from the CYP2J19 and GAPDH expression analysis. This data was used to calculate log fold change using the delta-delta ct method described in text.

The data set consists of three sheets:

1) gapdh_ct – this is the ct data using the GAPDH primers. This includes

• The well the sample was in
• The cq
• The sample number
• The phenotype of the sample. F is female, T-DM is testosterone-implanted unornamented male, DM is control unornamented male, and BM is ornamented male.
• The average GAPDH ct, calculated from the cq column as each sample was run in triplicate

2) cyp2j19_ct – this is the ct data using the CYP2J19 primers. This includes the same type of information as the gapdh_ct sheet

3) delta_delta_ct_math – this is all the math to calculate the log fold change using the delta delta ct method. This includes:

• The sample number
• The phenotype
• If the control unornamented male (DM) had a sham implant or no implant
• The average GAPDH ct from sheet 1
• The average CYP2J19 cqtfrom sheet 2
• The delta ct (CYP2J19 ct – GAPDH ct)
• The average delta ct for females. This is used as our “calibrator” sample, as described in the supplementary methods
• Delta delta ct (delta ct – average female delta ct)
• Fold change (2^delta delta ct)
• Log fold change. This value is used in the ANOVA and in figures.

4) Intra-assay CV - this is the math for how intra-assay coefficient of variation for the GAPDH and CYP2J19 qPCR assays