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.