Evidence of distant hybridization within Central Asian feather grasses (Poaceae: Stipa )

ecosystems of central Asia, we found challenging specimens of Stipa with an interesting combination of morphological characters suggesting their origination through hybridization between morphologically and phylogenetically distant species. To confirm our hypothesis, we applied a combination of classical morphological and genome-wide SNP genotyping methods. Using such an approach, we determined that the new taxon named Stipa × smelanskyi arose from crossing S. richteriana and S. drobovii and confirmed that it is an F1 hybrid. Moreover, we found a S. drobovii specimen with a minor admixture of S. richteriana loci that may indicate putative introgression events among these taxa.


Introduction
Distant hybridization refers to crosses between morphologically and phylogenetically distant species, in contrast to close hybridization that occurs between varieties of the same species, subspecies or closely related species (Chen, 2017). Distant hybridization constitutes the main facilitator of genome evolution and the formation of new species (Grant et al., 2005;Mallet, 2005Mallet, , 2007. For a long time, the importance of hybridization in evolution has been a controversial issue. However, over the past decades, thanks to applying the new molecular methods, it has been shown that hybrids are widespread in nature, essential in evolution, and can serve as a source of adaptive genetic variation. However, the role of interspecific hybridization in diversification is a matter of controversy. On the one hand, hybridization is a source of genetic variation, functional novelty, and the origin of speciation (Mayr, 1963;O'Brien, Mayr, 1991;Mallet, 1995;Allendorf et al., 2001). On the other hand, hybridization may lead to an evolutionary dead end (Mayr, 1942). The indisputable fact is that the gene flow between species via hybridization and its consequences for nature are definitely issues still worth studying for proper understanding of evolutionary processes in natural habitats.
The genus Stipa L. is represented by over 150 species distributed in warm temperate regions throughout Europe, Asia, and North Africa (Nobis et al., 2020). Currently, the taxonomic treatment of the Stipea Dumort. tribe is well developed. Morphological and molecular studies revealed 33 genera within the tribe and categorized Stipa as a monophyletic genus (Romaschenko et al., 2012;Peterson et al., 2019). Nonetheless, there is still no consensus regarding the sectional subdivision in Stipa (Roshevitz, 1934;Tzvelev, 1976Tzvelev, , 2012Freitag, 1985). Additionally, the issue of the phylogeny of the genus is greatly complicated by a large number of hybrids, many of which were previously recognized as separate species. Molecular studies have proven that hybridization and introgression processes exist in feather grasses (Nobis et al., 2019;Baiakhmetov et al., 2020Baiakhmetov et al., , 2021, and they are more common than previously believed. Kazakhstan is one of the largest and, in terms of biodiversity, richest countries in Eurasia. The flora of Kazakhstan contains 5658 species (Abdulina, 1998), 760 of which are endemic (Goloskokov, 1969). Poaceae is one of the most important families in the flora, including 538 species, 57 of which are endemic and geographically restricted to the mountain ranges of Kazakhstan (Kupriyanov et al., 2018). Within grasses, the genus Stipa is undoubtedly one of the largest and the most taxonomically difficult genus. In Central Asia, the genus Stipa comprises 72 species, including 23 of the putative hybrid origin (Nobis et al., 2020). Kazakhstan has the largest number of feather grass species among Central Asian countries and within the whole range of the genus (Bor, 1970;Freitag, 1985;Nobis et al., 2017Nobis et al., , 2020. In total, 42 species of Stipa are native to Kazakhstan, with 21 of them being endemics (Nobis et al., 2020). One of the essential characteristics of Kazakhstan's feather grasses is the richness of hybrids since the ranges of many western and eastern species overlap in its territory. To date, we recognized 13 nothospecies in Kazakhstan (Nobis et al., 2020).
During a floristic study of Kazakhstan, specimens of Stipa, which were initially proposed to represent a new nothotaxon, were found. Its morphological characters differ from other known species of Stipa and have intermediate morphological characters between S. richteriana Kar. et Kir. and S. drobovii (Tzvel.) Czerep. Here, we intend to test the hypothesis of the hybrid origin of the new taxon using a combination of classical morphological and genome-wide SNP genotyping methods.

Materials and Methods
The study is based on the herbarium materials deposited in the AA, ALTB, BM, GOET, KAS, KRA, KRAM, KUZ, LE, LECB, MW, TK, and our own collections. Herbarium acronyms follow Thiers (2021). Morphological measurements and 77 Turczaninowia 25, 2: 75-85 (2022) photos were made with a stereoscopic microscope SMZ800N (Nikon, Japan) and a scanner M1132 (HP LaserJet Pro, USA). In molecular analyses, we used samples from the localities, where S. richteriana and S. drobovii grow together with their putative hybrid as well as from those, where the two taxa grow separately from each other (Fig. 6, Supplementary Table S1). In total, 25 samples were used for the study; S. richteriana was represented by 10 individuals, S. drobovii -by 12 individuals, and the putative hybrid -by 3 individuals. Vouchers with the specimens used in the molecular analyses are preserved at TK, KRA and a personal collection of Gudkova P. D. (Supplementary Table S1).
The distribution map was visualized using ArcGIS Pro 2.7.1 (ESRI, Redlands, USA). The species distribution ranges were established based on the revision of the herbarium specimens (Fig. 6).

Molecular method
Genomic DNA was isolated from dried leaf tissues by Diversity Arrays Technology Pty Ltd. (Canberra, Australia) for the following genome complexity reduction using restriction enzymes and high-throughput polymorphism detection (Kilian et al., 2012) according to the previously reported procedures (Baiakhmetov et al., 2020(Baiakhmetov et al., , 2021. The resulted co-dominant single nucleotide polymorphism (SNP) markers were processed in the R-package dartR v.1.9.9.1 (Gruber et al., 2018) with the following parameters: (1) a scoring reproducibility of 100 %, (2) SNP loci with read depth < 5 or > 50 were removed, (3) at least 95 % loci called (the respective DNA fragment had been identified in greater than 95 % of all individuals), (4) monomorphic loci were removed, (5) SNPs that shared secondaries (had more than one sequence tag represented in the dataset) were randomly filtered out to keep only one random sequence tag.
Subsequently, four approaches were used to analyze the genetic structure of the studied taxa: (1) Unweighted Pair Group Method with Arithmetic Mean (UPGMA), (2) STRUCTURE analysis, (3) Principal Component Analysis (PCA), and (4) NewHybrids analysis. Firstly, the UPGMA cluster analysis based on a Hamming distance matrix with 1000 bootstrap replicates was performed in the R-package poppr v.2.9.1 (Gruber et al., 2021). The final UPGMA tree was visualized via iTOL v.6.3.2 (Letunic, Bork, 2021). Next, the genetic structure was investigated using STRUCTURE v.2.3.4 (Pritchard et al., 2000) via StrAuto v.1.0 (Chhatre, Emerson, 2017). A number of clusters (K-values) ranging from 1 to 5 were tested using ten replicate runs per dataset with a burn-in of 10,000 iterations followed by 100,000 Markov chain Monte Carlo (MCMC) iterations. The optimal K value was identified based on Evanno's method of ΔK statistics (Evanno et al., 2005) as implemented in Structure Harvester (Earl, Vonholdt, 2012). The R package pophelperShiny v.2.1.1 (Francis, 2017) was used to visualize the output matrix. We applied the threshold of 0.10 < q < 0.90 for the assessment of hybridization with q-values > 0.9 being pure species and 0.45 < q < 0.55 being F1 hybrids, while first-and second-generation backcrosses with one parent were considered at q-values 0.25 and 0.125, respectively. Then, PCA on a Euclidean distance matrix was performed using the R-package dartR and visualized with ggplot2 v.3.3.5 (Wickham, 2009). Finally, to assign individuals to a genetic category (pure, hybrid F1 or F2, backcross hybrid) we performed an analysis using NewHybrids v.1.1 (Anderson, Thompson, 2002) via the R-package dartR. Due to limitations of NewHybrids, only 200 loci were used per run. In the beginning, the first 200 loci ranked on information content were chosen via dartR. Then, we used ten different 200-SNP subsets selected in random. The Jeffreys prior was used for θ and π and a burn-in of 10,000 MCMC generations followed by 10,000 post burn-in sweeps. The resulting posterior probabilities were calculated based on 11 runs and a probability threshold > 0.9 was set for the assignment into a genetic category. The calculated posterior probabilities for the assigned categories were visualized using the R-package pophelperShiny.

Morphological analysis
A total of 25 specimens were measured across the 19 quantitative and 6 qualitative morphological characters used in our previous studies (Baiakhmetov et al., 2020(Baiakhmetov et al., , 2021Nobis et al., 2020; Table 1).
We used the R-package FactoMineR v.2.3 (Husson et al., 2014) to carry out a Factor Analysis of Mixed Data (FAMD; Pagès et al., 2021) aimed at characterizing the variation within and among groups of taxa without a priori taxonomic classification and at extracting the variables that best identify them. Cattell's scree test (Cattell, 1966) determined the number of principal components used in the analysis. R-packages factoextra v.1.0.6 (Kassambara, Mundt, 2021) and plotly v.4.9.2 (Sievert et al., 2020) were applied to visualize the first two and the first three principal components, respectively.
In addition, to estimate distributional relationships between each response variable and the Stipa × smelanskyi a new nothospecies taxa under consideration, notch plots and interactive box plots were constructed using R-packages ggplot2 v.3.3.0 and plotly v.4.9.2, respectively. The notched box plots visualize a confidence interval around the median that is normally based on the median ± 1.57 × interquartile range/square root of n. Such graphical method for data analysis proves that, if the notches of the two boxes do not overlap, there is "strong evidence" (95 % confidence) that their medians differ.

Results
The genetic structure of the studied specimens was estimated with 6316 SNP markers obtained using the DArTseq technique. Firstly, the analysis of genetic clustering with an unweighted pair group method (UPGMA) using arithmetic average revealed three major clades corresponding to species S. richteriana and S. drobovii, while the third one was represented by the hybrid S. drobovii × S. richteriana (Fig. 1). A congruent result was achieved by the STRUCTURE analysis with the best estimate of K = 2. The first cluster was represented by pure individuals of S. richteriana, while the second one was comprised by pure samples of S. drobovii. The remaining specimens were admixed between the pure clusters in the 50/50 ratio and were constituted by three individuals of S. drobovii × S. richteriana. Additionally, one individual of S. drobovii (002721a) had a minor proportion (0.02) of the first cluster representing S. richteriana. A similar result was obtained with a principal component analysis (PCA). The first two axes explained 90.4 % and 3.8 % of the total genetic divergence within the studied taxa, respectively. The pure individuals were grouped into two markedly differentiated groups, while F1 hybrids had intermediate positions between their parental species. Noteworthy, in accordance with the second axis, two specimens of S. richteriana (0452466 and 002707) were slightly distant to the remaining pure samples of the species. Additionally, to assess if a more complex pattern of hybridization is present within the admixed individuals, we performed an analysis in NewHybrids (Fig. 2). The Table 1 Morphological characters used in the present study  (Fig. 3, the interactive plot available at https://plot.ly/~eugenebayahmetov/44/). The first three dimensions explained 46.69 %, 12.62 %, and 8.17 % of the total variability, respectively. The first dimension is composed, in order of descending contribution, by the quantitative variables LG, LHCol I, LAnt, LHS, LHDor, WAwn, LBasCal, and qualitative variables AdSGL, AdSVL, AG, CharHAnt (Supplementary Table S2, for the character abbreviations see Table 1). The second dimension is composed, in order of descending contribution, by the quantitative variables LGL, LVL, DVL, DGL, and the qualitative variable AdSGL and AdSVL (Supplementary Table S2). The third dimension is composed, in order of descending contribution, by the quantitative variables LigV, DVL, DGL, LigGH, LVL, LGL, LAwn, and the qualitative variable AbSVL, AbSGL, AdSVL (Supplementary Table S2). The two-dimensional plot revealed a clear dispersal of all analyzed taxa (Fig. 3a).  (Table 2, Supplementary Fig. S1, S2, S4).
At the same time, Stipa × smelanskyi and S. richteriana share five characters with no significant differences between their means but differ with S. drobovii in the length of ligules of vegetative shoots, the length of hairs on the ligules of the vegetative shoots, the length of the callus base, the length of hairs on the dorsal line on the lemma, and the number of the awn geniculation. However, the three characters, i.e., the length of hairs on the top of the lemma, the length of the awn, and the type of hairs distributed on the lemma have no significant value in differentiation of Stipa × smelanskyi and S. drobovii but differ between Stipa × smelanskyi and S. richteriana (Table 2, Supplementary Fig. S1, S2, S3, S4).

Discussion
Recently, numerous hybrids, i.e., intersection hybrids, were described within the genus Stipa (Nobis, 2014;Nobis et al., 2019Nobis et al., , 2020. Some of them were confirmed with molecular analyses. For instance, it was evidenced that Stipa × heptapotamica Golosk. arose from hybridization of genetically closely related species S. richteriana and S. lessingiana Trin. et Rupr. (Nobis et al., 2019). Following the last sectional division of the genus Stipa, these species belong to sections Leiostipa Dumort. and Barbata Junge, respectively (Tzvelev, 2012). Another example is S. × lazkovii M. Nobis et A. Nowak that is an F1 hybrid of two genetically distant species, S. krylovii Roshev. and S. bungeana Trin. Here, we revealed a hybrid, Stipa × smelanskyi originated from morphologically and genetically distant species S. richteriana and S. drobovii that currently represent sections Leiostipa and Smirnovia Tzvelev, respectively.
The current study is the first molecular evidence for hybridization between species from the sections. Although previously we recorded the probability of hybridization between these sections in the north of Lake Issyk-Kul (Baiakhmetov et al., 2020), here we reliably confirm this assumption based on morphology and molecular data.   The molecular analyses confirmed that Stipa × smelanskyi is the F1 hybrid of S. richteriana and S. drobovii. Moreover, we found that one sample of S. drobovii has a minor admixture of S. richteriana that may indicate putative introgression events among these taxa. Besides, our study revealed two samples of S. richteriana (0452466 and 002707) deviating from the remaining specimens of the taxon, both morphologically and molecularly. This is connected to a high level of the genetic diversity of the taxon that was also shown in our other work (Nobis et al., 2019).

Diagnosis
Stipa × smelanskyi is similar to S. richteriana but differs in having 7 distinct lines of hairs on the lemma vs. lemma throughout pilose, hairs on the lower segment of the awn 0.7-1 mm vs. 0.2-0.3 mm long, hairs on the seta 1.5-2 mm vs. 0.3-0.9 mm long, respectively. Stipa × smelanskyi differs also from S. drobovii in having bigeniculate vs. unigeniculate awns, hairs on the dorsal part of the callus straight vs. falcate hairs, hairs on the lower segment of the awn 0.7-1 mm vs. (1.8)2-2.5 mm long, and on the seta 1.5-2 mm vs. 2.6-5.5 mm long, respectively.

Description
Plant perennial, densely tufted, with a few culms and numerous vegetative shoots; culms 40-75 cm, 3-noded (Fig. 4a), glabrous at the nodes and shortly pubescent below them (Fig. 4e). Leaves of vegetative shoots: sheaths glabrous or shortly pubescent with hairs 0.05 mm long, ciliate at margins with hairs 0.15-0.3 mm long; ligules rounded up to 0.25 mm long, densely ciliate with hairs up to 0.2-0.4 mm long (Fig. 4i); blades convolute, up to 20-30 cm long, 0.7-0.8 mm in diameter, abaxial surface glabrous, adaxial surface densely pubescent with mixture of short prickles and long hairs 0.1-0.3 mm (Fig. 4l). Cauline leaves: sheaths glabrous, usually shorter than internodes; ligules 0.25-0.4 mm, obtuse and ciliate at margins; blades glabrous, middle leaves up to 10 cm long, upper 5-8 cm long. Panicle 10-20 cm long contracted, branches erect. Glumes subequal, 18-23 mm long, narrowly lanceolate, tapering into long hyaline apex. Anthecium 5.75-8.25 mm long (with callus), ca. 0.7-1 mm wide; callus 0.5-1.1 mm long, densely long-pilose on ventral and dorsal surfaces, callus base acute, cuneate, scar elliptic to circular. Lemma with 7 lines of ascending hairs, hairs up to 0.5 mm, all of hairs line reach the top of lemma; top of lemma with hairs at apex up to 0.2-1.3 mm long (Fig. 4g). Palea equal to lemma in length. Awn 65-75 mm long, bigeniculate; lower segment of column 22-24 mm, twisted, ca. 0.5 mm wide near base, 82 Gudkova P. D. et al. Stipa × smelanskyi a new nothospecies shortly pilose, with hairs 0.7-1 mm, in the upper part gradually increasing in length toward second geniculation; seta 40-60 mm long, falcate, hairs in lower part of seta 1.5-1.8 mm, gradually decreasing in length toward apex (Fig. 4c).  Habitat and distribution Stipa × smelanskyi grows in open areas, in the steppes, semideserts, dry grasslands and rocky grasslands, at an elevation between 500 and 700 m. It is known only from Zhanaarkinsky district, in central Kazakhstan, however, it can be found within the overlapping ranges of the parent species (Fig. 6).

Phenology
Flowering and fruiting April to June.

Etymology
The specific epithet of the species is assigned in honors of the eminent Russian scientist, ecologist, and nature conservationist, Ilya Eduardovich Smelansky.