Recurrent allopolyploidizations diversify ecophysiological traits in marsh orchids (Dactylorhiza majalis s.l.)

Molecular Ecology. 2023;32:4777–4790. | 4777wileyonlinelibrary.com/journal/mec Received: 3 January 2023 | Revised: 29 June 2023 | Accepted: 3 July 2023 DOI: 10.1111/mec.17070 ORIGINAL ARTICLE Recurrent allopolyploidizations diversify ecophysiological traits in marsh orchids (Dactylorhiza majalis s.l.) Thomas M. Wolfe1,2,3 | Francisco Balao1,4 | Emiliano Trucchi1,5 | Gert Bachmann6 | Wenjia Gu7 | Juliane Baar1 | Mikael Hedrén8 | Wolfram Weckwerth6,9 | Andrew R. Leitch7 | Ovidiu Paun1 1Department of Botany and Biodiversity Research, University of Vienna, Vienna, Austria 2Vienna Graduate School of Population Genetics, Vienna, Austria 3Department of Forest and Soil Sciences, University of Natural Resources and Life Sciences, Vienna, Austria 4Departamento de Biologia Vegetal y Ecologia, University of Seville, Sevilla, Spain 5Marche Polytechnic University, Ancona, Italy 6Department of Functional and Evolutionary Ecology, Molecular Systems Biology (MOSYS), University of Vienna, Vienna, Austria 7School of Biological and Chemical Sciences, Queen Mary University of London, London, UK 8Department of Biology, Lund University, Lund, Sweden 9Vienna Metabolomics Center (VIME), University of Vienna, Vienna, Austria Correspondence Thomas M. Wolfe and Ovidiu Paun, Department of Botany and Biodiversity Research, University of Vienna, Rennweg 14, A- 1030 Vienna, Austria. Email: thomas.m.[email protected] and [email protected] Funding information Austrian Science Fund, Grant/Award Number: W1225- B20 and Y661- B16 Handling Editor: Richard J. Abbott Abstract Wholegenome duplication has shaped the evolution of angiosperms and other organisms, and is important for many crops. Structural reorganization of chromosomes and repatterning of gene expression are frequently observed in allopolyploids, with physiological and ecological consequences. Recurrent origins from different parental populations are widespread among polyploids, resulting in an array of lineages that provide excellent models to uncover mechanisms of adaptation to divergent environments in early phases of polyploid evolution. We integrate here transcriptomic and ecophysiological comparative studies to show that sibling allopolyploid marsh orchid species (Dactylorhiza, Orchidaceae) occur in different habitats (low nutrient fens vs. meadows with mesic soils) and are characterized by a complex suite of intertwined, pronounced ecophysiological differences between them. We uncover distinct features in leaf elemental chemistry, lightharvesting, photoprotection, nutrient transport and stomata activity of the two sibling allopolyploids, which appear to match their specific ecologies, in particular soil chemistry differences at their native sites. We argue that the phenotypic divergence between the sibling allopolyploids has a clear genetic basis, generating ecological barriers that maintain distinct, independent lineages, despite pervasive interspecific gene flow. This suggests that recurrent origins of polyploids bring about a longterm potential to trigger and maintain functional and ecological diversity in marsh orchids and other groups. KEYWORDS allopolyploidy, Dactylorhiza, differential expression, ecological differentiation, photosynthesis, soil chemistry This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2023 The Authors. Molecular Ecology published by John Wiley & Sons Ltd. 4778 | WOLFE et al. 1 | INTRODUCTION Wholegenome duplication is a central force in evolution (e.g. Van de Peer et al., 2017), and ploidy increase has been estimated to be directly associated with about one in seven angiosperm speciation events (Wood et al., 2009). In animals, polyploidy is rarer, but multiple cases of currently polyploid insects, fishes, amphibians and reptiles are known (Otto & Whitton, 2000). Moreover, genomic evidence suggests all vertebrates, flowering plants and some fungi descended from polyploid ancestors (Otto, 2007; Van de Peer et al., 2017). Polyploidy merges multiple entire genomes in one nucleus, and if polyploids are derived from divergent species (as for allopolyploids), or even from distinct individuals, it can trigger a range of transgressive transcriptomic responses and transacting effects (e.g. Adams & Wendel, 2005; Burns et al., 2021). Polyploids are often reproductively isolated from diploid relatives, causing population densitydependent fitness costs known as minority cytotype exclusion (Husband, 2000; Levin, 2002). To establish and spread, polyploids are believed to benefit from certain reproductive characteristics, such as assortative mating and a breakdown of selfincompatibility (Fowler & Levin, 1984, 2016; Soltis et al., 2014). Polyploidy appears to have significant physiological and ecological consequences, which may lead to establishment in divergent environments that would otherwise be considered stressful for the diploid parents (Novikova et al., 2018; Van de Peer et al., 2017). Furthermore, palaeopolyploidy events tend to cluster around past periods of unstable conditions, including mass extinction events, which also suggests that stressful environments could systematically induce higher rates of polyploidy (e.g. Lohaus & Van de Peer, 2016). Novel ecological properties may arise in earlygeneration polyploids as a byproduct of large- and smallscale adjustments of genome organization and gene regulation that can occur with polyploidy, which with meiotic segregation can lead to transgressive characters, especially in allopolyploids (Paun et al., 2007). Several lines of evidence suggest that polyploids can show novel morphological and physiological characteristics compared with their lower ploidy relatives, with indirect consequences on complex ecological interactions, even between different trophic levels (Bomblies, 2020; Ramsey, 2011; Segraves, 2017). Information on ecophysiological consequences of polyploidy is currently limited (Scheriau et al., 2017; Segraves, 2017; te Beest et al., 2012). Genome multiplication is expected to change the scaling between nucleus and other cellular components, and positive relationships are generally observed between genome size and cell size, guard cell length and epidermal cell area, whereas a negative relationship with stomatal density is observed (Beaulieu et al., 2008; Doyle & Coate, 2019). A main ecophysiological implication of polyploidy is thought to be an increased nutrient requirement for nitrate and phosphate, necessary for DNA replication and chromatin packing in polyploids compared to diploids, especially in habitats with limited nutrient resources (Glennon et al., 2014; Guignard et al., 2016, 2017; Jeyasingh et al., 2015; Ramsey, 2011). Systemspecific physiological consequences of polyploidy identified through comparisons between polyploids and related diploids have also been shown, for example, increased salt tolerance (Chao et al., 2013), alterations in wateruse efficiency (Garbutt & Bazzaz, 1983) and in photosynthetic activity (e.g. Coate et al., 2012). Such outcomes may allow neopolyploids to range more broadly over the adaptive landscape and access adaptive peaks that may otherwise be unattainable through incremental changes at the diploid level (Wang et al., 2022). Established polyploids can exhibit larger distributions than their lower ploidy progenitors, especially across habitats that feature less interspecific competition, for example, previously glaciated regions (Novikova et al., 2018; Rice et al., 2019; Stebbins, 1971). Disentangling the ecophysiological implications of polyploidization in a natural context, relative to adaptation during later evolution, is important for understanding how polyploids spread across niches and their effects on biodiversity (Abbott et al., 2013; Jurado et al., 2019; Soltis et al., 2010). Recurrence in polyploidization is commonly observed during polyploid evolution (Soltis & Soltis, 1999). Examples include up to 20 independent origins of the allopolyploid Tragopogon miscellus, 13 origins for Draba norvegica and 46 for tetraploid Galax urceola (Segraves, 2017). Within the same population, recurrent formation of polyploids from different parental genotypes can greatly enhance genetic diversity in the nascent polyploid (Soltis et al., 2010; Soltis & Soltis, 1999). However, independent wholegenome duplications can happen in different ecological contexts too, potentially leading to local adaptation and the formation of polyploid species with different evolutionary trajectories (Jurado et al., 2019). Subsequent gene flow, independent assortment and recombination between polyploids can produce additional genetic variation, and it is unclear to what extent and in which conditions such sibling polyploids can maintain distinctiveness in the face of gene flow (Novikova et al., 2020; Paun et al., 2007; Soltis et al., 2010). In this study, we focus on the ecophysiological diversity released by recurrent allopolyploidization events in European marsh orchids (Dactylorhiza) in order to understand the context that allows occupation of diverse environments, potentially leading to the evolution and maintenance of distinct allopolyploid species. We investigate two sibling allopolyploids, Dactylorhiza majalis and D. traunsteineri that are estimated to have formed at different times in the recent half of the Quaternary (Brandrud et al., 2020; Hawranek, 2021). They arose from independent, but unidirectional allopolyploidization between D. fuchsii (as maternal genome) and D. incarnata (Figure 1a; Brandrud et al., 2020; Pillon et al., 2007). With DIY ABC (Cornuet et al., 2014) and genomewide RADseq data, the maximum age of the sibling allopolyploids has been estimated as c. 1730 generations for D. majalis and c. 920 generations for D. traunsteineri (Brandrud, 2019). However, based on the same data and coalescent inference, Hawranek (2021) estimated a maximum age of c. 104,000 generations for D. majalis and of c. 74,000 generations for D. traunsteineri. Apart from applying different demographic inference frameworks, the two estimates used different approaches to separate the homoeologs of the allopolyploids. Due to the considerable difference in the age estimates, further studies are required to refine and validate these estimates. Independent on their exact age, it is clear that D. majalis is older (Brandrud et al., 2020) and has survived in one or few neighbouring refugia in Central Europe during the last ice age 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 4779 WOLFE et al. (Balao et al., 2016). Dactylorhiza traunsteineri is in turn younger and more heterogeneous (Balao et al., 2016; Brandrud et al., 2020). It occupies almost exclusively previously glaciated ranges, suggesting an origin around the last glacial cycle. However, its regional mosaic of plastid DNA variants, some not found in the presentday representatives of its diploid parents, could support a prelast glacial origin also for D. traunsteineri (Nordström & Hedrén, 2008, 2009; Pillon et al., 2007). Since their origin, both allopolyploids have experienced a considerable increase in genome size, mainly driven by a putative tandem repeat (Eriksson et al., 2022). However, the repetitive landscapes of either allopolyploid do not generally show signals of a major genomic shock following allopolyploidization. Each of these allotetraploids occupies wide distributions across Europe. Dactylorhiza majalis ranges from the Pyrenees to southern Scandinavia, whereas its younger sibling D. traunsteineri has a more disjunct distribution around previously glaciated areas: in the British Isles, in northern Europe, and also in the Alps where it can grow in sympatry with D. majalis (Figure 1b). Both allopolyploids are pollinated through food deceit by naïve bumblebees, which can promote interspecific crosspollination (Hedrén & Nordström, 2018). Indeed, based on a wide range of genetic data, frequent and ongoing gene flow has been documented between the two allotetraploid species in the Alps, but with pervasive traces affecting also populations outside strictly sympatric localities (Balao et al., 2016; Brandrud et al., 2020; Hawranek, 2021; Hedrén et al., 2008). A welldocumented example of a strongly admixed population is the type locality of D. traunsteineri close to Kitzbuhel, Austria (i.e. locality ALP9 in Balao et al., 2016 and, respectively, at8 in Hawranek, 2021), whose individuals exhibit characteristic morphological traits, but generally contain a large portion of D. majalis alleles. The two sibling allopolyploid species have been suggested to have distinctive ecologies, even at localities where the two allopolyploids are close neighbours (Paun et al., 2011). However, no quantitative data are available on the type and magnitude of ecological differences. Dactylorhiza majalis is represented by robust, broadleaved plants with large inflorescences and typically grows in mesic to moist meadows, whereas D. traunsteineri is slender with narrow leaves and a fewflowered inflorescence, and it generally occurs in marshes and fens with a continuous influx of underground water; both species can benefit from mowing or grazing (Djordjević et al., 2016; Janečková et al., 2006; Sletvold et al., 2010; Wotavová et al., 2004). Dactylorhiza traunsteineri is commonly thought to be a poor competitor and tends to grow in sparse vegetation, whereas D. majalis generally occupies denser communities where there are more resources (Paun et al., 2010). As established, widespread species (Figure 1b), these sibling allopolyploids will have gone through the filtering effects of natural selection. Using these sibling marsh orchids, we present a multilevel, integrated investigation of physiology, ecology and gene expression to understand the phenotypic differentiation among them, and begin to disentangle the basis of their adaptation to divergent habitats. We further study differential gene expression in an outdoors common garden experiment, in order to understand the molecular context that allows the two allopolyploids to occupy distinct niches and to maintain phenotypic distinctiveness, despite pervasive gene flow (Balao et al., 2016; Hawranek, 2021). 2 | MATERIALS AND METHODS 2.1 | Soil and leaf elemental analyses We first quantified 18 distinct soil elemental components (N– NO3, N– NH4, Al, Ca, Cd, Cr, Cu, Fe, K, Mg, Mn, Mo, Na, Ni, P, Pb, S, Zn) and three soil features (pH, N%, C%) at 18 different localities across FIGURE 1 (a) The sibling allotetraploids Dactylorhiza majalis and D. traunsteineri have each originated from allopolyploidizations involving the diploid D. incarnata acting as the paternal parent and D. fuchsii as the maternal. The bubbles give the chromosome number and genome configuration of each species. Plant illustrations modified from Nelson (1976). (b) Curated D. majalis (red) and D. traunsteineri localities (blue) extracted from GBIF (accessed February 2018). (a) D. fuchsii D. incarnata FF 2x = 40 II 2x = 40 FFII 4x = 80 FFII 4x = 80 D. traunsteineriD. majalis (b) 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4780 | WOLFE et al. Europe (Table S1). The soil samples were collected in spring 2017 for nine localities of D. majalis (five from the Alps, two from the Pyrenees and two from Scandinavia), and nine localities of D. traunsteineri (three from the Alps, three from the Britain and three from Scandinavia). To obtain an accurate characterization, soil measurements were performed for each locality for two composite soil samples, each including three subsamples in approximately equal proportion. The soil was dried, and the elemental composition was measured using a Truspec CN Elemental analyser (LECO Corporation). Available nitrate and ammonium were extracted in 1 M KCl and measured on a plate reader photometer (Bio- Rad). In addition, a Direct Soil pH Meter HI99121 (Hanna Instruments) was used to measure soil pH at c. 7 cm depth in the close proximity of plants at 22 European localities in multiple years from 2010 to 2015 (Table S1). Significance of distribution differences for soil characteristics between the two allotetraploids was tested using permutation tests with the permTS function in the R package ‘perm’ (Fay & Shaw, 2010). To test if the soil elemental profiles specific for each allotetraploid translate to a differential chemistry in the plant, we further analysed leaf chemical content (i.e. C, N and P) of wild plants at 19 localities across Europe (Table S1). Leaf tissue was collected in silica gel at nine D. majalis localities (five from the Alps, two from the Pyrenees and two from Scandinavia) and at 10 D. traunsteineri populations (three from the Alps, five from Britain and two from Scandinavia). To measure N and C in the leaf tissue, 1– 2 mg ovendried ground sample was mixed with a urea solution and the mixtures were quantified on a mass spectrometer (Integra 2). One blank and two case- Rs as reference were also loaded for every 10 samples. To measure P, 0.8– 2 mg ovendried, ground tissue was mixed with 0.15 g potassium persulphate and 1 mL 1 N sulphuric acid and extracted at 121°C for 40 min. The P quantification was finally performed on a Skalar machine (Segment Flow Analyser, San++, Skalar Analytical B.V.). Significance was assessed using permutation tests as explained above. 2.2 | Macroenvironmental climatic niche characterization We explored if macroenvironmental differences exist between the two allotetraploids. Occurrence information for the two allotetraploid species was collected from Global Biodiversity Information Facility (GBIF.org., 2020; https://www.gbif.org): A total of 958 D. majalis and 631 D. traunsteineri localities were initially extracted using the package ‘rgbif’ (Chamberlain et al., 2017) in R software v.3.4.2. The localities were manually curated based on expert knowledge, duplicates were removed and locally dense sampling was reduced by thinning the records to one per 10 km2 grid cell size, retaining 298 localities for D. majalis and 393 localities for D. traunsteineri (Figure 1b). We further extracted data for 36 environmental variables from the CHELSA database (Karger et al., 2017) and the ENVIREM database (Title & Bemmels, 2018) for these localities at ~1 km2 resolution. We used the variance inflation factor (VIF; threshold of 10) to remove highly correlated environmental variables in the data set. This retained only 15 of the 36 environmental variables, summarizing isothermality (BIO3), temperature seasonality (BIO4), mean temperature of the wettest quarter (BIO8), mean temperature of the driest quarter (BIO9), precipitation seasonality (BIO15), precipitation of warmest quarter (BIO18), precipitation of coldest quarter (BIO19), index of the degree of water deficit below water need (aridityIndex), sum of mean monthly temperature for months with mean temperature greater than 5°C multiplied by number of days (growingDegDays5), minimum temperature of the warmest month (minTempWarmest), mean monthly potential evapotranspiration (PET) of the coldest quarter (PETColdestQuarter), mean monthly PET of the driest quarter (PETDriestQuarter), monthly variability in PET (PETseasonality), mean monthly PET of the wettest quarter (PETWettestQuarter) and topographic wetness index (topoWet). To compare the macroenvironmental niche of the sibling allotetraploids, we calculated the kernelsmoothed densities of all occurrence data along environmental axes from a principal component analysis (PCA- env; Broennimann et al., 2012). For background, we extracted environmental data from 3000 spatially random localities within a buffer of 100 km of land around the occurrence localities. Then, we performed a niche equivalence test (Broennimann et al., 2012) to compare the observed environmental niche overlap (Schoener's D- statistic) with a null distribution based on 1000 replicates and an environmental grid resolution of 500 × 500 pixels in the R package ecospat v. 3.2. (Di Cola et al., 2017). In the niche divergence test, we used Student's ttests to evaluate the differences in PCA scores between the two allotetraploids. 2.3 | Photosynthesis and gas exchange To complement the transcriptomic and ecological investigations, we phenotyped accessions of both allopolyploids in the wild for leaf photosynthetic parameters using a MultispeQ v.1.0 (Kuhlgert et al., 2016). For a total of 38 individuals of D. traunsteineri and 20 individuals of D. majalis, we quantified photosynthetic parameters (i.e. SPAD relative chlorophyll; LEF, linear electron flow; PAR, photosynthetically active radiation (light intensity); Phi2, quantum yield of Photosystem II; NPQt, nonphotochemical quenching; PhiNPQ, ratio of incoming light that goes towards nonphotochemical quenching), against momentary air temperature and humidity. We took two replicate measures on different leaves for each accession in the field at two localities in the Alps, where D. majalis and D. traunsteineri are within a few hundred metres of each other (i.e. ALP9, ALP13; Table S1) and one population where only D. traunsteineri grows (ALP8). Over a period of three sunny end- of- May days, we randomized as much as possible the effect of the time of the day by alternating measurements between species. Because we assumed that our data could be clustered or correlated due to our experimental setup, we fitted the photosynthesis parameters measured on the MultispeQ device with mixed linear models implemented in the lmer function of the lme4 R package (Bates et al., 2015). We first tested the normality of the data using 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 4781 WOLFE et al. the descdist function of the fitdistrplus R package (Delignette- Muller & Dutang, 2015), and then fitted a null model with each measured photosynthesis variable as outcome and individuals, time of the day and date of the measurements as random effect explanatory factors. This gave us a null distribution for each photosynthesis outcome against which we tested the full distribution with the species factor as an explanatory variable fixed factor. We tested for the significance of species effect on the fitted data using the likelihood ratio test implemented in the anova R function (Bolker et al., 2009). We considered each photosynthesis response as independent to each other, so we have one model per factor. This might not be ideal for modelling, since different photosynthesis outcome could presumably be correlated. We did not explicitly test for such interactions since we aimed to obtain results that could reflect all facets of photosynthesis and could provide the context for interpreting gene expression and other phenotypic differences between the sibling allopolyploids. The detailed results of the ANOVA test between the null hypothesis distribution as described above, and the full distribution including species and an explanatory factor, can be found in the Table S4. Finally, we measured leaf stomatal conductance for water and net CO2 exchange in two replicates of each species in a controlled lab environment using a portable gas exchange fluorescence system GFS- 3000 (Heinz Walz). During the experiment, the conditions were set to 25°C during day and 23°C in the night, 50% air humidity, 12 h photoperiod at 400 μmoles m−2 s−1 with a LED simulated sunlight spectrum set to 5800 K, watering the pots daily to field capacity. 2.4 | Reference genome assembly and annotation To enable further transcriptomic analyses, we assembled de novo a draft genome for a diploid individual of D. incarnata from Austria (48°11.6′N, 16°29.0′ E) using 177.6 Gb PacBio long reads (i.e. c. 50.7× coverage). We chose to assemble a representative of the paternal species to the allopolyploids despite its larger genome size compared to the maternal D. fuchsii (3.6 pg vs. 2.9 pg; Aagaard et al., 2005; Eriksson et al., 2022), because of its marked reduced heterozygosity levels compared to D. fuchsii (Balao et al., 2017; Brandrud et al., 2020). PacBio library preparation and sequencing of 18 SMRT cells on a Sequel I instrument was performed at the sequencing facility of the Vienna BioCenter Core Facilities (VBCF; https://www.vienn abioc enter.org/). PacBio data are available from GenBank BioProject PRJNA934399. The assembly was performed with Canu v.1.8 (Koren et al., 2017) with default settings, including the recommended correctedErrorRate 0.045. Based on the PacBio reads, the raw assembly was scaffolded with SSPACE- LongRead v.1– 1 (Boetzer & Pirovano, 2014) and polished with Arrow v.2.3.3 (available from https://github.com/Pacif icBio scien ces/Genom icCon sensus). The genome was finally curated to minimize inclusion of both allelic contigs as haplotigs using minimp2 (Li, 2018) and the Purge Haplotigs pipeline (Roach et al., 2018). The draft D. incarnata genome v.1.0 (GenBank BioProject PRJNA934399) contains 3413 scaffolds and has an N50 size of 7.4 Mb, recovering a total length of 3.27 Gb. This corresponds to 94.3% of the estimated genome size (1C = 3.6 pg, Aagaard et al., 2005; Eriksson et al., 2022). Before gene annotation analyses, we generated a custom repeat library for repeat masking using RepeatModeler v.1.0.11 (http:// www.repea tmask er.org/Repea tMode ler/). De novo TE annotation was performed with RepeatMasker v.4.0.7 (https://www.repea t m a s k e r . o r g / R e p e a t M a s k e r / ) and the search engine set to NCBI (−e ncbi). This process masked as repeat sequences 83.7% of the genome reference. The genome was structurally annotated ab initio using Augustus (Stanke et al., 2006) and GeneMark- ET (Lomsadze et al., 2014), as implemented in BRAKER1 v.2.1.0 (Hoff et al., 2016) with the options - - softmasking = 1 - - filterOutShort. Mapped mRNA- seq data of the same accession were used to improve de novo gene finding. An additional source of extrinsic evidence was singlecopy protein sequences predicted by BUSCO v.3 (Simão et al., 2015) when run on the hard masked genome in genome mode with option - - long. The annotation was further improved using MAKER- P v.2.31.10 (Campbell et al., 2014) supplying gene models identified using BRAKER1. A transcriptome was assembled using Trinity v.2.4.0 (Haas et al., 2013) based on mRNA- seq data, and this has been used in MAKER- P as expressed sequence tags (EST). The annotation used as additional evidence (declared as altest) the transcriptome of Orchis italica (De Paolo et al., 2014) and protein sequences of Phalaenopsis equestris (Cai et al., 2015). Structural annotations identified a total of 52,665 proteincoding gene models. They were functionally annotated using Blast2GO v.5.2.5 (Götz et al., 2008) and the Viridiplantae database, recovering full annotations for 31,394 genes. Finally, we evaluated the completeness of the genome assembly by searching our gene models against the BUSCO v.3 Liliopsida odb10 data set (Simão et al., 2015). A total of 89.2% of the set of singlecopy conserved BUSCO genes were found within our annotated genes, out of which 15.3% genes were recorded as duplicated. 2.5 | RNA- seq analyses To examine the functional differentiation between the two sibling allotetraploids, we further searched for gene expression differences in a common garden setup. We transplanted 11 representative adult D. traunsteineri plants from nine localities throughout the distribution range and eight adult D. majalis plants from eight localities (see Figure S1 and Table S1), and grew them in an outdoors common garden setting for 2 years in Vienna, Austria, with the aim of removing carryover environmental differences. Every year, these perennial orchids store nutrients into a new tuber to support growth in the following year. Leaf tissue was fixed in RNAlater (Sigma) for all accessions in the morning of 14 May 2014 as quickly as possible and in a similar developmental stage. The samples fixed in RNAlater have been stored at −80°C until processing. The wet lab procedure followed the details given in a complementary study focusing on the diploid parents (Balao et al., 2017). The RNA- seq libraries were sequenced as directional, 100 bp pairedend reads with Illumina HiSeq at the Vienna 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4782 | WOLFE et al. BioCenter Core Facilities (VBCF; https://www.vienn abioc enter. org/). The RNA- seq data (Wolfe et al., 2023) are available from GenBank BioProject PRJNA317244. Read trimming was performed with Trimmomatic v.0.33 with LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50. For differential expression analyses, the qualitytrimmed RNA- seq reads were mapped to the D. incarnata v.1.0 reference genome with STAR v.2.5.2a (Dobin et al., 2013) using the gene models annotated. We checked for transcript 5′– 3′ coverage biases known to affect RNA- seq experiments with RSeQC (Wang et al., 2012). No substantial transcript coverage bias in our data was found (Figure S4). The table of counts was obtained using the featureCounts function in Rsubread R package using the full transcripts as features for which to count the reads (Liao et al., 2013). The RUVs normalization approach implemented in the RUVseq v.1.8.0 R package (Risso et al., 2014) was applied to the table of counts for removal of technical variance. RUVs use technical replicates as controls to perform factor analyses on the count matrix and return a linear model that is used to correct for unwanted variance in differential expression analyses. Three (k = 3) factors of unwanted variance were removed as this yielded the best overlap between technical replicates overall; k = 1, 2 and 4 were also tested. We checked these normalizations with Relative Log Expression plots (plotRLE function) and principal component analyses plots (plotPCA function) using the EDAseq v.2.8.0 R package (Risso et al., 2011). We then fitted the set of read counts using the glmQLFit function of EdgeR v.3.24.3 (Robinson et al., 2010) with a nonintercepting model for species and normalizing factors obtained from RUVs as fixed effect variables (McCarthy et al., 2012). The glmQLFTest function was applied to the fitted model and an FDR < 0.05 threshold for differential expression for the a priori species contrast was set. The genome annotation was then imported into R with topGO v.2.26.0 (Alexa & Rahnenfuhrer, 2018) and enrichments between pairs of species were done using a twosided Fisher's exact test against the set of retained transcripts with an FDR smaller than 0.05. Visualization and summary of overrepresented GO terms were performed using the REVIGO web program and results were exported into R. We also used the GOplot v.1.0.2 R package (Supek et al., 2011; Walter et al., 2015). In order to have the directionality of the differential expression for the enriched GO terms towards one species or another, we calculated the zscore for each GO category as zs c o r e = (#genesup − #genesdown)/#- genestotal. The table of counts and scripts used for differential expression analyses is available at https://github.com/twolf e/dacty lorhiza. 3 | RESULTS 3.1 | The sibling allopolyploid marsh orchids are ecologically distinct Using soil chemistry analyses and permutation tests, we determined that relative to D. majalis, D. traunsteineri (Figure 2b) occurs in soils depleted of several macro- and micronutrients (Figure 2a and Table S2), in particular characterized by extremely low levels of available nitrate (NO3 −), but also phosphate (P) and potassium (K). In contrast, D. majalis (Figure 2c) occurs in richer soils, albeit its populations show in general a wide variance in soil elemental content (Figure 2a). No significant difference was recorded for available ammonium (NH4+) between the specific soils for the two Dactylorhiza allotetraploids (Table S2). Soil pH monitored across multiple years across the distribution range was also found to be significantly different, with D. traunsteineri tending to grow in slightly more alkaline soils than D. majalis (Figure 2a and Table S2). Leaf elemental analyses (Figure 2d) showed that relative to D. majalis, the narrower leaves of D. traunsteineri (Figures 1 and 2) had significantly lower levels of N and P per weight unit, but not significantly less carbon (C), resulting in higher C:N and C:P ratios. Although overlapping, the climatic niches for the two allotetraploids were not identical (Niche equivalency test, Schoener's D = 0.32, p < .001; Figure 3). In particular, temperature seasonality (BIO4), mean temperature of the driest quarter (BIO9), precipitation seasonality (BIO15) and potential evapotranspiration of the driest quarter (PETDriestQuarter) appeared to contribute specificity to the climatic niches of each allotetraploid (Niche divergence test PC1, Student's t = 16.444, df = 962.31, p < .001; Figure 3). These results are in agreement with D. traunsteineri generally growing in more northern areas (Figure 1b), but also in wetter habitats than its sibling, D. majalis. 3.2 | The sibling allopolyploid marsh orchids are physiologically distinct Compared to D. majalis, the narrower leaves of D. traunsteineri (Figures 1 and 2) showed significantly less relative chlorophyll content in native conditions (Figure 4 and Table S4). D. traunsteineri used 82.0 ± 38.7 μmol photons m − 2×s − 1 more incoming light (400– 700 nm) to drive photosynthesis (PAR), compared to D. majalis. However, the quantum yield of PSII (Phi2), that is, the ratio of excited electrons that go into the PSII, was significantly lower in D. traunsteineri than in D. majalis. Finally, the amount of excess energy (i.e. excited electrons) regulated away from photosynthetic processes in order to reduce damage to the plant (i.e. nonphotochemical quenching, phiNPQ) was significantly higher in D. traunsteineri (Figure 4 and Table S4). Despite the proximity of the two species at the sympatric localities and our randomized sampling, the MultispeQ measurements uncovered a difference in the ambient conditions around the two allopolyploids. In accordance with its wet and exposed habitats (i.e. low surrounding vegetation), D. traunsteineri experienced on average 3.5% ± 1.1 (std. error) more ambient humidity (p = .0045) and a 1.4 ± 0.4°C higher temperature (p = 2e- 10) than D. majalis. Regarding leaf stomatal conductance for water and net CO2 exchange, we found that both allotetraploids maintained their stomata open at night, and hence, gas exchange through stomata occurred throughout the night (Figure 5), likely using 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 4783 WOLFE et al. evapotranspiration to maintain xylem transport to accommodate for nutrient scarcity. However, during the night, D. majalis showed a negative net CO2 (i.e. a release of CO2) in contrast to D. traunsteineri. In addition, during daytime D. traunsteineri transpires c. 20% less water, and this higher waterusage efficiency complements a c. 50% higher CO2 net assimilation rate compared to D. majalis. What this means for the enzyme apparatus and the light to energy conversion of the species needs to be elucidated in future research. As the transcriptomic aspects and gas exchange results were obtained in a common garden, it appears likely that these physiological features are genetically encoded. 3.3 | Distinct ecophysiologies are associated with a complex suite of transcriptomic differences Principal component analysis (PCA, Figure S2) based on the levels of gene expression produced partly overlapping clusters of the FIGURE 2 (a) Divergent preference for selected soil characteristics between two sibling Dactylorhiza allopolyploids: D. majalis (red) and D. traunsteineri (blue). Soil elemental profiling for available nitrate (NO3 −), phosphate (P), potassium (K), aluminium (Al), chromium (Cr), nickel (Ni) and soil pH across multiple European populations. (b) Dactylorhiza traunsteineri growing alongside carnivorous sundews (Drosera spp.) at a locality on Gotland island, Sweden. (c) Dactylorhiza majalis, growing in a typical meadow at a site near Örup, Skåne, Sweden. (d) Leaf elemental profiling for wild individuals of D. majalis (red) and D. traunsteineri (blue) at multiple European localities. Data for nitrogen (N), carbon (C), phosphorus (P) and their pairwise ratios in leaf tissues. For (a, d), the data have been normalized to a 0– 1 range using feature scaling only for visualization purposes, where xnorm = (x − xmin)/(xmax − xmin). The tabular form summarizes the raw data: Δ, mean difference; M, average; SD, standard deviation. The pvalues are for 1000 permutation tests. NO3 -PKAl Cr Ni pH Normalized value Component D. majalis (M±SD) D. traunsteineri(M±SD) Δp-value 1.2e-06 2.0e-03 6.4e-04 2.7e-02 3.5e-02 2.7e-02 NO 3 - (µg/g) P (µg/g) K(mg/g) Al (mg/g) Cr (µg/g) Ni (µg/g) pH 33.5 ± 32.5 983.5 ± 288.7 1.5± 0.8 12.4 ± 10.6 16.9 ± 13.4 19.3 ± 13.1 6.9± 0.5 1.5 ± 3.4 586.9 ± 288.5 0.8± 0.3 4.9± 4.7 8.4± 4.3 10.4 ± 6.4 7.2± 0.6 32.0 396.6 0.7 7.5 8.6 9.0 0.31.5e-06 1 0 (a) (b) NCP N:P C:N C:P Componen tΔ p-value 8.6e-03 2.0e-01 1.3e-02 N (mMol/mg) C (µMol/mg) P (mMol/mg) N:P C:N C:P 2075 ± 369 37.2 ± 1.3 125.0 ± 43.6 17.8 ± 4.6 18.5 ± 3.5 334.3 ± 125.7 1611 ± 346 36.4 ± 1.5 77.7 ± 36.3 23.5 ± 7.6 23.8 ± 6.5 575.5 ± 285.7 464.0 0.8 47.3 5.6 5.3 241.3 Normalized value 1 0 D. majalis(M±SD) D. traunsteineri(M±SD) 5.4e-02 2.6e-02 3.4e-02 (d) (c) 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4784 | WOLFE et al. two allopolyploids, and revealed some geographic signal, in particular within D. traunsteineri, likely reflecting its disjunct distribution (Figure 1b). Differential gene expression tests revealed that 316 genes (2.5% of those retained after countper- millionbased filtering) were significantly higher expressed in D. majalis, and 343 (2.7% of retained genes) significantly higher expressed in D. traunsteineri (Figure 6a and Figure S3). Photosynthesis (GO:0015979) and related processes were by far the most affected by differential FIGURE 3 Principal component analysis of nocollinear macroenvironmental variables performed for GBIF localities and background. The distribution of the selected variables loading on the two main axes is given on the left, with inertia explained by colour of the arrows according to the legend. Bioclim codes are explained in the text. On the right, density contour plots encompass occurrence points of the D. majalis (red) and D. traunsteineri (blue) in the 2D environmental space. The occurrence density curves for each axis showed greater environmental divergence in PC1 (Δ = 1.60, p < .001) compared to PC2 (Δ = 0.61, p < .001). 0 1 scaled −6 −3 0 3 −2.5 0.02.5 5.0 PC1 (34.8%) PC2 (31.4%) 01 scaled BIO3 BIO4 BIO8 BIO9 BIO15 BIO18 BIO19 aridityIndex growingDeg Days5 minTemp Warmest PETColdest Quarter PETDriest Quarter PET seasonality PETWettest Quarter topoWet 0 1 FIGURE 4 Photosynthetic characteristics for the sibling allotetraploids D. majalis (red) and D. traunsteineri (blue) at native localities in the Alps. LEF, linear electron flow; PAR, photosynthetically active radiation (light intensity); Phi2, quantum yield of photosystem II; NPQt, nonphotochemical quenching; PhiNPQ, ratio of incoming light that goes towards nonphotochemical quenching. The boxplots show normalized values only for visualization purposes calculated as in Figure 2. The tabular form gives the details of the likelihood ratio between the null model with time, date and individual measurements as random variables, and the full model with species added to the model as a fixed variable. SD, standard deviation. LEF PAR NPQtPhi2 PhiNPQ Chlorophyll Normalized value 1 0 Component p-value LEF PAR Phi2 PhiNPQ Chlorophyll NPQt Intercept (±SD) 50.1 ± 8.1 244.7 ± 39.8 0.48 ± 0.02 0.21 ± 0.03 42.9 ± 1.6 0.73 ± 0.2 Species effect (±SD) 10.8 ± 5.7 82.0 ± 38.7 -0.06 ± 0.02 0.08 ± 0.02 -5.3 ± 2.1 0.34 ± 0.13 3.0 4.2 9.8 8.8 6.1 χ 2 5.9 8.0e-02 4.1e-02 1.8e-03 2.9e-03 1.3e-02 1.5e-02 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 4785 WOLFE et al. expression between D. majalis and D. traunsteineri (Table S3). In particular, relative to D. majalis, D. traunsteineri shows an increased expression of the Lhca1– Lhca4 and Lhcb1– Lhcb4 antenna proteins that are part of the lightharvesting complex (LHC), which captures and delivers excitation to photosystems (Figure 6b; Li et al., 2004). Relatively higher expressed in D. traunsteineri are also several components of the oxygenevolving complex of the Photosystem II (PSII; Figure 6b) that are responsible for catalysing the cleavage of water to oxygen, protons and electrons (Raymond & Blankenship, 2004). Dactylorhiza traunsteineri has also higher expression in the Photosystem I (PSI) reaction centres V (psaG), X (psaK), N (psaN) and O (psaO) (Figure 6b), which are mediating the primary function of the PSI (i.e. electron transfer from plastocyanin to ferredoxin; Jensen et al., 2004). Notably, 10.9% of the differentially expressed genes (i.e. mostly higher expressed in D. traunsteineri) are involved in oxidation– reduction, a process that also shows one of the strongest overrepresentations in the enrichment test (Figure 6a). Also related to photosynthesis, several genes linked to protein– chromophore linkage and in the response to light stimulus are higher expressed in D. traunsteineri compared to D. majalis (Figure 4a). Chlorophyll biosynthesis processes appeared to be higher expressed in D. traunsteineri, whereas chlorophyll catabolic processes higher expressed in D. majalis. Several processes related to cell wall organization are found to be in general higher expressed in D. traunsteineri, including some linked to pectin catabolic processes and xyloglucan metabolism. Finally, some processes related to response to abiotic stress have been identified to be differentially expressed between the two species, including six genes involved in the response to heat and four linked to the response to cold. 4 | DISCUSSION Although the sibling species D. majalis and D. traunsteineri arose from recent unidirectional allopolyploidization events (Figure 1a; Brandrud et al., 2020; Hawranek, 2021; Pillon et al., 2007) and share part of their distribution area in Central Europe (Figure 1b), our investigations uncovered major and intertwined ecological, physiological and transcriptomic differences between them, despite not estimating directly the fitness of either species in the alternative habitats. In particular, whereas D. majalis generally occurs in mesic meadows, D. traunsteineri is found in environments with extremely poor mineral nutrition and low competition (Figure 2a; Paun et al., 2010), for example, often occurring alongside carnivorous plants such as sundews (Drosera spp.; Figure 2b). Importantly, the niches of the sibling allopolyploids are characterized by a marked difference in available NO3 − content (Figure 2a), but not NH4 +. Ammonium can provide a source of nitrogen (N) for plants, particularly in waterlogged soils, but the metabolic processes necessary to assimilate this form of N are different from a nitratebased metabolism (Xu et al., 2012). In contrast to NO3 − which is processed in the leaves, NH4 + is metabolized directly in the roots, a process that requires sugars which have to be transported from leaves to roots. Ammonium uptake can affect plant foraging of other nutrients, in particular, K, Ca and Mg. The specific nutrientpoor soils preferred by D. traunsteineri (Figure 2a) and the nitrogen deficiency quantified in its leaves in natural conditions (Figure 2d) are likely to be responsible for its low relative chlorophyll content (Figure 4), as nitrogen is required for chlorophyll synthesis (e.g. Allison et al., 1997). We observed a relatively higher expression of genes involved in light harvesting (Figure 6) in D. traunsteineri compared to D. majalis in an outdoors common garden setting, FIGURE 5 Examples of leaf stomatal conductance for water vapour (pink line, mmol m−2 s−1 on the rightside Y- axes) and net CO2 assimilation rate (blue line, μmoles m−2 s−1 on the leftside Y- axes) for an individual of D. majalis (a) and one of D. traunsteineri (b). The levels of stomatal conductance and CO2 exchange fluctuate with a 12- h day/night cycle. Stomatal conductance is lower in D. traunsteineri which shows that its leaves do not open stomata as much as D. majalis, resulting in less water loss (higher water use efficiency) for nutrient acquisition via xylem. At night, the net CO2 exchange is slightly negative for D. majalis (red arrow) meaning that it releases more CO2 through respiration during nighttime, in contrast to D. traunsteineri (blue arrow), which also features 50% higher CO2 assimilation rate during daytime. Plant illustrations modified from Nelson (1976). 1365294x, 2023, 17, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/mec.17070 by Universidad De Sevilla, Wiley Online Library on [13/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License