Parallelism in DE gene identity and expression direction between lineages
To ask whether parallel phenotypic changes across independent lineages rely on shared gene expression changes, we evaluated overlap in DE transcript identity and expression direction between the Aripo and Quare lineage datasets. We used chi-square tests of independence to test for greater than chance overlap between drainages in those gene sets differentially expressed based on rearing environment or population of origin within drainages, as well as to test for biases in concordant versus non-concordant expression between drainages. We performed two additional comparisons to address potential bias in our conclusions due to high false negative rates in gene expression data (Rice, Schork, & Rao, 2008). First, we performed the same analysis as before at two less stringent p-value cut-offs (p=0.1 and 0.2). Second, for comparisons of expression direction, we made an additional, more conservative comparison by including not only those genes differentially expressed in both drainages (i.e. the intersection of DE gene lists), but also those differentially expressed in eitherdrainage (i.e. the union of DE gene lists) in the analysis. We again used log-fold expression differences to call gene expression as being in the same/concordant or opposite/non-concordant direction between drainages.
We also performed one addition comparison to assess parallelism in a subset of the population DE genes likely to arise from selection rather than other evolutionary processes, such as drift. We identified to those genes mostly likely to have diverged by selection, rather than by drift, by calculating PST (a measure of phenotypic divergence between populations) from phenotypic variance components and comparing it to published estimates of FST (a measure of neutral genetic divergence between populations) for the focal populations, using the method of Leinonen, Cano, Mäkinen, & Merilä (2006). For these calculations, we made the conservative assumption that heritability (h2 ) of transcript abundance = 0.5 as we have done previously (Ghalambor et al., 2015), and obtained estimates of within and between population variance from general linear mixed models (see Supplemental Methods). We then restricted DE gene sets in both lineages to those in which PST > FST and compared the number of overlapping genes and their relative expression direction using chi-square tests as above. We used FST estimates for our focal populations (Aripo = 0.225; Quare = 0.340) reported by Willing et al. (2010).