Pooled whole‐genome sequencing of interspecific chestnut (Castanea) hybrids reveals loci associated with differences in caching behavior of fox squirrels (Sciurus niger L.)

2.1Seed collection

Seeds were collected in late September and early October from a planting of several hundred BC3 ([(Castanea mollissima × dentata) × dentata] × dentata) × dentata) chestnuts at Purdue University's Lugar Farm in Tippecanoe County, IN. The purpose of the planting is blight resistance screening. Most seed parents were 11 years old, but several were 4 years old in 2014 at the start of the study. “Clapper,” a BC1 tree, was the blight‐resistance donor and only source of Cm genetic material in this backcross population. Cm nuts were obtained from a pair of trees planted as blight‐resistant checks in the Lugar Farm orchards. Cd nuts were obtained from two adult trees growing at the Purdue Wildlife Area in Tippecanoe County, IN. BC3 seed parents were chosen based on seed size, with roughly equal numbers of large‐seeded, small‐seeded, and average‐seeded trees chosen, in order to capture a wide range of phenotypic variation. Seed parents were tagged with durable individual plastic nursery labels, but due to large annual variation in the size of seed crops and the loss of some seed parents due to chestnut blight, different seed parents were chosen each year of the study. Seeds were collected by knocking burrs off the parent tree using a ~2‐m wooden pole and manually removing seeds from the bur if necessary. Seeds were floated in water to determine viability; floating seeds were deemed nonviable and discarded.

2.2Seed measurements

Seeds were stored in a cooler (4.4°C) following de‐burring and floating and stratified in peat moss to maintain viability during cold storage. In October, at least 10 seeds from each seed parent were weighed on a digital scale to determine average seed mass. Length (from seed base to tip) and width (across the broadest part of the seed) were determined using digital calipers. In 2015, desiccation was also measured by weighing seeds immediately after collection and again 80 days following collection.

2.3Seed tagging

Tagging was carried out immediately before dispersal trials to avoid spoilage of seeds. A method similar to that employed by Xiao, Jansen, and Zhang (2006) and Hirsch, Kays, and Jansen (2012) was used. A hole was made in the proximal (wider) end of each seed using either a botanical dissecting needle or a small (~2 mm) drill bit. A piece of 24‐gauge green floral wire approximately 12 cm long was looped through the hole and twisted to secure it. A piece of brightly colored waterproof tape was attached to the end of the wire and labeled with a number designating the seed parent.

2.4Dispersal trials

Dispersal trials were conducted in late October, November, and December of each year at four feeding stations placed in and around the Lugar Farm chestnut plantings in Tippecanoe County, IN in a manner similar to the methods of Lichti, Steele, Zhang, and Swihart (2014). In 2016, a feeding station in Woodford County, IL adjacent to the campus of Eureka College was added. At both locations, fruiting chestnuts were present in addition to black walnut (Juglans nigra) and several oak species. Fox squirrels (Sciurus niger) were the only scatter‐hoarding squirrel species observed at either feeding station. Feeding stations were prebaited to acclimate local squirrels to the feeding locations in August‐September prior to dispersal trials. During dispersal trials, 10 (2016) or 25 (2014–15) seeds from 5 to 6 (2016) or 3 to 4 (2014‐15) parent trees were randomly distributed near a post at the center of each feeding site. Seeds were left out for 4‐5 days, and seed fates (cached, consumed, or left at feeding station) were recorded and dispersal distances measured with a forestry measuring tape attached to the post at the center of the feeding site. Intensive searches for seeds were conducted up to 20 m from the feeding site, although some seeds were found outside this distance due to the high visibility of the tags. Trials started in late October or early November and continued through December until the soil surface froze. Relationships between seed dimensions and dispersal parameters were statistically investigated using the lm and glm packages of R software version 3.2.3 (R Core Team 2015).

2.5DNA Isolation

DNA was isolated from BC3 seed parent trees following dispersal trials. Dormant twigs were collected for DNA extraction in early spring 2016 and 2017. Terminal sections (about 3–5 cm) of first‐year twigs were ground to a fine powder in liquid nitrogen using a mortar and pestle. The ground tissue was placed in 5 ml of heated (50°C) CTAB extraction buffer in a 15‐ml conical tube and incubated 4‐8 hr at 50°C. Following incubation, 1 ml of 20 mg/ml proteinase K solution was added and samples were incubated for an additional 15 min. Five milliliters of 25:24:1 phenol:chloroform solution was added, and samples were purified using a standard phenol:chloroform extraction (Doyle and Doyle 1987) followed by precipitation of DNA using 0.2 M sodium chloride and isopropanol. After pelletting and resuspending samples in TE buffer, contaminants were removed using Zymo Research OneStep PCR Inhibitor Removal kits (Zymo Research). Following purification, samples were quantified using a Nanodrop 8000 (ThermoFisher Scientific), and 2% agarose gel, pooled according to dispersal parameters, then submitted to the Purdue Genomics Core Facility for sequencing.

2.6DNA pooling and sequencing

Pools of samples were made for different phenotypic classes based on (a) mean dispersal distance for cached seeds and (b) frequency of caching. The strong‐dispersal pool (Pool A; eight samples) contained DNA from parents that produced seeds with a long dispersal distance (>5 m average dispersal distance for cached seeds) and high frequency of caching (5%–59%) of recovered seeds in caches), including one Chinese chestnut. A moderate‐dispersal pool (Pool B; seven samples) contained parents that produced seeds with a shorter dispersal distance (<5 m) and low frequency (4%–14%) of caching. The weak‐dispersal pool (Pool C; 10 samples) contained seed parents that produced seeds with a frequency of caching and dispersal distance near 0, including one American chestnut (Table 1). Twenty microliters of DNA at concentration 200 ng/μl from each individual sample was included in a pool, and the combined sample was submitted for library construction and sequencing; samples were sequenced as separate libraries on one Illumina HiSeq 2500 (Illumina Inc., San Diego, CA, USA) lane. Reads were paired‐end, 100 bp in length. The individual genomes of “Clapper,” several unrelated Chinese chestnuts, and two American chestnuts were sequenced separately with two samples per lane as part of a blight resistance study (LaBonte et al., in preparation).

2.7Genome assembly and SNP calling

Short reads were assembled to the draft Chinese chestnut reference genome v1.1 (Carlson et al., 2014; Pereira‐Lorenzo et al., 2016) using the Burrows‐Wheeler aligner (bwa) (Li and Durbin 2009). Alignments were processed and polymorphisms called for each pool of samples using Picard Tools and the Genome Analysis ToolKit (GATK) best practices workflow (DePristo et al. 2011; Van der Auwera et al. 2013), minus the quality score recalibration step. When calling SNPs using the HaplotypeCaller tool in the GATK, ploidy was set to twice the number of individuals in the pool.

2.8Analysis of SNP data

Several custom Perl scripts were used for processing of polyploid SNP data files generated by the GATK pipeline. The goal of these scripts was to discriminate between predicted genes that had two genotypes well‐represented in a pool (Cm/Cd sites) and genes that had a single genotype fixed or nearly fixed within a pool (Cd/Cd sites) for SNPs within the predicted gene sequence, informed by the observation (LaBonte et al. 2018, in preparation) that heterozygosity values in coding sequences of chestnut hybrids are usually much higher than either parent species. There are no Cm/Cm sites in a BC3 genome. To identify Cm/Cd sites in pooled BC3 genomes, the multi‐sample SNP file was filtered for polymorphisms occurring within predicted genes (AUGUSTUS gene prediction; Stanke et al. 2006) with strong (e‐value <0.001) alignments to the curated Uniprot/SwissProt protein database. For each SNP (minimum depth = 8) between the predicted transcription start and stop sites of each gene, a heuristic hybridity estimator (HE) was assigned to approximate the proportion of individuals that were heterozygotes (Cm/Cd) within a pool. If major allele frequency for a SNP within a pool was between 0.45 and 0.55, HE was assigned a value of 0.75 (most individuals heterozygous), if it was between 0.55 and 0.70, HE was 0.5, if it was between 0.70 and 0.85, HE was assigned the value 0.25, and if the major AF was >0.85, HE was assigned 0 (all or nearly all individuals homozygous). This estimate was averaged across all the SNPs in each predicted gene sequence. When more than 50% of the SNPs in a genotype's sequence could not be scored, that genotype‐locus combination was considered missing data. Finally, HE was averaged across 10‐gene bins for each pool and compared, by bin, among pools. Binning was carried out to detect linked Cm haplotype blocks and avoid detecting single outlier genes as false‐positives. Loci potentially contributing to interspecific differences in seed dispersal were identified as those 10‐gene bins that had a difference in average HE values >2 standard deviations greater than the average HE difference between the strong‐dispersal pool and the moderate‐ and weak‐dispersal pools. Genes within these regions were annotated using the UniProt entries for aligned proteins from the UniProt KB/ SwissProt database (Xenarios 2016). Predicted molecular interactions were analyzed using the STRING protein database.

All the BC3 trees in our sample inherited 100% of their Cm alleles from “Clapper,” but only regions segregating in “Clapper” (loci with Cd/Cm genotypes) were informative, that is, about ~50% of the genome of “Clapper” (a BC1 tree), the other half were Cd/Cd. Thus, our inference on the genomic basis of interspecific difference in seed dispersal was limited to those heterozygous regions. Of loci that were hybrid in “Clapper,” any given BC3 descendant of “Clapper” was expected to retain a Cm allele at one in four loci, with the rest acquiring a second Cd allele in two rounds of meiosis. Therefore, at a locus known to have a Cd/Cm genotype in “Clapper,” a random sample of “Clapper”‐derived BC3s is expected to have one Cm allele observed out of every eight (one Cm/Cd and three Cd/Cd). As we genotyped BC3s in pools, opportunities for random sampling error were present; a given individual's genotype might be over‐represented at a locus, biasing allele frequency. Over‐representation of an individual could be due to differences in DNA quality, inaccurate estimates of DNA concentration prior to pooling or random inclusion of more DNA fragments from one individual during the high‐throughput sequencing process. We developed a Perl script to estimate the likelihood that more Cm alleles than expected by chance alone were present at a given SNP locus in the pooled data.

First, a panel of eight Cm genomes with no evidence of hybrid background, two Cd, and “Clapper” whole‐genome sequences (LaBonte et al. 2018) were used to filter a the pooled whole‐genome SNP file for loci with one allele fixed in Cd, one allele fixed in Cm, and a Cm/Cd genotype in “Clapper.” The coordinates of these loci were recorded as informative SNPs because markers at those loci allowed us to make inferences about the effects of Cm alleles on seed dispersal. Only informative SNPs were kept from the pooled genome SNP file for the analysis.

Next, the program made random draws from arrays of 100 binary values (0 for Cd, 1 for Cm) set to represent the expected species allele frequency at a given SNP locus for the strong‐dispersal, moderate‐dispersal, and weak‐dispersal pools. As the strong‐dispersal pool contained one Cm individual, the expected frequency of Cm alleles at a locus that was hybrid in “Clapper” was 3/8 rather than 1/8; therefore, the array of potential alleles contained 38 “1” values and 62 “0” values. The moderate‐dispersal pool only contained BC3s, so 1/8 was the expected fraction of Cm alleles. As the weak‐dispersal pool contained one Cd individual, the expected fraction of Cm alleles was slightly lower (1/10). To simulate the process of pooled DNA assembly, random draws were made from this distribution up to a simulated read depth of 8, and the number of Cm alleles in the sample was tallied. This process was repeated 1,000,000 times for each pool to create null distributions for Cm allele frequencies in BC3 pooled genomes at “Clapper” hybrid loci and at an assembly depth approximately equal to our actual assemblies.

Subsequently, a p‐value was assigned for each informative SNP in a pool, based on the percent of simulated SNP genotypes that had a count of Cm alleles greater than or equal to the observed number of Cm alleles at that SNP locus. If this percentile‐based p‐value was lower than 0.05, the null hypothesis that Cm alleles were randomly distributed at the locus in a pool was rejected; low p‐values were interpreted as evidence that more Cm alleles were present in the pool than expected by chance alone. For each predicted gene in the genome that contained informative SNPs, an average p‐value was computed using all informative SNPs within the predicted gene sequence. Predicted genes where the null hypothesis was rejected in the strong‐dispersal, but not in the moderate‐ or weak‐dispersal pools were included as potential candidates for influencing seed dispersal. p‐Values assigned to loci using this method were also used to validate candidates identified by the HE heuristic.

2.9Validation of predicted genes

To validate predicted genes from the whole‐genome analysis, cDNA data for a number of species in the order Fagales were aligned to predicted proteins from the Castanea mollissima genome (Carlson et al., 2014; Staton et al., 2014; Appendix S1). cDNA contig consensus sequences were aligned to a database of predicted Castanea protein sequences using the Diamond sequence aligner (Buchfink et al. 2015). A predicted gene was counted as having transcript support if at least one cDNA contig had the predicted gene's protein sequence as its best alignment. The Arabidopsis best hits for each predicted chestnut peptide in genome regions determined to be associated with seed dispersal by the HE method were submitted to gene ontology analysis using g:prolifer (Reimand et al. 2016).

3.1Seed phenotypes

Dispersal trials were conducted for 13 BC3, one American, and one Chinese chestnut in 2014; 11 BC3, one American, and one Chinese chestnut in 2015; and 12 BC3 in 2016 (Table 1). The average mass (mean ± SD) of BC3 seed over the three years was 3.51 ± 1.47 g, ranging between 1.12 and 7.78 g. The average for American chestnut was 3.05 ± 0.17 g; for Chinese chestnut, the average was 7.82 ± 1.01 g. Average seed length of BC3 in 2014 and 2015 was 22.12 ± 2.43 mm, with a range of 17.73–25.86 mm; length for American chestnut was 20.49 ± 0.16, and for Chinese chestnut average length was 24.27 ± 0.69. Width across the wider axis of the attachment‐scar end of the nut was 20.63 ± 3.88 mm for BC3, ranging between 14.42 and 26.4 mm; for American chestnut the mean was 20.04 ± 0.70 mm and for Chinese chestnut 27.29 ± 0.46 mm. In 2015, Cd seeds lost more of their mass through drying (15.85%) than Cm (10.26%) over 2.5 months of cold‐room storage (October 28 to January 14). The individual half‐sib seed lots with the highest rate of caching (68%, 55%, and 25% of seeds recovered in caches rather than recovered eaten) lost moisture at rates similar to Chinese chestnut (8.99%, 8.55%, and 11.61% of moisture lost, respectively) while seeds that were less likely to be dispersed and cached had highly variable (5.41%–31.64%) loss of mass due to drying and an average rate of moisture loss (17.11%) closer to American than to Chinese chestnut.

3.2Seed dispersal

Average recovery rate (% of tagged seeds recovered after 4–5 days) was 66% in 2014, 42% in 2015, and 46% in 2016. Of seeds that were recovered, in 2014, 36.5% were eaten without being moved away from the feeding site, 49.3% were moved and eaten, and 20.5% were moved and cached. In 2015, these numbers were 40.8%, 36.1%, and 23%, respectively; in 2016, they were 66.1%, 29.0%, and 4.8% cached, respectively. In 2016, the apparent shift in proportions was driven by a low caching rate at the Indiana site rather than the addition of the Illinois site. Average dispersal distance for individual BC3s with more than one dispersal event ranged from 4.92 to 9.08 m, which was less than the average for Chinese chestnut (10.49 m) and greater than the single American chestnut that was cached (1.85 m). Chinese chestnuts were dispersed farther and cached more frequently (29.6% of the time) than American chestnuts (3.2% of the time) and most BC3s (averaged over 14 families with at least one dispersed seed: 16.5%). Pools created for genotyping reflected the wide range of variation in dispersal. The average dispersal distance for seeds of trees placed in the strong‐dispersal pool (seven BC3 and one Cm) was 9.08 ± 4.12 m and average caching frequency (no. cached / no. found) 25.3% (Table 1). For the moderate‐dispersal pool (seven BC3), the average distance dispersed was 4.12 m ± 0.78 m, with a caching frequency of 8.2%. In the weak‐dispersal pool (nine BC3 and one Cd), the average distance dispersed was 0.185 ± 0.59 m, and the caching frequency was 0.5%. Mean individual seed size was a statistically significant predictor of mean individual caching distance in a simple linear regression where individuals with average seed dispersal distance 0 (i.e., seeds that were only recovered eaten at the feeding platform) were excluded (t_1,25 = 4.43, p = 0.0002, adjusted r² = 0.42) and seeds cached/total number of seeds recovered (t_1,25 = 2.26, p = 0.03, r² = 0.14) (Figure 1). In a binomial regression, mean seed mass was not a significant predictor of whether an individual mother tree had at least one seed recovered in a cache (z value = 1.394, p = 0.163).

Figure 1

Scatterplot with simple linear regression line of average distance to caching (m) over average seed mass (g) for 25BC3, 2 Castanea dentata, and 2 C. mollissima mother trees measured 2014–2016

3.3Pooled genome SNP genotyping

Enough 100‐bp paired‐end reads (57–67 million) were obtained for each pool to cover the ~800 Mb chestnut genome between 7.2 and 8.5 times, so that each individual tree in each pool was represented by about one read at any locus in the genome. A small fraction of total bases (~2%) were removed from each sample by Trimmomatic due to low read quality prior to analysis. In the strong‐dispersal pool, 341363 informative SNPs with coverage >8 were identified; 177,884 were identified in the moderate‐dispersal pool, and 215,590 were identified in the weak‐dispersal pool. As expected, “Clapper” had a Cd/Cm genotype at 50% of the loci with one allele fixed in Cm and another in Cd, and a Cd/Cd genotype at the other 50%.

3.4Analysis of hybrid regions among pools

The mean value of the heuristic hybridity estimator (HE) over all SNPs in predicted genes with coverage ≥8 was highest for the strong‐dispersal pool (0.44 ± 0.123) (Figure 2) and lowest for the weak‐dispersal (0.294 ± 0.164) (Figure 2). For the moderate‐dispersal pool, the mean value of HE was 0.313 ± 0.174 over all predicted genes (Figure 2). When windows of 10 genes were used, mean difference in HE among windows was greatest between the high‐ and weak‐dispersal pool (0.155 ± 0.088), but the difference between the strong‐ and moderate‐dispersal pools was similar (0.137 ± 0.101) and both were much larger than the average difference between weak‐ and moderate‐dispersal pools (0.019 ± 0.085). Of 2,714 bins of ten predicted genes, there was one region where the difference in HE between the strong‐dispersal pool and the weak‐dispersal pool was >3 standard deviations greater than the mean difference, and 53 bins >2 standard deviations above the mean. There were two bins for which the difference in HE between the strong‐dispersal pool and the moderate‐dispersal pool was >3 SD above the mean, and 58 where the difference was >2 SD above the mean (Table 2).

Figure 2

Histograms of heterozygosity estimates forSNPs in the pooled genome sequences of (a) sevenBC3 and one Chinese chestnut with more frequent and longer‐distance nut dispersal (b) sevenBC3 chestnut with intermediate dispersal distance and low caching frequency and (c) eightBC3 and one American chestnut with low caching frequency and short dispersal distance

3.5Annotations of genes within hybrid regions

Candidate genes for differences in seed dispersal were analyzed for 18 bins with the largest deviations from the mean difference in the heterozygosity estimate between the strong‐dispersal pool and the moderate‐ and weak‐dispersal pools. Fourteen of these bins had a large difference in heterozygosity between the strong‐dispersal pool (pool A) and the others (pools B and C); three were identified based on the difference between the strong‐ and moderate‐dispersal pool; and one was identified based on the difference between strong‐ and weak‐dispersal pools while the strong and moderate‐dispersal pools showed no difference (Table 2). Of these 14 genome regions, seven that had additional support from the simulation‐based estimation of significance were chosen as the most likely candidates (Table 2). Additional individual candidate genes (rather than regions) were identified based on simulation‐based evidence of Cm alleles in the strong‐dispersal pool (Table 3). Examining annotations of predicted genes in these regions revealed several that have plausible roles in seed development and subsequent seed handling and dispersal by squirrels (Figure 3). Many of the predicted genes in these regions aligned to cDNA sequences from chestnuts and other nut‐bearing species in the order Fagales (Table 5). Gene ontology terms that were enriched in candidate regions, as determined by the difference in HE between strong‐dispersal and weak‐dispersal pools, included “regulation of cellular localization” (p = 0.0181), “membrane‐bounded organelle” (p = 0.0428), and “plasma membrane” (p = 0.0181).

Figure 3

Depiction of the role of candidate genes in Castanea nut development (upper right) and the perceptions of nut phenotypes by squirrels (Sciurus) that are hypothesized to cause someBC3 chestnuts to be dispersed farther and cached more frequently than others

4.1Dispersal trials

Previous studies have indicated that seed‐caching rodents are more likely to cache (rather than eat) relatively large seeds, and carry larger seeds farther before caching (Jansen et al., 2002; Moore et al., 2007; Tamura & Hayashi, 2008; Xiao et al., 2005). The dormancy status of seeds is also likely a factor in caching decisions (Moore et al., 2007; Smallwood et al., 2001; Xiao, Gao, Jiang, & Zhang, 2009; Xiao, Gao, Steele, & Zhang, 2009): squirrels are more likely to eat seeds perceived as nearing germination and to cache seeds perceived as dormant. The value of squirrel caches to nut survival (fitness) has been amply demonstrated (Lichti, Steele, & Swihart, 2017), so seed phenotypes that influence squirrel caching have evolutionary significance for nut‐bearing trees. In our study, seed size, as measured by seed mass, was associated with both dispersal distance (Figure 1) and likelihood of caching vs. consumption. Seed mass was highly variable among the backcrossed trees in the study, which we deliberately chose based on seed size (Table 1) from a much larger population of BC3 chestnuts. Seeds with greater mass tended to exhibit greater dispersal distance and greater likelihood of caching vs. immediate consumption by squirrels. Chinese chestnuts in this study were much larger than American chestnuts (Table 1), so BC3 hybrids with larger seeds (similar to Chinese chestnut) were more likely to be cached. Variation in seed dispersal distance and caching likelihood among BC3 trees, however, was not fully explained by seed size (Figure 1). There is no documented difference in seed dormancy between American and Chinese chestnut—both species must undergo a dormant period of several months to germinate (Saielli et al. 2012), and in both species, the seed is metabolically active during its dormant phase because chestnuts are recalcitrant seeds (Leprince, Buitnik, & Hoekstra, 1999; Roach et al., 2009). The sugar content of chestnuts under cold storage increases while starches diminish (Ertan, Erdal, Gulsum, & Algul, 2015). Differences in genes that regulate dormancy or signals that communicate dormancy to squirrels, or binding sites for regulatory molecules, are therefore less likely to be false‐positives than structural or housekeeping genes.

4.2Genome scan for loci involved in caching decisions

Interspecific hybrid phenotypes are not always intermediate between parents (Woeste et al. 1998), and intermediate genotypes do not simply correspond to intermediate phenotypes. In hybrid chestnuts, for example, gene interactions can cause unpredictable dormancy‐related phenotypes: hybrid seeds derived from a cross with Allegheny chinkapin (Castanea pumila) as pollen parent and Cm as the seed parent exhibited reduced seed dormancy (Jaynes 1963; Metaxas 2013). No precocious germination of seeds was observed in the course of our experiment, but subtle phenotypes may have been present. In red oaks (Quercus section Lobatae), which are closely related to chestnuts, dormancy and germination appear to be controlled primarily by the pericarp, the dry fruit structure that makes up the hard outer “shell” of both oak acorns and chestnuts (Smallwood et al., 2001). There is evidence that squirrels make use of changes in the pericarp to sense impending germination in oaks and chestnuts. These changes include degradation of pericarp waxes and the release of low molecular weight volatile compounds from inside the pericarp (Paulsen et al., 2014; Sundaram et al. 2015, Sundaram, 2016). By demonstrating that a germinating white oak embryo inside a “dormant” red oak shell is perceived as dormant by squirrels, Steele et al. (2001) showed that signals at the pericarp surface may be more important than signals from the kernel. If some chestnut hybrids have a thicker pericarp wax layer than American chestnut, squirrels might perceive the seeds as reliably dormant and cache them more frequently.

By sequencing pools of chestnuts with different seed dispersal phenotypes, we attempted to identify regions of the genomes of BC3 hybrids where Chinese chestnut allele frequencies were higher than expected in the most frequently dispersed trees. The HE statistic seems to have captured the elevated heterozygosity that is characteristic of Cm/Cd hybrid gene loci. The inference space for our study was limited to the phenotypic effects, in a BC3 population, of loci where “Clapper” may have contributed a Cm allele, which only included half of the genome. Our ability to find genomic regions that were hybrid in one pool but not the others was impeded by the uncertainty associated with estimating heterozygosity in pooled sequence data. By comparing results from pooled data with individual chestnut genome sequences, we determined that the individual “Clapper” genome was significantly more heterozygous than individual Cd genomes at SNP loci in most of the seed dispersal candidate regions (Table 4), which indicated that these regions were plausible candidates. Differences in Cm allele frequencies among pools (Figure 2) were accounted for by the simulation‐based method for identifying outliers; only a small fraction of the many genes with divergent allele frequencies between Cd and Cm were included in candidate regions. The study design limited our inference to maternal effects, but as squirrel caching decisions appear to be influenced strongly by characteristics of the maternally derived pericarp (Steele et al. 2001), the paternal contribution to differences in dispersal is likely to be small. Finally, the small number of genotypes utilized (three pools derived from 24 individual trees) limits the strength of conclusions drawn from this study because the number of false‐positive candidate loci is inversely related to sample size. The importance of candidate genes that remained after statistical validation was rendered more plausible, however, because their predicted function often corresponded to factors known to influence squirrel caching decisions—seed size and the perception of dormancy. These remaining candidates point to new hypotheses on seed/seed‐disperser coevolution in hardwood trees. Seed size, the most obvious dispersal‐associated phenotype that distinguishes Cd and Cm, is a likely a complex trait in chestnut as it is in other plants (e.g., Gnan, Priest, & Kover, 2014). Several candidate loci identified in this study had annotations that point to a potential role in seed development and seed size. The EMBRYONIC FLOWER 2‐like (EMF) gene on LGC (c.g6050; Table 2) could directly influence seed size by regulation of development of female flower parts. EMF2 in Arabidopsis encodes a Polycomb group protein (Yoshida et al., 2001) that regulates vegetative growth and development by suppressing the flower‐development program, Yoshida et al., 2001). The predicted EMF2 gene in chestnut had strong transcript support from C. mollissima and C. dentata and was one of 5 EMF2‐like genes predicted (by AUGUSTUS) in the entire chestnut genome (Carlson et al., 2014). A candidate gene on LGL (g6184) was similar to GIF2 of Arabidopsis, which regulates the expansion of cotyledons (Kim & Kende, 2004). Cotyledons make up the majority of the mass of a chestnut seed. Several other candidates, including a LATERAL ROOT PRIMORDIUM 1 (Kuusk, Sohlberg, Magnus, & Sundberg, 2006) homolog on LGH, a VERNALIZATION‐3‐like gene (VIL1) on LGA from (Table 2), a gene similar to FRIGIDA‐like 4 of Arabidopsis on LGK, also may function in the regulation of flower development. The latter two loci are both involved in the FLOWERING LOCUS C (FLC) regulatory pathway in Arabidopsis (Greb et al., 2007; Michaels, Bezerra, & Amasino, 2004). Whether homologs of these flowering regulatory loci have effects on the development of floral parts, and thereby influence seed size in chestnut, is uncertain. It is also possible that they affect dispersal by modifying seed dormancy, given that FLC and its interactors have documented pleiotropic effects on the regulation of seed dormancy and germination (Chiang et al. 2009). If such pleiotropic effects exist, changes in seed size due to natural selection could also lead to changes in dormancy.

4.3Genomic loci associated with differences in seed dispersal: pericarp‐mediated dormancy

In red oak acorns, which are anatomically similar to chestnuts, the pericarp prevents absorption of water by the embryo and allows germination only after the pericarp's permeability has increased following a period of cold storage (Peterson 1983; Steele et al., 2001). The pericarp is derived from the ovary walls of female chestnut flowers and consists of several layers of lignified cells with a waxy coating on the outermost layer. Both the breakdown of this waxy layer and the subsequent release of volatile compounds from the pericarp may serve as olfactory cues for squirrels that a seed is approaching germination and is therefore more perishable and a poor candidate for caching. The candidate genes we identified include some that may be involved in the formation of pericarp cells, some that may influence the composition of waxes on the pericarp surface, and others that may influence the release of volatile compounds from the pericarp. While none of the transcriptomic data from Fagales trees we aligned to the Cm predicted gene set was seed specific, several dispersal candidate genes were supported by cDNA contigs from several species (Table 5).

Our analysis identified several predicted genes that may have a role in cell wall modification during nut development and ripening in chestnut. These include an extensin‐like predicted gene on LGA (LGA.sd03; g10648), which has a similar protein sequence to an LRX‐family gene in Arabidopsis (Baumberger et al., 2003); the LRX family has been implicated in the modification of plant cell walls (Draeger et al., 2015). Predicted genes similar to pectin methylesterases were identified on LGA (g11556) and LGL (g7556). The latter predicted gene (LGL.Sd17) was similar to a gene in Arabidopsis expressed in developing siliques (seed pods) (Louvet et al., 2006) and could be involved in the formation of the lignified chestnut pericarp. Pecinesterase genes are active in maturing (lignifying) wood in poplar (Mellerowicz, Baucher, Sundberg, & Boerjan, 2001). Conversely, these genes could be directly involved in the germination process: in yellow cedar (Chamaecyparis nootkaensis), pectinesterases are active in germinating seeds (Ren & Kermode, 2000), and Arabidopsis lines with overexpression of a PME inhibitor showed more rapid germination (Müller et al. 2005). Given the importance of the pericarp in squirrel perceptions of seed perishability, these cell‐wall modification genes may influence seed dispersal by acting in the maturing pericarp, rather than in the germinating embryo. Other candidate genes that could influence formation of the pericarp include a cellulose synthase (CESA2)—like gene (Beeckman et al., 2002; Li et al. 2009). Adjacent to the pectinesterase at SD17 on LGL and a gene similar to a sugar transport carrier in castor bean (Ricinus communis) (STC_RICCO; LGE_g7467). Several candidate genes we identified appear to have a role in lipid metabolism that may be related to the formation and/or degradation of pericarp wax layers (Pollard, Beisson, Li, & Ohlrogge, 2008). The potential importance of fatty‐acid metabolic processes in regulating squirrel dispersal was explored by Sundaram (2016), who showed that differences in the outer wax layer of the pericarp influence squirrel perception of seed dormancy. Nonspecific lipid transfer proteins (LGC.Sd15) in Arabidopsis are involved in the formation of suberin in crown galls (Deeken et al., 2016), various tissues of tomato in response to drought stress (Trevino & O'Connell, 1998), and the surface wax of broccoli leaves (Pyee, Yu, & Kolattukudy, 1994). Another lipid‐modifying gene, a cytochrome p450 oxidase (g3304), occurs at the same locus as the NLTL‐like predicted gene on LGC, and a second cytochrome p450 oxidase (LGA_g3246) was identified as a candidate gene on LGA. Based on alignments to the UniProt database, these genes are both similar to a gene from tare (Vicia sativa) that catalyzes hydroxylation of myristic acid and other fatty acids (Le Bouquin, Pinot, Benveniste, Salaun, & Durst, 1999).

Squirrels perceive volatile compounds from seeds as cues of metabolic activity and impending germination (Sundaram, 2016). Volatile compounds are thought to escape the pericarp as it becomes more porous and germination approaches. One particularly interesting locus appeared to contain a cluster of four volatile terpene synthase genes, which are most similar to terpene synthesis genes highly expressed in the fruits of strawberry (Aharoni et al., 2004) that are thought to influence the fruit's flavor and aroma profile. Nerolidol is a sesquiterpene compound found in many plants (Chan, Tan, Chan, Lee, & Goh, 2016). Sundaram (2016) found that release of beta‐amyrin, a triterpene, was associated with germination of chestnuts. While these compounds are distantly related, their synthesis may be metabolically linked by production of the intermediate squalene. Nerolidol synthase converts farnesyl diphosphate (FPP) to nerolidol. In a yeast study, overexpressing FPP synthase and squalene synthase greatly increased beta‐amyrin production (Zhang et al., 2015). Beta‐amyrin has been associated with wax degradation in other plants (Buschhaus & Jetter, 2012), so it could degrade the cuticular waxes of the outer pericarp as it is released, preparing the seed for germination (Sundaram, 2016). The genes found here do not directly influence beta‐amyrin production, but could influence production of substrate molecules or divert carbon away from beta‐amyrin production. Interestingly, of the three nerolidol synthase‐like genes at this locus, only one showed evidence of expression in Chinese chestnut and two others showed evidence of expression in American chestnut and oaks, but not Chinese chestnut (Table 4), possibly indicating interspecific differences in the expression of these genes. There was little evidence of their expression in the non animal‐dispersed taxa examined (Alnus, Betula) nor in Fagus. If the expression of multiple copies of nerolidol synthase in American chestnut leads to an increase in the activity of volatile organic compounds that degrade pericarp waxes, the result in BC3 seeds that express the Cd alleles could be a signal to squirrels to eat rather than cache these seeds. Our results support previous studies (e.g., Smallwood et al., 2001) that point to seed dormancy and germination as a primary influence on squirrel caching decisions. In particular, they support the notion that pericarp waxes and volatile compounds are important for conditioning squirrels’ perceptions of seed dormancy in nut‐bearing trees. The potential evolutionary role of loci with pleiotropic effects on flower development, seed size, and germination in nut‐bearing trees merits further study. The possibility that lipid‐ and secondary metabolite‐synthesis genes expressed in developing pericarp tissues are ultimately important for seed dispersal phenotypes should be investigated in squirrel‐dispersed tree species and interspecific hybrids. The results of the present study need to be validated and clarified using larger numbers of plants and individual, rather than pooled, genotypes. We hypothesize that genes controlling differences in seed dispersal are primarily expressed during flower and seed development and maturation (the formation of the pericarp) rather than during dormancy, when dispersal takes place.