Seven conifer species were sampled, namely Picea glauca, Picea mariana, Pinus strobus, Pinus banksiana, Abies balsamea, Larix laricina, and Thuja occidentalis. For each species, seeds from ten provenances were obtained from the National Tree Seed Center (Fredericton, New-Brunswick, Canada), paying special attention to avoid provenances located within sympatric or paratric zones in species known to spontaneously hybridize with related taxa ( Figure 1 ; Supplementary Methods S1 ). For each species, between 15 and 18 diploid embryos representing ten distinct provenances (two seeds per provenance on average) were extracted and flash frozen in separate tubes. All species considered, a total of 118 individuals from 67 populations were sequenced. In addition, four provenances were randomly selected per species and one seed per provenance was dissected to extract the haploid megagametophyte, which was flash frozen in liquid nitrogen in separate tubes. Megagametophytes were barcoded individually, while embryos were pooled at equimolar concentration and barcoded as a single library, prior to the sequencing step ( Supplementary Methods S1 ). All SNPs identified subsequently in the pool of embryos and in at least one megagametophyte were discarded so to filter out paralogous non-mendelian SNPs (see next section 'SNP calling').

Total RNA was extracted using the MasterPure™ Plant RNA Purification kit (Epicenter, Madison, WI, USA). RNAs were sequenced in paired-end mode (2×125 bp) with an Illumina HiSeq 2500 (Methods S2). Raw sequencing data (reads) were deposited in the public database ENA (European Nucleotide Archive, https://www.ebi.ac.uk/ena/browser/home, accessions ERS16017105-ERS16017139 and ERS16049778-ERS16049791) and vcf files containing variants identified in each species were deposited in DRYAD ('https://doi.org/10.5061/dryad.p8cz8w9w1). We assessed the good representativity of the analyzed transcriptomes based on sequence similarity searches (Methods S3).

After sequence quality controls and filtering ( Supplementary Methods S4 ), reads were aligned to the reference transcriptomes of each species previously published (Van Ghelder et al., 2019). SNPs were called using HaplotypeCaller version 3.4, from the GATK tool kit (McKenna et al., 2010; DePristo et al., 2011) and subsequently quality-filtered ( Supplementary Methods S5 ). Since the megagametophyte is a haploid tissue in all seven conifer species investigated herein, SNPs identified within megagametophyte libraries were likely indicative of variations between paralogous gene sequences also occurring in embryos. Thus, these SNPs were considered as false-positives and were subtracted from SNPs identified in pools of embryos before subsequent analyses were carried out. In total, above 1.4 million raw SNPs were identified ( Supplementary Methods S5 , Supplementary Table S5.1 ). After removal of paralogous SNPs identified in haploid megagametophytes, around 867K SNPs remained ( Supplementary Table S5.1 ). Among them, ~398K high-quality SNPs ( Supplementary Methods S5 ), were located in coding sequences representing almost ~97K Open Reading Frames in total (~16K transcripts per species on average) ( Supplementary Table S1 ). All species considered and on average, ~82% of the transcripts and ~70% of the ORFs carried SNPs ( Supplementary Table S1 ).

Transcripts that were included in the comparison of molecular genetic diversity across species and the analysis of the functional annotations of polymorphic gene sequences had to exhibit an average coverage of 10 reads or more, at least 300 nucleotides with a depth of 10 or more reads, and contain coding sequences.

The length and read depth of transcripts were significantly heterogenous across the seven conifer species investigated ( Supplementary Methods S6 ). It is essential to control for such effects before analyzing SNP abundance differences across species (Eo and DeWoody, 2012). In this purpose, we applied a regression model assuming that the number of SNPs within transcripts follows a negative binomial distribution (Eo and DeWoody, 2012).

The model corrected efficiently for variations among transcripts depth and length, thus enabling a more rigorous comparison of the SNP abundance across species ( Supplementary Methods S7 ). Then, SNP rate heterogeneity among species was tested using a Kruskall-Wallis test. To group species based on their level of total SNP diversity, Kolmogorov-Smirnov and Cramer-von Mises tests were conducted.

Based on the premise that nonsynonymous substitutions are predicted to contribute more to adaptive evolution than synonymous substitutions (Stern and Orgogozo, 2008), one way to study molecular convergence is to compare the ratio of substitution rates at nonsynonymous (Ka) versus synonymous (Ks) sites in orthologous protein-coding sequences betwen species. Similar inferences can be drawn within taxa from gene SNP A/S ratios, since both ratios have been shown to be strongly positively correlated (Liu et al., 2008).

The SNP A/S ratio was then calculated for each gene as the number of SNPs per nonsynonymous site (A) divided by the number of SNPs per synonymous site (S). An adjusted SNP A/S ratio was used to include genes with no synonymous SNPs following the empirical logit principle (Agresti, 2013):

The SNP A/S ratio was calculated over the longest open reading frame predicted for each transcript in each species. In order to identify putative positively selected genes (PSGs), we first retained those with A/S values exceeding 1, the threshold usually considered as evidence for positive selection (Kimura, 1983), and applied a subsequent filtering step to ensure that the probability that a transcript has an A/S ratios exceeding 1 by chance alone was lower than 5%, thus reducing much the number of false positives (see in Supplementary Methods S8 ). The resulting sets of PSGs were then used to analyze and compare the sequences and functional annotations of these deemed PSGs across the seven conifer species considered in the study.

Predicted protein sequences were clustered into orthogroups with OrthoFinder version 2.3.8 (Emms and Kelly, 2019) run with default settings. Functional annotations of ORFs were derived from sequence similarity searches conducted with blastp version 2.13.0 against Uniprot (E-value <e^-15) and PFAM (El-Gebali et al., 2019). Sequences were also assigned to Gene Ontology (GO) classes by using the mapping between the UniprotKB sequences and the GO terms. The heatmaps were generated using the pheatmap R package (Kolde, 2019).

Enrichment tests were conducted with the R (version 4.0.2) package topGO (version 2.42.0; Alexa et al., 2006; https://bioconductor.org/packages/release/bioc/html/topGO.html), in order to identify GO terms enriched among annotations of the genes with high SNP A/S values ( Supplementary Methods S9 ).

PSGs were identified in Brassica (Guo et al., 2017). Their Arabidopsis orthologs were retrieved (https://www.arabidopsis.org/) for a total of 621 sequences. PSGs were also identified in poplar (Lin et al., 2018). The Populus trichocarpa sequences were retrieved from PopGenIE.org. Sequences of the 2,047 conifer genes under positive selection were then compared at the protein level to poplar and Arabidopsis proteins. Overall, pairs of homologous sequences between these dicots and conifers were identified following a blastp search (E-value<1E-30). When one dicot gene sequence was found homologous to several conifer gene sequences, or when one conifer gene sequence was homologous to several dicot gene sequences, only the best match was selected.

This study enabled the identification of nearly 1.5 million of SNPs across the transcriptomes of seven conifers (Picea glauca and Picea mariana, Pinus strobus and Pinus banksiana, Larix laricina, Abies balsamea, Thuja occidentalis). However, we retained the ~867K SNPs with highest quality (see Methods S5) to conduct the subsequent analyses. Among them, ~398K high-quality SNPs were located in coding sequences representing almost ~97K Open Reading Frames (equivalent to ~16K transcripts per species, on average) ( Table 1 ). All species considered, ~82% of the transcripts and ~70% of the ORFs carried SNPs. In such transcriptome sequencing endeavor, it appeared important to minimize the effects caused by sequencing depth before undertaking any analysis of SNP diversity. We carefully adjusted the SNP diversity by both sequence length and sequencing depth before proceeding to data analysis and comparison across genes and across species (see Methods). After these adjustments, SNP diversity was estimated for each species and it was found significantly heterogeneous among the seven conifer taxa. A Kruskal-Wallis test revealed significant differences in rates of synonymous, nonsynonymous and total SNPs across the seven species ( Supplementary Table S2 ). Based on overall SNP diversity, three groups were delineated based on the results of Kolmogorov-Smirnov and Cramer-von Mises tests ( Figure 1 ; Table 1 ; Supplementary Methods S7 ). The group of species with the highest level of overall SNP diversity included the two Picea species, the group with the lowest diversity included Pinus strobus and Thuja occidentalis, while the three remaining species, Abies balsamea, Larix laricina, and Pinus banksiana, had intermediate overall SNP diversity ( Table 1 ; Supplementary Table S1 ).

Synonymous and nonsynonymous sites were identified and rates of synonymous (S) and nonsynonymous SNPs (A) were estimated to calculate the gene SNP A/S ratio (see Methods). A ratio above 1 is indicative of positive or balancing selection related to adaptive evolution (Kimura, 1983; Fay et al., 2001). Ratios above 1 were found in around 19% and 25% of genes depending on the species ( Supplementary Table S8.2 ). Moreover, an excess of nonsynonymous SNPs was significant in 2,047 genes (around 2% of genes within each species) ( Supplementary Table S8.2 ) that were considered as genes under putative positive selection (PSGs, Positively Selected Genes) for subsequent analyses. In PSGs, the mean A/S values were in the range of 2.70-2.89 ( Supplementary Table S8.2 ). The proportion of PSGs was highly correlated with the level of overall SNP diversity detected within each species (R² = 0.93; p-value < 0.01; Figure 2 ). It should however be noted that while high A/S ratios are indicative of positive selection, further investigation would be required to validate the role of positive selection.

Out of the 2,047 putative conifer PSGs, 73.5% (1,504 genes) had a significant match (blastp E-value <1E-15) with a SwissProt-Uniprot protein, a proportion consistent with other studies in conifers (Hart et al., 2020). Moreover, 932 genes (45.5%) had a match with 434 PFAM families (match E-value <1E-15). These annotated genes had a wide variety of functions ( Supplementary Figure S1 ). Among GO terms, 1,219 Biological Processes (BPs) were assigned to these 2,047 genes. The most represented processes were directly related to signal transduction, responses to biotic stresses, as well as various related processes ( Supplementary Figure S2 ). For instance, six terms describing plant responses to pathogens were found 576 times. Unsurprisingly, signal transduction, which is a common denominator of cell response to a stimulus, was the most predominant term (10.3% of the genes), along with defense response (7.7% of the genes) ( Supplementary Figures S1 - S2 ). The nicotinamide adenine dinucleotide (NAD) catabolic process was also highly represented (4.3% of the genes), which is consistent with the central role of NAD in plant defense responses (Pétriacq et al., 2013).

The complete dataset of gene sequences was successfully clustered into orthogroups, demonstrating both a high clustering capability at the intraspecific level (with few unassigned genes) and at the interspecific level (with few species-specific orthogroups) ( Supplementary Table S3 ). Altogether, 94.7% of the gene sequences were assigned to an orthogroup and the remainder were orphans, as previously observed in conifers (Gagalova et al., 2022). Out of the 16,982 orthogroups delimited in total, 8,647 contained gene sequences from the seven species, and 1,034 others contained gene sequences from all species except the more phylogenetically distant Cupresseae taxon Thuja occidentalis ( Figure 3A ). The number of species-specific genes in our dataset (hereafter referred to as ‘species-specific genes’ to simplify terminology) was low in Pinaceae taxa (between 3% and 9%) and higher (19%) inThuja occidentalis, which was expected given that this taxon belongs to the more divergent Cupressaceae family ( Figure 3A ). These species-specific gene sequences could either not be assigned to any orthogroup, or represented species-specific orthogroups in our dataset ( Supplementary Table S3 ). In contrast to the trend observed in the complete dataset, orthogroups derived from PSGs in each species showed a much lower overlap among species ( Figure 3B ). The vast majority of them were species-specific (total of 699; 68.7%), with only eight orthogroups shared across all seven species ( Figure 3B ; Table 2 ). No particular trend between gene convergence and phylogenetic relatedness among species was detected. For instance, shared PSGs were not more predominant within genera than among genera ( Supplementary Table S4 ).

Similar trends were observed at the gene family level. In the overall dataset, only 102 PFAM accessions (3.2%) were species-specific, indicating that protein families or domains were predominantly shared among the seven conifer species ( Figure 3C ). In contrast, in PSGs, a majority of PFAM accessions (407 accessions, 68.8%) were species-specific ( Figure 3D ). Nevertheless, the proportion of shared families among species increased as compared to that for orthogroups ( Table 2 ; Supplementary Figure S3 ).

Among GO terms, 4,504 BPs and 814 BPs were associated with the overall gene sequence dataset and 814 BPs for the dataset of PSGs ( Figure 3E ; Supplementary Figure S4 ). Species-specific BPs were few in the overall sequence dataset (total of 389, 8.6%) but relatively more abundant (333, 40.9%) in PSGs, although at a much lower rate than that observed for orthogroups or PFAM families. In spite of more convergence observed at the level of BPs, these results are indicative of a high level of functional diversity in PSGs ( Figure 3F ; Table 2 ).

The eight orthogroups of PSGs shared by all conifer species ( Figure 3B ) were homologous to sequences of known functions except one. They included three disease resistance genes, a member of pentatricopeptide repeat family, a C-type lectin receptor-like tyrosine-protein kinase, a cytochrome P450 and an ABC transporter family member.

The shared PFAM accessions included the NB-ARC domain (PF00931), the protein kinase domain (PF00069) and the TIR domain (PF01582), which are frequently found in combination in proteins involved in defense responses, as well as the cytochrome P450 family (PF00067), and UDP-glucoronosyl and UDP-glucosyl transferase family (PF00201) ( Figure 3D ; Supplementary Figure S3 ). The most abundantly represented and shared families by at least two species are detailed in Supplementary Figure S3 .

Regarding the gene ontologies of PSGs, 29 BPs (2.38% of the total number of BPs observed in PSGs) were shared by all conifer species ( Figure 3F ) and were categorized into eight groups ( Supplementary Figure S4 ). Eleven BPs were related to responses to biotic (like response to fungi, hypersensitive response) or abiotic stresses (like oxidative stress, cold) ( Supplementary Figure S4 ). Four terms were related to protein modification (like protein ubiquitination) and three others were related to growth and development (like pollen development). Main shared BPs also encompassed other mechanism important for plant life and survival such as seed development and germination or response to light ( Figure 4 and Supplementary Figure S5 ).

PSGs were enriched in several GO classes, including 16 BPs, 15 molecular functions and 5 cellular components ( Figure 5 ). Among them, five terms were enriched in the seven conifer species and are involved in defense against pathogens. Half of the enriched BPs were related to stress responses ( Figure 4 ). The enriched molecular functions were involved in several enzymatic activities or in binding ( Figure 5 ). Noticeably, terms involving nucleotidases, ADP or ATP binding were enriched in almost all seven species ( Figure 5 ).

Several gene families representing PSGs were involved in disease resistance. Among resistance genes, there were 72 genes homologous to run1 and 41 homologous to rpv1 conferring resistance to mildew ( Supplementary Table S5 ). Several gene families were involved in resistance against Pseudomonas syringae (rps2, rps5, rfl1, tao1, rpp3) or against viruses (eds1, tmvrn) ( Supplementary Table S5 , Figure 4 ). Among families involving PSGs shared by all species, the C-Lectin (CLEC) and the DRL28 protein are also involved in disease resistance. CLEC proteins have a diverse range of functions including cell-cell adhesion, immune response to pathogens and apoptosis.

Several genes belonged to the chitin pathway involved in fungal wall degradation ( Supplementary Table S6 ). They included several chitinases, an homologue of the LYK5 chitin receptor (Cao et al., 2014), and an homologue of CERK1 (Chitin Elicitor Receptor Kinase 1) required as a cell surface receptor for chitin elicitor signaling leading to innate immunity in response to biotic stresses (Hu et al., 2021).

Also, several PSGs were found among gene families implicated in the secondary metabolite pathways ( Figure 4 , Supplementary Table S7 ). They included proteins with roles in wood formation and defense response against insects such as the cytochrome P450 76T24 involved in the monoterpenoid synthetic pathway (Miettinen et al., 2014), the abietadienol/abietadienal oxidase (Ro et al., 2005), and the delta-selinene synthase (Steele et al., 1998).

Conifer PSGs had 54 homologs also positively selected in both poplar and Brassica. This set of 54 PSGs matched a total of 17 distinct Arabidopsis genes including five transcription factors, three glucosyltransferases, two peroxidases and a range of other gene families ( Table 3 ).

In addition, 384 conifer PSGs were homologous to PSGs in either poplar or Brassica. Their annotations were diverse ( Figure 6 ) and related to a variety of BPs ( Supplementary Table S8 ). The most represented BPs in this group were related to stress responses (164 genes) and metabolism (127 genes) ( Figure 6 ). Several genes were homologous to disease resistance genes including several members of the rps family and homologs to summ2 (SUPPRESSOR OF MKK1 MKK2 2) conferring resistance against Pseudomonas syringae, as well as lrks4 and lrks7 (lectin-domain containing receptor kinase involved in resistance response to the pathogenic oomycetes Phytophthora) ( Supplementary Table S5 ).

Several lines of evidence indicated that levels of total molecular genetic diversity differed substantially among the seven conifers ( Table 1 ; Supplementary Methods S7 ; Supplementary Table S2 ). Pinus strobus and Thuja occidentalis had the lowest overall SNP diversity in their transcriptome, the two spruces were the most diverse, and the three remaining species were intermediate ( Figure 1 ). Several factors related to the experimental design may result in an underestimation of the intraspecific diversity estimates reported herein. For instance, rare alleles may be missed out due to their low frequencies in sequencing pools. Moreover, diversity estimates rely only on gene expressed in the embryo, which nonetheless represents roughly 60% of the transcriptome (Rigault et al., 2011). However, our comparative analyses should not be biased by these limitations, as they would affect all species equally. Interpreting transcriptome-wide genetic diversity patterns is not straightforward because intraspecific variation results from the complex interplay between mutation rate, effective population size (long-term Ne, which depends itself on historical and demographic factors), and linked selection (the molecular genetic diversity-reducing effect of selective sweeps on neutral loci in linkage with loci under natural selection) (Ellegren and Galtier, 2016).

The seven conifer species analyzed here are not expected to have significantly different mutation rates, as they are all long-lived woody perennials (Petit and Hampe, 2006; Sung et al., 2012) and considering that our genetic diversity estimates derive from transcriptome-wide SNP data rather than a limited number of genes. However, in relation to the neutral theory of evolution which assumes that much of the standing genetic variation derives from neutral or nearly neutral mutations (e.g. Kimura, 1983; Ohta, 1992), part of the observed interspecific differences in molecular genetic diversity likely relates to historical effective population sizes (Bousquet et al., 1992). Indeed, the minimum historical population size (Ne) of Picea glauca and Picea mariana, the most diverse species group in this study, was estimated at ~100,000 or more individuals (Bouillé and Bousquet, 2005; Chen et al., 2010), while that of Pinus strobus, which belongs to the low diversity group, has been estimated to be an order of magnitude lower (5000 to 10,000 individuals; Zinck and Rajora, 2016). Likewise, Thuja occidentalis, the less genetically diverse species studied herein, also harbors a low N _e (Pandey and Rajora, 2012). There was also an apparent relation between the level of intraspecific SNP diversity and the geographical extent of the species range ( Figure 1 ), though the study of more species would be needed to confirm this trend. For instance, the two most diverse species have wide transcontinental distributions, while the two less diverse species are the most geographically restricted. As the number of glacial lineages of the Pleistocene era is usually positively related with range size in North American tree taxa (Jaramillo-Correa et al., 2009), widely-distributed species are likely to have retained larger historical population size and standing genetic variation than species with currently more restricted natural ranges.

Since our molecular genetic diversity estimates originate from transcriptomic data, it is very likely that selection also played a role in shaping molecular genetic diversity because the transcriptome mostly encodes functional information. Indeed, natural selection can constrain intraspecific molecular genetic diversity through selective sweeps (Smith and Haigh, 1974; Hermisson and Pennings, 2005). Although selective sweeps are presumably uncommon in plants (Wright and Gaut, 2005; Grossman et al., 2010) including in conifers (Palme et al., 2008; Pavy et al., 2012b; Eckert et al., 2013), evidence for selective sweeps has been reported previously in some conifers (Eckert et al., 2009; Namroud et al., 2010; Wang et al., 2020; De La Torre et al., 2021; Gagalova et al., 2022). The pool sequencing data obtained for each species in our study does not allow to test directly for the existence of such sweeps. However, it is unlikely that selective sweeps are a major determinant of transcriptome-wide genetic diversity in conifers because, on average, linkage disequilibrium decays rapidly within gene limits in conifer natural populations (e.g. Pavy et al., 2012b; De La Torre et al., 2017), therefore restricting the possible loss of neutral diversity surrounding selected loci. Hence, although both N _e and linked selection may have contributed to shaping molecular genetic diversity at the intraspecific level, the former is more likely to have been the main driver of differences in overall SNP diversity observed herein among the transcriptomes of the studied species.

The positive relationship observed between the overall SNP diversity of the transcriptome of each species and their proportions of PSGs ( Figure 2 ) suggests that standing genetic variation can constrain variation of more adaptive nature. It is well established that the fixation probability of favorable alleles in a finite population increases along with the effective population size by reducing the strength of genetic drift and therefore limiting the loss of beneficial alleles as well as the fixation of deleterious ones (Charlesworth, 2009). Consistently, our results show that species of presumably larger historical population sizes carry the most adaptive molecular variation. This is also in agreement with more common local adaptation in large plant populations than in small ones (Leimu and Fischer, 2008), which has been reported for several species investigated herein, namely white spruce (Namroud et al., 2008; Hornoy et al., 2015; Depardieu et al., 2021), black spruce (Prunier et al., 2011; 2012), jack pine (Cullingham et al., 2014), and eastern white pine (Nadeau et al., 2016). Because environmental adaptation is highly polygenic in conifers and involves heterogeneous gene responses (e.g. Hornoy et al., 2015; Yeaman et al., 2016; Depardieu et al., 2021), high standing genetic variation associated with large historical population size likely improves species adaptative potential by increasing the number of possible genetic trajectories to achieve adaptation. Hence, such flexibility may allow species to cope with a wider range of environmental conditions (i.e. gain the ability to colonize larger natural range and/or increase their ecological amplitude) and selective pressures (i.e. biotic and abiotic pressures encountered across their range).

Despite the high overlap among gene sequences of the seven conifer species ( Figure 3A ), we only found limited molecular convergence among their genes under positive selection ( Figure 3B ). Convergence appeared also limited at the protein family level ( Figure 3D ), a result consistent with the pattern of species-specific expansion of large paralogous gene families (reviewed by De La Torre et al., 2020) and high functional redundancy in conifers (Guillet-Claude et al., 2004; Bedon et al., 2010; Stival Sena et al., 2018; Van Ghelder et al., 2019). The extent of molecular genetic convergence among taxa would also be expected to increase with their phylogenetic proximity, as a result of shared ancestry (Losos, 2011; Storz, 2016). However, no pattern related to phylogenetic relatedness among taxa was evident, with no sign of increased convergence among the two pairs of congeneric taxa that would have diverged the most recently (divergence between Picea glauca and Picea mariana ~10 Mya (Bouillé and Bousquet, 2005); divergence between Pinus strobus and Pinus banksiana ~85 Mya (Leslie et al., 2018) ( Supplementary Table S4 ). Similarly, convergence was only marginally higher among the Pinaceae taxa that between the Pinaceae taxa and the Cupressaceae taxon despite their more recent divergence (Leslie et al., 2018) ( Supplementary Table S4 ). These observations suggest that each species followed a largely distinct adaptive path, and that adaptive convergence at the molecular genetic level appears to be limited in such reproductively isolated and phylogenetically distant conifers.

There is also strong evidence that the marked difference in sets of PSGs among species would be primarily driven by natural selection, rather than by stochastic processes such as mutation or genetic drift (e.g. Losos, 2011; Mosca et al., 2012; Storz, 2016). The low levels of convergence observed among the species gene sets under positive selection either indicate that gene functional redundancy would allow species to cope with similar selective pressures using alternative genes, and/or that species experienced heterogeneous selective pressures throughout their historical and extent natural ranges. Our data support the first hypothesis, as functional convergence among biological processes associated with genes under positive selection was quite higher that molecular convergence among these genes ( Figure 3 ). This partly decoupled pattern indicates that the studied species would have had sufficient metabolic and gene network flexibility to evolve alternative responses to the various selective pressures they faced under temperate and boreal climate regimes. It is also consistent with gene family expansions in conifers, implying some redundancy in gene functions (Guillet-Claude et al., 2004; Bedon et al., 2010; Pavy et al., 2012a; Stival Sena et al., 2018; Van Ghelder et al., 2019). This redundancy at the functional level may have had significant evolutionary implications for the persistence of these northern conifer species during millions of years, in the face of geological climate instability and in spite of demographic fluctuations. For instance, with the multiple glaciation cycles of the Pleistocene era in northeastern North America, signatures of demographic fluctuations such as bottleneck effects or founder events have been detected in various conifer species from this region, which would have implied more or less important losses of overall genetic diversity (Perron et al., 2000; Gamache et al., 2003; Jaramillo-Correa and Bousquet, 2003; Godbout et al., 2010; Namroud et al., 2010).

In addition to the gene redundancy hypothesis, it is possible that the low convergence observed among conifer PSG sets also reflects that the species have coped with specific long-term selective pressures throughout their large natural ranges. Some species such as the two Picea spp. have large ecological amplitude and a transcontinental range across which they may encounter a variety of selective pressures related to biotic and abiotic stresses (Nienstaedt and Zasada, 1990; Viereck and Johnston, 1990). Also, Picea glauca would prefer mesic sites (Hornoy et al., 2015) while Picea mariana could adapt to a larger variety of site conditions (Lo et al., 2024). Others have more distinct preferred habitats, such as wetter sites for Larix laricina (Cheliak et al., 1988), dryer sites for Pinus banksiana (Rudolph and Laidly, 1990), or both for the more extremophile Thuja occidentalis (Matthes-Sears and Larson, 1991), which can trigger specific adaptive responses and explain the finite extent of adaptive convergence in gene sets under positive selection among the conifer species studied.

Hence, taken together, our results suggest that the adaptive trajectories of these conifer species were likely shaped by the interplay of gene redundancy and heterogeneous selection landscape, and that these two drivers likely contributed to the low convergence observed in terms of gene sets accumulating nonsynonymous SNPs, but higher convergence at the functional level.

This pattern aligns well with other reports of low estimates of molecular genetic convergence among conifers. A study of adaptive traits in four alpine conifers from the Pinus, Abies and Larix genera identified only seven climate-associated genes shared by two or more species out of several hundreds of sequences analyzed (Mosca et al., 2012). Another study uncovered only 47 convergent genes (representing between 10% and 18% of all genes putatively under selection) involved in local adaptation in a spruce and a pine taxon that diverged ~150 Mya (Yeaman et al., 2016). Also, a transcriptome-wide survey of genetic variation in the two quite closely-related but ecologically contrasted Picea glauca and the coastal Picea sitchensis from the Pacific Northwest, revealed only 15 shared genes out of hundreds of genes showing selection footprints (Gagalova et al., 2022). Similar modest molecular genetic convergence was also reported in Angiosperms such as in Arabidopsis (Guggisberg et al., 2018; Preite et al., 2019), or among taxa from the Brassicaceae family (Rellstab et al., 2020).

Functional annotations of genes under positive selection were associated with a great variety of molecular functions and biological processes ( Supplementary Figure S1 ), in agreement with the polygenic nature of adaptive traits and the many empirical studies reporting a wide range of genes and functions underlying them in conifers (e.g. Prunier et al., 2011; Mosca et al., 2012; Hornoy et al., 2015; Depardieu et al., 2021).

Despite the large functional diversity observed, the core set of functions and processes shared by all conifer species analyzed revealed a clear pattern of shared adaptive evolution at the functional level. Indeed, half of the 29 shared BPs by all conifer taxa examined herein were related to environmental stress responses, and mechanisms related to defense against pathogens (responses to biotic stress, programmed cell death) were widely represented and enriched in genes under positive selection ( Figure 4 ; Figure 5 ). We found many homologs of disease resistance genes with the NB-ARC domain ( Supplementary Figure S4 ), several gene families involved in resistance to Xanthomonas, Pseudomonas or rusts ( Figure 4 ; Supplementary Table S5 ), and several genes encoding enzymes from the secondary metabolite pathways that are known to be involved in defense mechanisms as well as wood formation ( Supplementary Table S7 ). This indicates that selection pressures exerted by pathogens are likely ubiquitous in conifers and play a prominent role in their adaptation to changing environments. Similarly, several shared BPs were linked to abiotic stress response, and more specifically to water and oxygen stimuli (i.e. cellular response to hypoxia, response to oxidative stress, or response to water deprivation), suggesting that drought and flooding could also be drivers of adaptive evolution at the functional level in conifers.

Mechanisms with crucial roles in stress response such as RNA modification and regulatory mechanisms were also found in conifer genes under selection ( Figure 4 ). Regarding RNA modification, we identified 59 sequences encoding pentatricopeptide repeat (PPR) proteins, among which two were also reported as convergent adaptive genes in pine and spruce taxa ( Figure 4 , Supplementary Figure S4 ; Yeaman et al., 2016). The PPR proteins have fundamental roles in organelle biogenesis and function, being involved in photosynthesis, respiration, development and environmental responses (Barkan and Small, 2014). Positive regulation of programmed cell death, a process known to be involved in response to biotic and abiotic stresses in plants (Daneva et al., 2016), was also enriched in all Pinaceae species ( Figure 5 ). Transcriptional regulators included several MYB and WRKY gene sequences. Homologs of genes involved in seasonal transitions in Arabidopsis are likely prime targets of natural selection given that they contribute to the adaptation of plants to their environment, assuming that their functions are conserved across seed plants. Indeed, a FRIGIDA-like protein (accession FRL3_ARATH) was found among conifer genes under positive selection ( Figure 4 ). In Arabidopsis, it regulates phase transition during shoot, flower and seed development such that FRIGIDA gene sequences are required for the winter-annual habit (Michaels et al., 2004). The identification of several homologous regulators involved in survival in Angiosperms also suggests a possible key role in conifers and makes them a relevant class of genes to target in future molecular functional studies.

We uncovered 17 PSGs shared by conifers and two Angiosperm taxa (Brassicaceae and poplar) ( Table 3 ). This level of molecular convergence was higher than expected, given that Angiosperms and Gymnosperms (to which belong conifers) diverged ~350 Mya (Li et al., 2019). Disease resistance genes against pathogens are known to evolve rapidly in flowering plants (Meyers et al., 2005). Although well represented among conifer PSGs, they were not predominant among the core set of convergent genes between conifers and Angiosperms. In contrast, several genes encoding enzymes (transferases, peptidase, peroxidases), one PPR gene involved in RNA editing, and several transcription factors such as MYBs showed up in this set of genes ( Table 3 ). To our knowledge, such widely spread molecular signatures of adaptation have not been reported to date and may be interpreted as a sign that adaptive convergence at the molecular level, though limited, can take place at a very broad taxonomic level. Consequently, these genes represent valuable candidates for future evolutionary studies aiming to characterize molecular and functional convergence among seed plants.