Moths passing in the night: Phenological and genomic divergences within a forest pest complex

Abstract Temporal separation of reproductive timing can contribute to species diversification both through allochronic speciation and later reinforcement of species boundaries. Such phenological differences are an enigmatic component of evolutionary divergence between two major forest defoliator species of the spruce budworm complex: Choristoneura fumiferana and C. occidentalis. While these species interbreed freely in laboratory settings, natural hybridization rates have not been reliably quantified due to their indistinguishable morphology. To assess whether temporal isolation is contributing to reproductive isolation, we collected adult individuals throughout their expected zone of sympatry in western Canada at 10‐day intervals over two successive years, assigning taxonomic identities using thousands of single nucleotide polymorphisms. We found unexpectedly broad sympatry between C. fumiferana and C. occidentalis biennis and substantial overlap of regional flight periods. However, flight period divergence was much more apparent on a location‐by‐location basis, highlighting the importance of considering spatial scale in these analyses. Phenological comparisons were further complicated by the biennial life cycle of C. o. biennis, the main subspecies of C. occidentalis in the region, and the occasional occurrence of the annually breeding subspecies C. o. occidentalis. Nonetheless, we demonstrate that biennialism is not a likely contributor to reproductive isolation within the species complex. Overall, interspecific F1 hybrids comprised 2.9% of sequenced individuals, confirming the genomic distinctiveness of C. fumiferana and C. occidentalis, while also showing incomplete reproductive isolation of lineages. Finally, we used F ST‐based outlier and genotype–environment association analyses to identify several genomic regions under putative divergent selection. These regions were disproportionately located on the Z linkage region of C. fumiferana, and contained genes, particularly antifreeze proteins, that are likely to be associated with overwintering success and diapause. In addition to temporal isolation, we conclude that other mechanisms, including ecologically mediated selection, are contributing to evolutionary divergence within the spruce budworm species complex.


| INTRODUC TI ON
Temporal isolation of potentially interbreeding evolutionary lineages can be a significant contributor to species diversification (Hendry & Day, 2005;Nosil, 2012;Rundle & Nosil, 2005). Such phenological divergences may occur at different temporal scales, whether within a day (e.g., Devries et al., 2008;Schöfl et al., 2009;Ueno et al., 2006), seasonally (e.g., Adamski et al., 2018;Bell et al., 2017;Yamamoto & Sota, 2009), or between years (Cooley et al., 2003;Gradish et al., 2015). Multiple ecological and physiological processes may lead to temporal isolation, and it is almost always difficult to determine whether present-day phenological divergence was the initial cause or a later reinforcing consequence of reproductive isolation (Egan et al., 2015;Taylor & Friesen, 2017). In any case, the essential first step to inferring whether phenological divergences are contributing to contemporary reproductive isolation among taxa is to rigorously document the extent and nature of their temporal separation and associated characteristics.

Choristoneura fumiferana ranges across the northeastern United
States and much of Canada's boreal forest, feeding mainly on white as C. o. biennis Freeman . While C. o. occidentalis and C. o. biennis are parapatric in British Columbia and Alberta, Canada, each shares a substantial zone of sympatry with C. fumiferana in western Alberta (Dupuis et al., 2017;Lumley & Sperling, 2011a, 2011b. Intriguingly, pairings of C. fumiferana and C. occidentalis hybridize successfully in laboratory settings (Harvey, 1997;Nealis, 2005;Sanders et al., 1977), but few genetically admixed individuals indicative of hybridization have been observed in nature Lumley & Sperling, 2011a). This suggests that C. fumiferana and C. occidentalis are reproductively isolated, but specific mechanisms that underlie their isolation are not well understood.
Seasonal flight peaks of C. fumiferana and C. occidentalis have been inferred to differ within their zone of sympatry in Alberta (Lumley & Sperling, 2011a). Assuming flight peaks are adequate proxies of reproductive timing for each species, divergence of flight peaks suggests that temporal isolation may be contributing to their continued reinforcement as distinct species .
However, although flight peaks may differ between C. fumiferana and C. occidentalis, any degree of flight period overlap would suggest that gene flow among these sympatric species is possible, barring other pre-or post-zygotic reproductive incompatibilities. Therefore, explicit quantification of flight period overlap is needed and requires surveys repeated at regular intervals spanning the zone of C. fumiferana/C. occidentalis sympatry. Another important consideration is that assigning species identifications to Choristoneura individuals is generally unreliable without genomic data (French et al., 2020), weakening inferences of phenological divergence based on historical data, morphology, sex pheromones, or life-history traits (e.g., Powell & De Benedictus, 1995;Smith, 1954;Volney et al., 1984). Targeted phenological studies using genomic data are therefore required to resolve mechanisms that contribute reproductive isolation within this species complex.
In addition to seasonal divergences in the flight periods of C. fumiferana and C. occidentalis, it has also been hypothesized that biennialism contributes to their reproductive isolation (Shepherd et al., 1995). While C. fumiferana and C. o. occidentalis are univoltine, larvae of C. o. biennis undergo an obligate second diapause and overwinter a second time before completing a two-year (semivoltine) life cycle (Harvey, 1967;Nealis, 2005). This biennialism results in adult C. o. biennis flying in high abundance every second year, separated by a year of very few adults (Nealis & Turnquist, 2003;Shepherd et al., 1995;Zhang & Alfaro, 2002). However, biennialism on its own cannot account for complete reproductive isolation among sympatric taxa. If sex ratios, abundances, and flight periods of two sympatric taxa are very similar and they do not exhibit prezygotic barriers to reproduction, interspecific mating will not occur in years in which only one taxon reproduces but will be substantial for years in which both taxa reproduce. In other words, reproductive isolation would not be maintained in years in which both semivoltine C. o. biennis and univoltine C. fumiferana emerge as adults. Nonetheless, variation in voltinism may still indirectly contribute to within-year reproductive isolation if C. o. biennis flies earlier than C. fumiferana due to a phenological "head start" in years when both emerge as adults. This hypothesis would be supported if C. o. biennis populations consistently fly before sympatric C. fumiferana populations, but must be tested at fine temporal and spatial grains using genomic data to consistently assign taxon identities to individuals.
In this study, our principal aims were to (1)  of sympatry and surrounding areas over two successive years and genotyped using next-generation sequencing to assign taxon identities to sampled individuals. We also identified genomic regions under putative divergent and ecologically mediated selection, using a combination of F ST outlier and genotype-environment association analyses, to identify additional factors that may be contributing to genomic variation and reproductive isolation among these taxa.

| Specimen collection
We sampled adult Choristoneura at 10 locations in Alberta and British Columbia at 10-day intervals during their 2017 and 2018 flight seasons ( Figure 1, Table S1). Six collection locations in Alberta (d-i, Figure 1) were within the documented zone of sympatry of C. fumiferana and C. o. biennis . Three locations were chosen in British Columbia (a-c) outside of the known range of C. fumiferana (Dupuis et al., 2017) based on aerial surveys depicting severe C. o. biennis defoliation in 2015 and 2016 (Westfall & Ebata, 2016. The 10 th location (Elk Island, Alberta) was geographically distant from the documented range of C. o. biennis and therefore expected to contain only C. fumiferana Lumley & Sperling, 2011b).
Adult moths were collected during 13 June-25 August 2017 and 1 June-3 September 2018 (Tables S2 and S3). At each of the 10 locations, four green Unitraps (Contech, Delta, British Columbia, Canada) were hung two meters above ground level on mature white spruce, Engelmann spruce, or Douglas fir, with at least 40 meters between each trap to reduce interference (Allen et al., 1986;Sanders, 1996).
Three Unitraps were baited with rubber septum lures containing 100 μg of synthetic 95:5 E, Fredericton,New Brunswick,Canada), the dominant component of the C. fumiferana sex pheromone (Silk & Kuenen, 1988 (Sanders et al., 1977;e.g., Lumley & Sperling, 2011a, 2011bBrunet et al., 2017), although the optimum E/Z ratio for C. o. biennis is unknown (Sanders, 1971;Sanders et al., 1974). The fourth trap did not contain a lure, acting as a control. We included one strip of 10% dichlorvos in each trap as a killing agent (Hercon Environmental, Emigsville, Pennsylvania, USA). Trap catch was collected every 10 days and specimens were then stored at −20°C. When we began our study in 2017, aerial reports from 2015 and 2016 did not resolve whether even-or odd-year defoliation was dominant, so we could not be sure which year C. o. biennis was likely to fly in high abundances (Westfall & Ebata, 2016. We therefore contin-

| DNA extraction and ddRADseq
A minimum of three adult individuals were randomly subsampled for genotyping within each collection period at each collection location when available. Up to 10 individuals were genotyped from collection intervals if samples were numerous. To maximize coverage of total flight period, we also selected singletons that were trapped earliest or latest in each collection season. Whenever F I G U R E 1 Collection locations (a) and phenology (b) of adult Choristoneura across the study area. Moths were pheromone-trapped at each location during ~10 days before collection. Histograms show numbers of adult Choristoneura trapped in each collection interval, with traps being emptied on: 1 = 11-15 June, 2 = 21-27 June, 3 = 1-7 July, 4 = 11-17 July, 5 = 21-27 July, 6 = 31 July -6 August, 7 = 10-17 August, 8 = 20-24 August, 9 = 30 August -2 September (Tables S2 and S3 show days since 1 January for each location). Collection intervals without number digits from 2017 indicate traps were not deployed, while numbered collection intervals without histogram bars indicate no moths were trapped. Inset map of North America shows rectangular study area in Alberta and British Columbia with C. fumiferana (green) and C. occidentalis (blue) ranges (adapted from Dupuis et al., 2017). Map was prepared in QGIS version 3.6.0 using base maps obtained from https://github.com/stame n/terra in-classic and https://github.com/Carto DB/basem ap-styles possible, we selected both fresh and worn moths to cover the breadth of variability present. Individuals that were wet upon collection were not selected for genotyping due to the possibility of degraded genomic DNA.
DNA was extracted using DNeasy Blood & Tissue DNA Kits (QIAGEN, Hilden, Germany), with the addition of bovine pancreatic ribonuclease A treatment (RNaseA, 4 μl at 100 mg/ml; Sigma-Aldrich Canada Co., Canada). DNA of each specimen was ethanol precipitated and resuspended in millipore water for further con-

| DNA data processing
Single-end, 75 bp Illumina reads were demultiplexed using the process_radtags program in Stacks 2 version 2.3e (Rochette et al., 2019). We discarded any reads with Phred scores <20 over 15% of the read length and any reads that failed the Illumina purity filter.
Retained reads were trimmed from 75 bp to 67 bp after removing the 8 bp index sequence. Some sequencing error was identified in the PstI site, so another 5 bp were removed from the 5' end with Cutadapt version 1.18 (Martin, 2011), resulting in a final read length of 62 bp.
We aligned these 62 bp reads to a 30,339 scaffold C. fumiferana draft genome, bw6 (Larroque et al., 2019), using Burrows-Wheeler Aligner version 0.7.17-r1188 (Li & Durbin, 2009). Resulting SAM files were converted to BAM format for use in the Stacks 2 pipeline with the ref_map program stipulating a single population containing all individuals (Rochette et al., 2019). Stacks-filtered single nucleotide polymorphisms (SNPs) were further filtered in VCFtools version 0.1.16, using missing data across SNPs and minor allele frequency thresholds of 95% and 0.10, respectively (Danecek et al., 2011). We also thinned SNPs that were within 10k bp of one another, to minimize physical linkage (MacDonald et al., 2020).

| Taxonomic identifications and hybrid individuals
Taxon identities were determined with a combination of principal component analysis (PCA), using the R package adegenet version 2.1.1 (Jombart & Ahmed, 2011), and model-based clustering analyses, using the programs structure version 2.3.4 (Falush et al., 2003;Pritchard et al., 2000) and NewHybrids version 1.1 beta3 (Anderson & Thompson, 2002). For structure analyses, we ran 10 iterations for each of k = 1-10 with a burn-in period of 100,000 followed by 1,000,000 Markov chain Monte Carlo repetitions. The locpriors parameter, with values corresponding to the 10 collection locations, was used to better resolve spatial genetic structure without biasing the model (Porras-Hurtado et al., 2013). We used StructureSelector (Li & Liu, 2018) to determine the number of genetic clusters with the highest support using both LnP(k) (Pritchard et al., 2000) and Δk (Evanno et al., 2005). CLUMPAK version 1.1 (Kopelman et al., 2015) was used to average the 10 replicates of each k and create matrices of admixture coefficients (Jakobsson & Rosenberg, 2007). For NewHybrids analysis, we averaged the results of 5 iterations of the k = 2 dataset. Each iteration ran with a burn-in period of 100,000 followed by 1,000,000 sweeps.
Because NewHybrids experiences underflow errors with large datasets (Elliot & Russello, 2018;Hu et al., 2021), we chose the first 400 SNPs in the dataset for easier reproducibility (Gramlich et al., 2018). We also ran a second analysis using 400 randomly chosen SNPs to check if results were congruent (Gramlich et al., 2018). In both NewHybrids analyses, we assigned individuals to the two parental classes (parental 1, parental 2) using >99% membership to the two genetic clusters found by structure (Bell et al., 2015).
either genetic cluster based on structure admixture coefficients for k = 2. We used NewHybrids to assign hybrid class membership to each individual (parental 1, parental 2, F 1 , F 2 , backcross to parental 1, and backcross to parental 2), executed in EasyParallel version 1.0 (Zhao et al., 2020). We then reran structure to resolve substructure within each of the two species (primary genetic clusters). F 1 hybrids, as determined by NewHybrids, were excluded in the second runs.
VCFtools was used to subset the genomic dataset for each species and refilter SNPs in accordance with the filtering parameters detailed above. Relationship between C. o. biennis and C. o. occidentalis was inferred based on admixture coefficients >50% Brunet et al., 2017).

| Evaluation of temporal isolation
We estimated temporal overlap of C. fumiferana and C. occidentalis flight periods using the taxonomic identities of sequenced individuals and two co-occurrence indices. The first index was calculated as the number of collection periods in which both C. fumiferana and C. occidentalis were present divided by the total collection periods that either species was present. This was first completed on a location-by-location basis for all sympatric collection locations. Resulting values were averaged across these locations to give an overall measure of temporal overlap. We also calculated this index by first pooling all individuals from all sympatric collection locations and then dividing the number of collection periods in which both species were present by the total collection periods that either species was present. Comparison of index values calculated on a location-by-location vs. pooled-location basis allowed us to infer whether combining species occurrence records across a broad spatial extent overestimates the actual temporal overlap of sympatric populations. We also used a second index using species' abundances. This index was calculated as the total number of C. fumiferana and C. occidentalis individuals that were collected in the same, pooled-location, collection interval divided by the total number sequenced individuals across all collection intervals at all locations. This was again completed on a location-by-location basis.
F 1 hybrids were not included in either index.
To allow more controlled comparisons among locations under varied climatic conditions, we related phenology to local thermal accumulation using degree-day estimations (Zalom et al., 1983).
We recorded degree-day accumulation using HOBO Pendant Temperature/Light Data Loggers (Onset Computer Corporation, Bourne, Massachusetts, USA) during the 2018 collection season, with data loggers hung two meters above ground in the shade at each location. Resulting temperature data were then calibrated against records for nearest weather stations (Table S4) Wilcoxon-Mann-Whitney tests were used to evaluate differences between logger and station temperatures (Table S5). We applied median difference between logger and station temperature as a correction factor applied to all daily high and low temperatures recorded by weather stations between 1 January and trap removal. Régnière (1990) showed that photoperiod has minimal effect on postdiapause emergence of C. fumiferana, so we did not use this index in our analyses.
We used degree-day accumulation to compare phenology between disparate locations. This was calculated using single triangulation (Lindsey & Newman, 1956;Zalom et al., 1983), which is more accurate than alternatives when temperatures are far below the lower developmental threshold of an insect (Roltsch et al., 1999). We used 5.5°C as the lower development threshold for both C. fumiferana and C. o. biennis based on previous work (e.g., Kemp et al., 1986;Shepherd, 1961). We did not use an upper developmental threshold, as temperatures only exceeded the C. fumiferana threshold of 38°C (Régnière et al., 2012) during one-hour-long period at Red Lodge (location i). To summarize results, we plotted the number of genotyped individuals collected each 10-day interval against the number of within-location accumulated degree-days for each location.

| Genomic regions with high divergence
We used BayeScan version 2.1 (Foll & Gaggiotti, 2006) to identify specific genomic regions that were highly diverged among Choristoneura taxa. BayeScan has been recognized as an effective method for identification of F ST outlier loci between discrete genetic clusters that may warrant further investigation using whole-genome sequence data (De Mita et al., 2013;Lotterhos & Whitlock, 2014;Narum & Hess, 2011 Assignment of individuals to either C. fumiferana or C. occidentalis was based on admixture coefficients greater than 0.9 in structure analyses. We assigned C. occidentalis subspecies identities using admixture coefficients greater or less than 0.5. For all Bayescan runs, prior odds were set to 10, thinning interval to 10, number of pilot runs to 20, length of pilot runs to 5000, burn-in length to 50,000, and number of output iterations to 10,000. We assessed significance of F ST outliers using the false discovery rate (FDR) criterion and a q-value threshold of 0.05 (−log10 q-value ~ 1.3) (Benjamini & Hochberg, 1995

| Genotype-environment associations
To determine whether variation in environmental conditions is related to genomic variation within the spruce budworm species complex, we used latent factor mixed modeling (LFMM) to assess correlations between allele frequencies and environmental conditions across sampling locations. It is possible that both C. fumiferana and C. occidentalis share ancestral adaptations to environmental conditions and may therefore exhibit the same genotype-environment associations. Alternatively, these species may be adapted to variation in environmental conditions in different ways and exhibit different genotype-environment associations. We therefore quantified genotype-environment associations for all Choristoneura individuals was again used to control for false positives associated with multiple tests with a q-value threshold of 0.05 (Benjamini & Hochberg, 1995). For loci with multiple significant environmental associations (likely due to spatial correlation of environmental variables), we determined the strongest association based on median |z|-scores (MacDonald et al., 2020;Martins et al., 2018).

| Genomic locations and identities of SNPs
We performed BLASTN searches using flanking sequences around candidate/outlier loci derived by BayeScan and LFMM 2 to determine putative SNP identity. We used the Larroque et al. (2019) draft genome to retrieve 5 kb on either side of each SNP, or all bases between the SNP position and the end of the scaffold if it was within 5 kb. We then searched NCBI's nonredundant sequence databases using BLASTN version 2.11.0 to find nucleotide matches for these regions (Madden, 2002). We inferred the function of the top sequence matches using the UniProt database (The Uniprot Consortium, 2019). Additionally, we determined the location of these SNPs on the C. fumiferana linkage map (Picq et al., 2018) to evaluate whether major genetic differences between and within these taxa are restricted to linkage groups. We filtered raw Illumina-generated reads from 267 Choristoneura specimens in a preliminary run of Stacks 2. We removed five individuals with >30% missing data and reran Stacks 2 with 262 individuals.

| Specimen collections and DNA data processing
Two more individuals were removed after re-running Stacks 2; one with >30% missing data in the re-run, and the other with divergent placement in PCA, which we suspect represented another third spe-  We found greatest support for two genetic clusters among the 260 individuals using both the Δk method (Evanno et al., 2005) (Table 1) based on the proportion of F 1 hybrid individuals relative to those in the two parental genetic clusters (locations d-j, Figure 1). Table S7.

Records of all genotyped individuals are listed in
Hierarchical structure analyses suggested the existence of two indistinct genetic clusters within C. occidentalis (Figures S3-S5). Individuals assigned to the largest cluster (n = 147) were identified as C. o. biennis in accordance with Brunet et al. (2017). The six remaining individuals F I G U R E 2 Principal component analysis, structure barplot, and NewHybrids barplot of SNPs for 260 genotyped Choristoneura. Individuals with <90% assignment to either cluster in the structure analysis are labeled "F 1 hybrids" based on NewHybrids assignment were identified as C. o. occidentalis, which were primarily collected at Clearwater in British Columbia (location a). There was no apparent substructure within C. fumiferana based on either PCA or hierarchical structure analysis ( Figure S6). This was supported by incongruence between the LnP(k) and Δk methods, in addition to low support for any one Δk value (Figures S7 and S8; structure plot not shown).

| Temporal separation between and within species
The median flight dates of genotyped C. fumiferana and C. o. biennis

| Genomic regions with high divergence
BayeScan identified 96 loci (3.4% of 2831 SNPs) with F ST values that were significantly higher than background divergence between C. fumiferana and C. occidentalis (Table S8) (Table S9). Mean F ST of these loci was 0.17 (SD = 0.035), and overall F ST between the subspecies was 0.01 (SD = 0.08).

| Genotype-environment associations
LFMM 2 identified two loci that were significantly associated with mean precipitation of warmest quarter (all Choristoneura, Figure   S9, Table S8) and mean annual precipitation (C. fumiferana independently, Figure S10, Table S10). Neither of these loci were identified as significant F ST outliers by BayeScan (Tables S8 and S10 (Table S9). Of these, 2 were associated with genes encoding antifreeze proteins. One locus flanking sequence had a match on the C. fumiferana linkage map on linkage group 16 (Picq et al., 2018). In the C. fumiferana dataset, 1 locus was identified by LFMM 2 as having a significant environmental association (Table S10). Its flanking sequence had no match on the C. fumiferana linkage map (Picq et al., 2018).

| Temporal separation and isolation between species
Our study demonstrates substantial, but incomplete, phenological divergence between sympatric C. fumiferana and C. occidentalis.
This trend is apparent in either pooled-location or location-bylocation evaluations of temporal relationships among closely related taxa. When pooling individuals from all locations, temporal overlap between species occurred in most (63.6%) collection intervals over two years (Figure 4). In contrast, temporal overlap was observed in far fewer intervals on a location-by-location basis (12.4%; Figure 5). Similar estimates of temporal overlap were found using our second co-occurrence index that accounted for the abundance of each species. Together, these results suggest that combining species occurrence records across a broad spatial extent overestimates the actual temporal overlap of sympatric populations. The extent of temporal contact between these species also varied substantially between years, highlighting the need to consider variation in voltinism in these investigations. In  (Tables S2 and S3). Traps were not deployed for 2017 collection intervals 1 and 9 We observed two instances wherein C. o. biennis was detected in its region of sympatry with C. fumiferana at the beginning of the collection season in late June, much earlier than the taxon's regional flight peak ( Figure 5, Obed Lake and Red Lodge). It is possible that these two individuals eclosed in British Columbia in late June and were blown into Alberta via prevailing winds, as their calendarday phenology matches that of C. o. biennis from Clearwater and Vavenby (see Lumley et al., 2020). Such dispersal may also explain why the frequency of genotyped C. o. biennis was unusually high at Wildhorse Lake and Pembina Forks locations, which are closer to the main range of this species (Figure 1). Choristoneura are known to disperse great distances; Dobesberger et al. (1983) found that C. fumiferana can travel more than 450 km under certain conditions and that these long-distance dispersal events are common during outbreaks of C. fumiferana in eastern North America (Greenbank et al., 1980;Sturtevant et al., 2013). Results from structure reinforce this prevailing-wind hypothesis and suggest that gene flow between these species is asymmetric; more C. fumiferana individuals exhibit some membership to the C. occidentalis genetic cluster than C. occidentalis individuals to the C. fumiferana cluster (Figure 2; and see Brunet et al., 2017).
Our repeated sampling across consecutive years helps clarify temporal and spatial relationships between C. fumiferana and C. occidentalis. Previously, range-wide sympatry of these species has been explored without explicit consideration of their temporal interactions Dupuis et al., 2017;Lumley & Sperling, 2011b), despite the known biennialism of C. o. biennis.
We found that C. o. biennis extends at least to Elk Island, Alberta  (Tables S2 and S3). Collection intervals without number digits from 2017 indicate traps were not deployed, while numbered collection intervals without histogram bars indicate no moths were trapped and environmental conditions (Silva-Brandão et al., 2015). Indeed, adult C. fumiferana collected in Unitraps are not always from the same population as larvae feeding at the same location (James et al., 2015). Likewise, the C. o. occidentalis collected in both years from Clearwater, British Columbia, were outside of their known range but within plausible dispersal distance from suitable habitat.
We find that adult occurrences of both C. o. biennis and C. o. occidentalis span a larger geographic area than previously documented.
However, we cannot be sure whether these are vagrant individuals or representatives of local, reproducing populations. We therefore recommend that future studies should compare larval and adult genotypes in the zone of sympatry between C. fumiferana and C. occidentalis to better understand the reproductive range of each species, and their interaction with host conifers.

| Extent of hybridization between species
Several studies have demonstrated that Choristoneura taxa can hybridize freely in laboratory settings (Harvey, 1997;Nealis, 2005;Sanders et al., 1977). However, our SNP-based identifications of wild specimens allowed us to confirm that hybrid adults are rare in their zone of sympatry, despite some degree of reproductive overlap. Using NewHybrids, we estimated a hybridization rate of 2.9% for C. fumiferana × C. o. biennis (Table 1). Previous studies have found hybridization rates of 0.6 and 7% between C. fumiferana and various subspecies of C. occidentalis, based on proportions of hybrids detected using genomic data (Table 1; Brunet et al., 2013Brunet et al., , 2017Lumley & Sperling, 2011a). However, these studies did not sample at temporal or spatial scales sufficient to confidently assess hybridization across entire flight periods, as each was focused on resolving taxonomic relationships within Choristoneura, rather than the ecology and evolution of component species.
Selection against hybrid individuals exhibiting intermediate developmental rates could explain why few hybrid adults were detected in this study and why they are commonly collected near the peak flight of either parental taxon, rather than between flight period peaks. Hybrid larvae with intermediate developmental rates may miss their feeding "window of opportunity" (sensu Lawrence et al., 1997;Fuentealba et al., 2018). Past research has shown that larvae of conifer-feeding Choristoneura are in tight temporal synchrony with the annual bud flush of their host plant (Lawrence et al., 1997;Marshall & Roe, 2021;Nealis, 2012;Volney & Cerezke, 1992). If larvae develop too slowly, their likelihood of achieving adulthood drops precipitously (Lawrence et al., 1997 (Coates et al., 1994), putting intermediate larvae at a disadvantage at any one location. F 1 hybrids may also be infertile.
All hybrid individuals were assigned to the F 1 hybrid class with >99% certainty, with no F 2 or backcrossed offspring detected. This may be evidence that hybrid breakdown between these species is underway (Sperling, 1994). Hybridization of these Choristoneura species in laboratory settings will be required to resolve whether adapt to cooler environments, as this subspecies is typically found at higher elevations and latitudes than the other taxa (Dupuis et al., 2017;Shepherd et al., 1995).
Despite the reduced-representation methods applied here, wherein only a small fraction of individuals' genomes are sequenced, we still detected a few significant associations between genomic variation and environmental conditions within the spruce budworm species complex. The LFMM 2 candidate SNP in the C. fumiferana + C. occidentalis dataset was strongly associated with both mean temperature and precipitation of the warmest quarter of the year, during which Choristoneura diapause is initiated (Han & Bauce, 1998). Together, Bayescan and LFMM 2 results suggest that ecologically mediated selection and temporal isolation may be synergistically contributing to evolutionary divergences within the spruce budworm species complex and that either mechanism in isolation cannot account for reproductive isolation. We suggest that future work, using whole-genome sequencing techniques, should be completed to identify additional genomic regions that may be implicated in divergent selection among taxa and local adaptation to environmental conditions.

| Genomic architecture of Choristoneura
Our finding that 26 of 37 F ST outliers between C. fumiferana + C. occidentalis had linkage group matches on the Z sex chromosome (Table S8) suggests that adaptive divergence within Choristoneura is disproportionately sex-linked (Sperling, 1994). Sex-linked genes are often greater contributors to species-level divergence than autosomal genes (Baiz et al., 2020;Janz, 2003;Presgraves, 2018) and are known to evolve more quickly due to their hemizygous state in one sex ("faster X evolution", Bachtrog et al., 2009;Charlesworth et al., 1987). Further, several of the Z-linked candidate loci were associated with antifreeze genes that may contribute to overwintering success. Sex chromosomes of other Lepidoptera have been found to house many adaptive loci associated with diapause (Fu et al., 2015;Pruisscher et al., 2018Pruisscher et al., , 2020, demonstrating their importance for rapid evolution of the trait. Indeed, diapause frequency is the only consistent phenotypic difference that separates C. o. biennis from other conifer-feeding Choristoneura. However, it is also possible that higher rates of genetic drift associated with lower effective population sizes and the absence of recombination could account for high F ST values observed on the Z sex chromosome (Charlesworth et al., 1987;Ravinet et al., 2017). Further investigation of divergence of coding regions on the Z sex chromosome, such as those associated with diapause, are required to advance our understanding adaptive divergence in the C. fumiferana species complex.

| Applications to forest management
Our study effectively quantified temporal and spatial relationships between C. fumiferana and C. occidentalis across their zone of sympatry, demonstrating that phenological divergences are sufficiently large to contribute to reproductive isolation between the species.
However, some degree of phenological overlap was still observed between the species, providing opportunity for hybridization. Both C. fumiferana and C. occidentalis are destructive defoliators of conifers throughout North America, and several phenological models have been developed to detect and manage their outbreaks (e.g., Régnière, 1990;Régnière & Nealis, 2018;Régnière et al., 2012). Choristoneura fumiferana models have generally been parametrized using occurrence and phenological records from Ontario and Quebec, Canada, and extrapolated to the western edge of the species' range, but with reduced predictive power (Candau et al., 2019;Régnière et al., 2012).
Similarly, C. occidentalis models have been parameterized far from the zone of C. fumiferana-C. occidentalis sympatry (Nealis & Régnière, 2014;Régnière & Nealis, 2018). Outbreak models would benefit from the inclusion of occurrence and phenological data from the western and eastern edges of ranges of C. fumiferana and C. occidentalis, respectively, where the two species are sympatric. Our findings show that, due to overlap of flight periods within this zone of sympatry, species identity cannot be assigned based on collection time.
Individuals will need to be identified using genomic data to effectively parameterize outbreak models. Additionally, given that C. fumiferana and C. occidentalis are likely to increase their zone of sympatry in the northern boreal forest by mid-century (Régnière & Nealis, 2019;Régnière et al., 2012), understanding their potential for introgression and mechanisms of reproductive isolation gives insight into how they will interact at higher latitudes. Hybridization is likely to occur across an increasingly larger zone of sympatry as climatic conditions change, but we predict that these species will maintain their genomic integrity as they have in their present-day zone of sympatry.

ACK N OWLED G M ENTS
We greatly appreciate the loan of Unitraps and data loggers by Maya

CO N FLI C T O F I NTE R E S T
We declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences are available in fastq format in the National Center