Hybridization and sex determination

Salix is of interest for the study of sex-determination systems (SDSs), and these may relate to the phylogenetic relationships in the following way, making it important to clarify relationships among willows. Populus and Salix are sister genera in Salicaceae, and both are dioecious. It has therefore been suggested that this state was present in the common ancestor of these genera (Dai et al., 2014) before they diverged from each other about 40–45 Mya (Boucher et al., 2003; Wu et al., 2015). Despite the long period during which dioecy could have existed, the chromosomes carrying the sex-determination region (SDR) are homomorphic (van Buijtenen and Einspahr, 1959). In both genera the SDR systems are based on a single factor located on either the female (W) chromosome or the male (Y) chromosome (Renner and Müller, 2021). This factor is either female-specific expressed on the W chromosome or dominantly repressed by the male Y chromosome. Changes in the position of the SDR can occur relatively easily by translocation or new mutations (Renner and Müller, 2021). It is therefore possible that the SDR has evolved to some extent independently in these two genera, and that turnover events could have occurred, in which an established sex-determining system was replaced by a different one. Indeed, both male and female heterogamety are now known within the genus Salix. Studies of SDSs have mainly been based on members of subg. Vetrix, including Salix polyclona, S. suchowensis, S. viminalis, S. purpurea and S. triandra, which have female heterogamety (ZW) with physically extensive sex-linked regions (SLRs) on chromosome 15 in all species (reviewed in He et al., 2021b; Table 1). However, recent studies have demonstrated male heterogamety (XY) and SLRs on chromosome 7 in two species of subg. Salix, namely S. dunnii (He et al., 2021b) and S. nigra (Sanderson et al., 2021). Thus, the different SDSs may support a biologically important subdivision of Salix that could be used as a character in subgeneric classification (although sex determination has so far been studied in only a few species of Salix).

Table 1.

Summary of current information about sex determination systems in Salix, adapted from He et al. (2021b)

Species	Male or female heterogamety	Chromosome carrying the sex-determining locus	References
Salix clade
S. dunnii	Male (XX/XY)	7	He et al. (2021b)
S. nigra	Male (XX/XY)	7	Sanderson et al. (2021)
Vetrix clade
S. triandra	Female (ZW/ZZ)	15	Li et al. (2020)
S. purpurea	Female (ZW/ZZ)	15	Zhou et al. (2020)
S. suchowensis	Female (ZW/ZZ)	15	Hou et al. (2015)
S. viminalis	Female (ZW/ZZ)	15	Almeida et al. (2020)
S. polyclona	Female (ZW/ZZ)	15	L. He et al. (in prep.)

Open in a separate window

Recent studies have tended to decrease the number of subgenera, predominantly agreeing on the recognition of only two major clades, the Salix clade consisting of subg. Salix, Longifoliae and Protitea, and the Chamaetia/Vetrix clade which includes subg. Chamaetia and Vetrix along with Salix arbutifolia (Chosenia), S. cardiophylla (Toisusu) and sect. Amygdalinae. Furthermore, interspecific hybridization occurs mostly within the Vetrix or Salix clades (Wagner et al., 2021a), raising the question whether the two main clades with different SDSs have reproductive barriers.

Thus, in the present study we aimed: (1) to use WGR data to reconstruct the phylogeny of sampled specimens (this approach has never been applied to the reconstruction of willow phylogeny across the whole genus); (2) to test whether phylogenies based on autosomes, SLRs and chloroplast sequences are consistent with one another; and (3) to determine whether differences in SDSs between species of the two clades offer a reliable character for their subgeneric division, or not.

Phylogenetic analysis

Since different SDSs were identified in the two main diploid clades, we extracted the SNPs in autosomal regions from both nuclear datasets (excluding the entire chromosomes 7 and 15). A python script (https://github.com/zhangrengang/degeneracy) was used to generate all four-fold degenerate (FF-D) sites in the genomes of Salix dunnii and S. purpurea, based on their gene annotation (GFF3 files, Zhou et al., 2020; He et al., 2021b). We extracted the SNPs at these FF-D sites from the two nuclear autosomal datasets. Polyploidy can confuse phylogenetic analysis, especially if interspecific hybridization was involved in duplication of the genome, which can lead to incongruent reconstructions as well as reticulate evolutionary patterns (Alix et al., 2017). We therefore obtained 233 684 FF-D SNPs for all individuals, and 207 155 FF-D SNPs for diploids, from the autosome regions AR-dun (SNPs in autosome regions, using Salix dunnii as reference genome); similarly, we obtained 258 908 FF-D SNPs for all individuals and 230 084 for diploids from the AR-pur (SNPs in autosome regions, using S. pupurea as reference genome).

Phylogenetic relationships were inferred by a maximum-likelihood approach using RAXML v0.8.2.4 (Stamatakis, 2014) based on the concatenated sequences of six datasets: AR-dun (all individuals), AR-Di-dun (diploids only), CP-dun, AR-pur (all individuals), AR-Di-pur (diploids only) and CP-pur. Support values were calculated using 100 rapid bootstrap replicates (-f a option) based on the GTR + GAMMA nucleotide substitution model as used by Wagner et al. (2020) and He et al. (2021a) for the phylogenomics of willows.

The SDSs are not known for any polyploid willows. We therefore extracted SNPs from the SLRs only for the diploids. For the dataset using the S. dunnii reference, this region is X-LR, chr07:5675000–8880000 of chromosome 7, and for that using the S. purpurea reference it is the Z-linked region of chromosome 15 (Z-LR, Chr15Z: 2341099–6715814) (Zhou et al., 2020; He et al., 2021b). Because recombination suppression of SLRs may lead to gene duplications, we used only single-copy genes (SCGs) in the chromosome 7X and 15Z regions, identified using OrthoFinder (Emms and Kelly, 2019). Because of assembly quality and availability of phased SLRs, we used only the genomes of S. brachista, S. dunnii and S. purpurea to identify the chromosome 7X SCGs, and only S. dunnii and S. purpurea (excluding Chr15W) to identify the chromosome 15Z SCGs. Then we used the two SCG datasets to extract the SNPs in the X-LR and Z-LR. This yielded 1165 SNPs in 32 SCGs in X-LR, and 663 SNPs in 15 of 74 SCGs in Z-LR. We developed a custom script to obtain the sequences of each SCG in these datasets, and used ModelFinder (Kalyaanamoorthy et al., 2017) to select the best model under the Bayesian information criterion, and IQ-tree based on the best model (Minh et al., 2020) to construct individual gene trees. Species trees of the two datasets were estimated using ASTRAL (Zhang et al., 2018a).

Phylogenetic trees based on nuclear and chloroplast data

As we hypothesized that the species of the two main clades of Salix have different SDSs (see Introduction), we chose Salix dunnii and S. purpurea from the two clades as reference genomes in order to analyse SLRs separately from autosomal ones, and used a total of eight datasets to reconstruct phylogenies and species trees of Salix species (see Methods; Table 2, Figs 1–4; and Supplementary Data Figs S2–S5).

Table 2.

Details of the eight datasets used to conduct phylogenetic analyses

Dataset
Nuclear genome
Four-fold degenerate SNPs of autosomal region
AR-dun	233 684
AR-pur	258 908
Four-fold degenerate SNPs of autosomal region (diploids only)
AR-Di-dun	207 155
AR-Di-pur	230 084
SNPs of single-copy genes in sex-linked region (diploid only)
SCG-X-LR	1165
SCG-Z-LR	663
Chloroplast genome
CP-dun	416
CP-pur	988

Open in a separate window

Dun and X-LR used the Salix dunnii genome as reference; X-LR represents the X-linked region on chromosome 7 of Salix dunnii. Pur and Z-LR used the Salix purpurea genome as reference; Z-LR represents the Z-linked region on chromosome 15 of Salix purpurea. AR, CP and SCG represent autosomal region, chloroplast genome and single copy gene, respectively.

Open in a separate window

Fig. 1.

Phylogeny inferred for 70 taxa and 91 samples of the genus Salix and the outgroup Populus euphratica based on maximum-likelihood analyses of the AR-dun (using the genome of Salix dunnii as reference, excluding chromosomes 7 and 15) dataset (233 684 four-fold degenerate SNPs) using RAxML. The subgenera are according to Skvortsov (1999), Ohashi (2006) and Argus (2010). Fang et al. (1999) divided Salix species only at the sectional level, but not all of those sections are monophyletic groups (He et al., 2021a). Thus, only species that were confirmed to belong to particular subgenera, based on Skvortsov (1999), Ohashi (2006) and Argus (2010), were coloured, accordingly. The scale bar refers the mean number of nucleotide substitution per site. Numbers above branches refer to bootstrap support values.

Open in a separate window

Fig. 4.

Phylogeny inferred for 70 taxa and 91 samples of the genus Salix and the outgroup Populus euphratica based on maximum-likelihood analyses of the CP-dun (using the chloroplast genome of Salix dunnii as reference) dataset (416 SNPs) using RAxML. The scale bar refers the mean number of nucleotide substitution per site. Numbers above branches refer to bootstrap support values.

Five major clades and one paraphyletic group were found by analyses of putatively autosomal sequences, using both the AR-dun and AR-pur datasets, which include both diploid and polyploid species (Fig. 1; Supplementary Data Fig S2). Clade Vetrix includes four subclades: C1, C2, C3 and C4. Salix includes one subclade (C6) and a small polyploid group, C5, consisting of Salix shihtsuanensis, S. chienii and S. matsudana, plus three branches with S. cf. fragilis, S. pentandra, S. paraplesia and S. lucida, paraphyletic to the Vetrix clade. Individuals of the same species are well supported as monophyletic in the two trees, except a few species or complex with multiple ploidy levels (S. polyclona, S. shihtsuanensis), allopolyploid origin (S. opsimantha) or unclear morphological boundaries (S. luctuosa and S. hypoleuca; S. fargesii and S. moupinensis).

Compared to the complete autosome datasets, the data from just the diploid species [AR-Di-dun, AR-Di-pur, and SLR datasets (SCG-X-LR and SCG-Z-LR)] are expected to be less affected by the problems outlined above for polyploid species, such as paraphyly and reticulate evolution. All analyses of the diploid data support separate monophyletic groups of Vetrix and Salix species (Figs 2 and ​and3;3; Supplementary Data Figs S3 and S4). Due to lower numbers of SNPs, the SCG-X-LR and SCG-Z-LR species trees did not distinguish C1 and C2, but revealed C3 as sister to the C1 and C2 groups. The AR-Di-dun and AR-Di-pur analyses revealed a similar topology of C1, C2, C3, and C4 of clade Vetrix as with the whole autosome datasets.

Open in a separate window

Fig. 2.

Phylogeny inferred for 38 diploid taxa and 49 samples of the genus Salix and the outgroup Populus euphratica based on maximum-likelihood analyses of the AR-Di-dun (using the genome of Salix dunnii as reference, excluding chromosomes 7 and 15) dataset (207 155 four-fold degenerate SNPs) using RAxML. The scale bar refers the mean number of nucleotide substitution per site. Numbers above branches refer to bootstrap support values.

Open in a separate window

Fig. 3.

Phylogeny inferred for 38 taxa and 49 diploid samples of the genus Salix and the outgroup Populus euphratica based on the SCG-X-LR dataset (using the genome of Salix dunnii as reference, 32 single-copy genes, 1165 SNPs) analysed with ASTRAL species tree methods. Values above branches refer to local posterior probabilities.

Both trees based on chloroplast genomes, CP-dun and CP-pur, again included two major well-supported clades, Salix and Vetrix (Fig. 4; Supplementary Data Fig. S5). Bootstrap support values for most of the internal nodes within the major clades were low, and no subclades were clearly defined within Vetrix, with the notable exception of subclade P1, Salix triandra. Representatives of subg. Chamaetia and Vetrix were intermixed in one highly supported clade. In contrast, the Salix clade was fully supported as distinct; all species previously recognized as members of subg. Salix (including S. exigua and S. interior) fell within it, in two well-supported subclades.

Reproductive barriers via SDSs in willows

In recent years, both male and female heterogamety have been found within Salix (see Table 1). Based on our analyses, we hypothesize that the two major groups, Salix and Vetrix, might have different male and female heterogamety, respectively. One exception is Salix triandra (female heterogamety; Li et al., 2020), whose taxonomic position remains uncertain (see ‘Placement of Salix triandra’ above). We next discuss whether the sex-determination difference between the two major clades in Salix might act as a reproductive barrier.

Various mechanisms can limit plant hybridization, and both pre- and postzygotic barriers are known (Baack et al., 2015); and prezygotic reproductive barriers such as different elevation preferences (Wagner et al., 2021a) have been documented between willow species in the subgenera Vetrix and Salix.Mosseler and Papadopol (1989) concluded that common prezygotic barriers, such as ecological or spatial isolation, did not act, since a wide range of hybrids was found in natural populations within a 100-km area around Toronto. On the other hand, the same study recorded the flowering times of seven species from the subgenera Vetrix and Salix, and revealed two phenological groups: early-flowering Vetrix species and late-flowering willows belonging to Salix, suggesting that this could be a prezygotic barrier preventing interspecific hybridization between these two groups of species. Moreover, Mosseler (1989) demonstrated pollen–pistil incongruity between species of Vetrix and Salix, as successful pollinations were rare between representatives of these subgenera.

Given the impact of polyploidy on hybridization, in general, different ploidy levels produce strong crossing barriers in willows (Wagner et al., 2021a). On the other hand, homoploid crosses of species within the same subgenus usually yield viable seeds (Argus, 1974; Mosseler, 1990; Choudhary et al., 2013; Gramlich and Hörandl, 2016). For postzygotic barriers it is difficult to disentangle crossing barriers produced by different ploidy levels from those produced by putative different SDSs, and hence we focus here on crosses on the same ploidy level between the groups.

Data concerning postzygotic barriers and viability of seeds and hybrids obtained from interspecific crosses are scarce. Argus (1974) carried out pollinations between four Salix species. Among them, in crosses between tetraploid S. discolor from subg. Vetrix and S. lucida from subg. Salix (4x; Dorn, 1976), no seeds were produced. Diploid Salix exigua (sect. Longifoliae, subg. Salix, Argus, 2010) expressed neither pre- nor postzygotic barriers with diploid representatives of subg. Vetrix (Mosseler and Papadopol, 1989; Mosseler, 1989, 1990), suggesting that it may belong to this subgenus (which is consistent with the trees based on nuclear sequences in Figs 1–3 and Supplementary Data Figs S2–S4, but not with the chloroplast result in Fig. 4 and Fig. S5). Salix humboldtiana, a diploid South American species from subg. Salix (Argus, 1997), not studied here, was used in breeding experiments for developing multipurpose willows. When used as a pistillate parent, it produced progeny, though weak, when pollinated by diploid S. viminalis, S. purpurea and S. daphnoides (also from the Vetrix clade; not studied here) (Argus, 2010; Bubner et al., 2018; Förster et al., 2021). Finally, the allopolyploid origin of the C5 grade in our phylogeny could have resulted from parental lineages of the diploid Vetrix and Salix clades.

Salix fits a single-factor model of the SDS (He et al., 2021b; Renner and Müller 2021; Sanderson et al., 2021). According to Renner and Müller (2021), females in the ZW system are heterozygous for a dominant W-linked femaleness factor. In the XY system, a Y-linked specific factor in the heterozygous males dominantly suppresses female functions. A hypothetical scheme in Supplementary Data Figure S6 shows the four offspring genotype classes when a zW female is inter-crossed with an xY male. Only two of them, with W (15, dominant) and x (7, recessive), and z (15, recessive) and Y (7, dominant) have no conflict between the dominance of the SD factors on the different chromosomes, and (assuming that presence of a chromosome 7 has no effect on sex in the species whose SD locus is on chromosome 15 unless it carries a dominant factor, and vice versa) might be expected to be female and male, respectively. The zx lacks the W femaleness-specific factor, but also lacks the Y-specific factor, so it could be female. Although WY were observed in fish (Xiphophorus maculatus) (Kallman 1984), its YW individuals can be female or male. However, SDSs in plants and animals evolved in different ways (Mank, 2022). The outcome cannot be predicted with certainty for the WY of willows, too. In the reciprocal cross (zz male × xx female), all offspring are xz. Hence, even if these hybrids were viable, fertility would be very low, especially in competition with the parental population. Even if fertile, there would probably be an excess of female genotypes. Eventually, polyploidy could overcome expression bias of the diploid hybrid by gene duplication, or two locus systems could evolve; such mechanisms could have helped the allopolyploid C5 group to establish and to evolve fertile species, but this needs to be studied.

The rate of sex chromosome divergence is neither unidirectional nor correlated with time (reviewed by Mank, 2022). In willows, homomorphism of sex chromosomes does not mean that they have recently evolved (He et al., 2021b; Mank, 2022). Although turnover events are possible within each clade of Salix (Renner and Müller, 2021), it is reasonable to assume that species within each clade share the same ancestral SDS on the same chromosome (Almeida et al., 2020). SDSs of willows may have contributed to isolation of the two major groups. A comprehensive dated tree of willows and divergence time estimation of X (Z) and Y (W) with a broader sampling should further aid our understanding of the correlation between SDSs and Salix diversification.