Modeling the effect of sex on gene expression
We used a mixed modeling approach with the R package EMMREML to quantify the effect of sex on gene expression while controlling for sample collection year, RNA extraction date, RNA concentration, and RNA quality (Akdemir & Godfrey, 2015). We used an identity matrix as the known covariance structure which is required for the EMMREML mixed modeling approach. To focus on gene expression of this putatively sexually selected trait, we only included adults in this model (n =14 adult males, n =14 adult females). We calculated the false discovery rate (FDR) for each gene using the R package qvalue(Storey et al., 2022). Genes that passed a threshold of a FDR of 20% were considered differentially expressed between the sexes: “male-biased” if they were more highly expressed in males and “female-biased” if they were more highly expressed in females. We then Z-transformed expression values for this set of male-biased and female-biased genes across all 36 individuals (n= 14 adult males, n =14 adult females, n =6 subadult males, n =2 subadult females) and averaged the Z-transformed expression levels of these genes per individual to obtain a composite sex-biased gene expression score for each individual. Finally, we ran a linear regression model with expression score as the outcome variable and age category (adult male, adult female, subadult male, subadult female) as the predictor variable to investigate whether subadults differed from adults in sex-biased gene expression.