Estimation of maternal expression trait values
Prior to conducting further analyses using the initial 47,651 transcripts obtained after a series of standard filtering steps (Methods S1), we normalised read count data for varying library sizes using the trimmed mean of M values (TMM) approach as implemented with thecalcNormFactors function from the edgeR v.3.14.0 package (Robinson et al ., 2010) in R. We further filtered transcripts that were not expressed across all three siblings nested within the same maternal tree in either garden. We implemented a linear mixed effects model for each of the remaining 27,441 transcripts (Eq. 1):
\(Y_{\text{ijklmn}}=\ \mu+\text{Batch}_{j}+\ \text{SY}_{k}+\text{Garden}_{n}|\left(\frac{\text{Pop}_{l}}{\text{Fam}_{m}}\right)+\ \epsilon_{\text{ijklnm}}\)(1)
Here, \(Y_{\text{ijklmn}}\) represents the normalised expression value of a transcript for the i th seedling in thej th sequencing batch sowed in thek th year from the m thmaternal tree nested within the l th population at n th garden. The global mean across all seedings for a transcript is given by \(\mu\). The terms Batchand SY are fixed effects and represent the date of sequencing and year of seed sowing, respectively. The term Fam represents the maternal tree that is nested within Pop , which represents the source population of the seeds for a maternal tree. Both Fam andPop are treated as random effects that are allowed to vary by the fixed effect of Garden . This linear mixed effect model was fitted using the fitVarPartModel function with log-transformed normalised counts and precision weights, which are ideal for RNAseq datasets (Law et al ., 2014), obtained from thevoomWithDreamWeights function in the variancePartition v.1.16.1 package (Hoffman & Schadt, 2016) in R. Using the fitted model in equation (1), we obtained garden-specific genetic values for maternal trees and populations for each transcript.