Molecular and cytological evidences denied the immediate-hybrid hypothesis for Saxifraga yuparensis (sect. Bronchiales , Saxifragaceae) endemic to Mt. Yubari in Hokkaido, northern Japan

An alpine plant Saxifraga yuparensis is endemic to a scree consisting of greenschist of Mt. Yubari in Hokkaido, Japan and it has been proposed as an immediate hybrid derived from two species of the same section Bronchiales based on morphological intermediacy: namely S. nishidae , a diploid species endemic to a nearby cliff composed of greenschist and tetraploid S. rebunshirensis comparatively broadly distributed in Japan and Russian Far East. Saxifraga yuparensis is red-listed and it is crucial for conservation planning to clarify whether this is an immediate hybrid and lacks a unique gene pool. The im-mediate-hybrid hypothesis was tested by molecular and cytological data. In nuclear ribosomal and chloroplast DNA trees based on maximum parsimony and Bayesian criteria, S. yuparensis and S. rebunshirensis formed a clade with several other congeners while S. nishidae formed another distinct clade. Genome-wide SNP data clearly separated these three species in principal coordinate space, placing S. yuparensis not in-between of S. rebunshirensis and S. nishidae . Chromosome observation indicated that S. yuparensis is tetraploid, not triploid directly derived from diploid-tetraploid crossing. Additionally, observation of herbarium specimens revealed that leaf apex shape of S. yuparensis fell within the variation of S. rebunshirensis . These results indicate that S. yuparensis is not an immediate hybrid of S. rebunshirensis and S. nishidae but a distinct lineage and an extremely narrow endemic species, that deserves for intensive conservation. hybrid, Japan, MIG-seq, Russia, Saxifraga , tetraploid


Introduction
The genus Saxifraga L. (1753: 398) (Saxifragaceae, e.g., Soltis et al. 2001) encompasses more than 440 species and is composed of at least 13 sections and nine subsections (Tkach et al. 2015).One of the Saxifraga sections, sect.Bronchiales DeChaine E.G. (2014: 27) was subdivided from sect.Trachyphyllum Koch W.D.J. (1836: 270) based on molecular phylogenetic analyses (DeChaine et al. 2013, DeChaine 2014).This section has approximately 20 species (DeChaine et al. 2013, DeChaine 2014) mostly in northern high latitudes, and among them, three species are distributed in Hokkaido, northern Japan: S. yuparensis Nosaka (1974: 149), S. rebunshirensis Sipliv. (1971: 155), and S. nishidae  Miyabe & Kudô (1917: 170) (Fig. 1 & 2).Saxifraga rebunshirensis is broadly distributed in Japan (Hokkaido and high mountains in Honshu) and Russian Far East (Sakhalin and the Kuril Islands), while S. yuparensis is endemic to a scree consisting of greenschist at Mt. Yubari in Hokkaido and S. nishidae is endemic to a cliff composed of greenschist at Mt. Yubari and a cliff of schalstein at Mt. Ashibetsu, Hokkaido (Shimizu 1983, Charkevicz 1989, Sato 2007, Nakagawa & Sato 2015).Mountains composed of diabasic rocks (including greenschist, shalestein, etc.) are known to house plants MOLECULAR AND CYTOLOGICAL EVIDENCES Phytotaxa 373 (1) © 2018 Magnolia Press • 55 adapted to infertile soils and fragile geological features, as is the case in several endemics in Mt.Yubari (Watanabe 1971, Nosaka 1974).Diabasic plants are typically narrow endemics and grow on unstable soils (Watanabe 1971), and thereby include many endangered species.The number of individuals and patches of S. yuparensis are decreasing due to collapse of rocks and collection for horticultural purposes, and it is designated as an endangered species: CR (critically endangered) at the national level (Ministry of the environment of Japan 2017) and En (Endangered) at the provincial level (Hokkaido 2001).
The diagnostic characters of these three species are the color of petal spots, presence/absence of a conspicuous claw of petal base, and the shape of leaf apex.Saxifraga yuparensis has yellow petal spots, petal claws, and slightly tricuspidate leaves; S. rebunshirensis has yellow and red (irregularly only yellow) petal spots, claw-less petals, and cuspidate leaves; S. nishidae has yellow and red petal spots, claw-less petals, and tricuspidate leaves (Miyabe & Kudô 1917, Hara 1952, Murata 1961, Nosaka 1974, Nosaka 1980, Shimizu 1983, Ohba 1999, Nakagawa & Sato 2015).However, several authors considered that S. yuparensis is an immediate hybrid derived from crossing between S. rebunshirensis and S. nishidae based on morphological intermediacy in the shape of leaf apex (Toyokuni 1988, Nishikawa et al.1992, Iwatsuki & Kato 1994, Umezawa 2004, Shimizu et al. 2014, Takahashi 2015a).Additionally, seeds of S. yuparensis rarely germinate (Shimizu 1983) and this suggests that the species is an infertile triploid hybrid, because S. nishidae is diploid (2n = 26, Funamoto & Nakamura 1993) while S. rebunshirensis is tetraploid (2n = 48 using samples with no locality information, Sakai, 1935; 2n = 50 with samples from Nagano of central Japan, Funamoto & Nakamura 1990), although the chromosome number of S. yuparensis has not been reported.The hybridorigin hypothesis, however, has not been tested by molecular analyses.On the other hand, S. yuparensis is sometimes treated as a synonym or variety of S. rebunshirensis (Shimizu 1983, DeChaine 2014, Takahashi 2015a, Okuyama 2016; note the last one is the latest Japanese flora).For planning the conservation of S. yuparensis, it is crucial to clarify whether this is an immediate hybrid and lacks a unique gene pool.In a preceding molecular phylogenetic study of Saxifraga, the species relationships were partially revealed, but the relationship among the three species remained unclear (DeChaine et al. 2013).The aim of this study is to elucidate the species relationships and test the immediate-hybridity of S. yuparensis using molecular analyses and cytological observation.In molecular analyses, DNA sequencing and genome-wide single nucleotide polymorphism (SNP) detection were employed.Collaterally, the diagnostic character of leaf apex shape was investigated using specimens of the three species to reevaluate the morphological intermediacy of S. yuparensis.

Molecular analyses
Taxon sampling Molecular phylogenetic analyses were conducted using samples of our collection and DNA sequence data from GenBank (Table 1).Saxifraga yuparensis grows on a rocky scree of Mt.Yubari and has only two patches (0.3 x 0.2 m and 0.5 x 1.5 m), being located only ca. 5 m away from each other.One sample from each of the two patches (two samples in total) was collected to avoid collecting the same clones because this species can propagate vegetatively (authors' observation).DeChaine et al. (2013) sequenced a sample of "S.nishidae" (RBGE-E00295524) but we identified this specimen as S. yuparensis based on leaf apex morphology and found that the specimen was confused and mislabeled with genuine S. nishidae (RBGE-E00295525) collected on the same day on Mt.Yubari.The data of "S.nishidae" (RBGE-E00295524) was used for S. yuparensis.Saxifraga nishidae has more than 100 patches on a nearby rocky cliff of greenschist of Mt.Yubari.Six samples were collected from different patches with at least 5 m interval.In Mt.Yubari, S. rebunshirensis (Ken Sato 85.0207 in SAPS, Appendix) had been collected more than 30 years ago from Rosoku-iwa, that is several kilometers away from the rocky place where S. yuparensis and S. nishidae grow; but this time we could not collect the species there.We did not conduct DNA extraction/PCR with the herbarium specimen because the specimen was old and we were afraid to damage the valuable specimen in vain.Instead, two samples each from two localities in Hokkaido, Rebun Island and Mt.Nishi-Kumaneshiri, were used.Also, two plants from Mt. Hakuba in central Japan were hired.GenBank data for S. rebunshirensis from Moneron Island to the south of Sakhalin (UBC-V164570) were also used.To elucidate the phylogenetic positions of these three species, we incorporated allied 16 species of sect.Bronchiales (DeChaine et al. 2013) in the analyses.In Russian Far East, we collected S. ascoldica Sipliv.(1971: 156) in Primorsky Krai and S. cherlerioides D. Don (1822: 382) and S. funstonii (Small) Fedde (1905: 613) in Kamchatka.For the other species, DNA sequence data were obtained from GenBank.Distribution of the 19 species covered the whole distribution of sect.Bronchiales.For outgroups, two species of sect.Gymnopera D. Don (1822: 343) and one species of sect.Cymbalaria Griseb.(1843: 336) were selected (Table 1), following the result of Tkach et al. (2015).

DNA extraction and sequencing
Total genomic DNA was extracted from fresh leaves using the cetyltrimethyl ammonium bromide (CTAB) method (Doyle & Doyle, 1987) with some modifications.Preceding phylogenetic studies on Saxifraga employed internal transcribed spacer (ITS) of nuclear ribosomal DNA (nrDNA) and trnL-F intergenic spacer of chloroplast DNA (cpDNA) (DeChaine et al. 2013, Tkach et al. 2015).We used primers following these preceding studies: trnL(UAA) and trnF(GAA) for trnL-F region (Taberlet et al. 1991), and 17SE and 26SE for the entire ITS1, 5.8S and ITS2 region (Sun et al. 1994).PCR was performed in 25 μl total volume with the following reagents: about 10 ng of genomic DNA, 1 unit of Taq DNA polymerase master mix (Ampliqon, Rødovre, Denmark), 0.4 μM of each primer, and 4 % DMSO.After an initial heating step at 94°C for 3 min, samples were incubated for 25 cycles of 94°C for 1 min, 60°C (ITS) or 52°C (trnL-F) for 1 min, 72°C for 2 min, with final extension at 72°C for 5 min.The cycle sequencing reaction was carried out with a Big Dye Terminator Cycle Sequencing Kit ver.3.1 (Applied Biosystems, Foster City, CA, USA) with the same primers used in the PCR.Direct sequencing was performed on an ABI Prism 3130 DNA analyzer (Applied Biosystems).The sequence data were deposited in DDBJ (DNA Data Bank of Japan) database (2017, Table 1).
In the Bayesian phylogenetic analysis, MrModeltest 2.3 (Nylander 2008) was used to estimate the appropriate evolutionary model of nucleotide substitution based on the Akaike Information Criterion (AIC).Based on the model selected, two separate runs of Metropolis coupled Markov chain Monte Carlo (MCMCMC) analyses were performed, each with a random starting tree and four chains (one cold and three heated).The MCMCMC length was one million generations, and the chain was sampled every one hundredth generation from the cold chain.The mixing and convergence of the MCMC chains of the two runs was assessed by inspection of the trace plots of parameters using Tracer ver.1.6 (Rambaut et al. 2013).The first 250 sample trees (2.5% of the total 10,000 sample trees) were discarded as burn-in.After the burn-in, the effective sample sizes (ESS) of all parameters were more than 500, indicating that MOLECULAR AND CYTOLOGICAL EVIDENCES  the analyses sampled the posterior distributions of each parameter satisfactorily, and the values of Average Standard Deviation of Split Frequency (ASDSF) were below 0.008.The 50% majority rule consensus tree and Bayesian posterior probabilities (PP) of all the post-burn-in trees was generated using FigTree ver.1.4.2(Rambaut 2014).
In the MP phylogenetic analysis, indels were treated as missing data.The characters were treated as unordered, and the character transformations were weighted equally.The branch collapse option was set to collapse at a minimum length of zero.A heuristic parsimony search was performed with 1000 replicates of random additions of sequences with ACCTRAN character optimization, tree bisection-reconnection (TBR) branch swapping, and MULTREES and STEEPEST DESCENT options on.Statistical support for each clade was assessed by bootstrap analysis (Felsenstein 1985).One thousand replicates of heuristic searches, with the TBR branch swapping switched on and MULTREES options off, were performed to calculate bootstrap percentages (BP).

Genome-wide SNP analysis
To test the immediate hybridity and species relationships based on multiple independent nuclear loci, genome-wide SNP detection was conducted with the focal three species.The samples of our own collection used in the sequencing analysis were employed.Preceding studies (Rutledge et al. 2015, Karimi et al. 2016) based on genome-wide SNP data (observed and simulated data) successfully indicated that hybrids were plotted at the intermediate between the parental species in a principal coordinate analysis (PCA) (e.g., F1 and F2 individuals occurred at intermediate between the two parental species, and backcross individuals occurred between the F1 and the one parental species; Rutledge et al. 2015).For SNP detection, we employed a method MIG-seq [i.e., multiplexed inter simple sequence repeats (ISSRs) Genotyping by sequencing; Suyama & Matsuki 2015].MIG-seq is a type of reduced representation sequencing, where ISSRs are amplified by multiplex PCR.Thereby the method is more applicable to low-quality genomic DNA than RADseq.Most of SNPs obtained by MIG-seq are putatively selectively neutral because these are immediately adjacent to two repeat regions (Takahashi et al. 2016).Highly reduced representation libraries were constructed following Suyama & Matsuki (2015) and sequenced on an Illumina MiSeq sequencer (Illumina, San Diego, CA, USA).Both ends of the fragments were read by paired-end sequencing.
Sequences of reads were quality-filtered using FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/),with the settings of q = 30 and p = 40.The quality-filtered reads were then bundled together to form putative loci (stacks) with the software Stacks v1.15 (Catchen et al. 2013) using the 'ustacks' option, with the settings of minimum depth of coverage required (m) = 20, maximum distance between stacks (M) = 2, and the deleveraging (d) and removal (r) algorithms enabled.Using the 'cstacks' option, a catalog was created for all possible loci and alleles with the parameter 'number of allowed mismatches between samples (n)' = 4.The stacks created by 'ustacks' were then searched against the catalogue using the 'sstacks' option.The data set of all the samples was considered as a single population and SNPs were retrieved using the 'populations' option, with the 'write_single_snp' option to select only the first SNP per locus to avoid the linkage between SNPs.Only loci that were shared among at least 75 % of the samples were considered.If a locus was not recovered in a sample, the genotype of the sample was treated as missing data.PCA of SNP genotypes was conducted using GENODIVE (Meirmans & Van Tienderen 2004) with the default parameter settings.

Morphological observation
A diagnostic character of leaf apex shape was observed in the three species.Typically, S. rebunshirensis has cuspidate leaves while S. nishidae has tricuspidate leaves; and S. yuparensis has slightly tricuspidate leaves (Shimizu 1983 Wakabayashi 2001, Shimizu 2014, Okuyama 2016).Herbarium specimens of the three species (Appendix) were examined in the herbaria of Hokkaido University Museum (SAPS), Botanic Garden, Hokkaido University (SAPT) and the National Museum of Nature and Science (TNS) to elucidate the variability of the diagnostic character especially within S. rebunshirensis, that has large species range.

Phylogenetic relationships based on ITS
The aligned length of the ITS data was 822 bp.In S. rebunshirensis, five sequence types were recognized among the seven samples, and S. nishidae had two sequence types in the six samples.The two samples of S. yuparensis of our collection (SAPT-042183, 042184) showed one type of sequence with a double-peak signal (adenin and guanin) at one site (position 423 in the registered sequences) in the sequence chromatogram, while all the samples of S. rebunshirensis and S. nishidae had only adenin but not guanin at the same site.This site was coded as R in the two S. yuparensis samples.In the Bayesian analysis, the AIC selected the SYM+I+G model.The 50% majority rule consensus tree of all the post-burn-in trees is depicted (Fig 3).In the MP analysis, 201 nucleotide substitutions were found in 57 variable sites and 114 sites were parsimony informative among them.15,104 equally parsimonious trees of 314 steps were obtained with a consistency index (CI) = 0.78, a retention index (RI) = 0.91, and a resealed consistency index (RC) = 0.71.The topology of the strict consensus tree was the same as that of the Bayesian tree, and therefore BPs of the MP analysis are plotted on the Bayesian tree (Fig 3).In the following, only clades with PP ≥ 0.95 and/or BP ≥ 70% were considered adequately supported.The phylogenetic tree had moderate resolution and supported monophyly for several species, including S. nishidae (PP = 1.00 / BP = 100).On the other hand, monophyly was not supported for each S. yuparensis and S. rebunshirensis; the two species were included in a clade (1.00 / 55.4) with a portion of samples (three samples) of S. funstonii and S. caulescens Sipliv. (1971: 151).This clade was nested within a clade with the other five samples of S. funstonii and one S. codyana (1.00 / 66.3).Thereby, S. yuparensis plus S. rebunshirensis were phylogenetically distinct from S. nishidae.

Phylogenetic relationships based on trnL-F
The aligned length of the trnL-F data was 446 bp.All the samples of S. rebunshirensis and S. yuparensis had the same sequence, that was shared with other 11 species in the same clade (below).In the six samples of S. nishidae, only one type of sequence was recovered.In the Bayesian analysis, the AIC selected the GTR model.The 50% majority rule consensus tree of all the post-burn-in trees is depicted (Fig 4).In the MP analysis, 28 nucleotide substitutions were found in 12 variable sites and 13 sites were parsimony informative among them.574,318 equally parsimonious trees of 30 steps were obtained with CI = 1.The topology of the strict consensus tree was the same as that of the Bayesian tree, and therefore BPs of the MP analysis are plotted on the Bayesian tree (Fig 4).While the phylogenetic tree largely had low resolution, but nevertheless S. nishidae formed a well-supported clade (1.00 / 97.2), and this clade was sister to a portion of samples (three samples) of S. austromontana (0.97 / 63.0).In contrast, each S. yuparensis and S. rebunshirensis was not monophyletic; these two species were recovered in the same polyphyletic clade (0.98 / 61.1) with other 11 species.

Species relationships based on genome-wide SNPs
The resultant data set comprised 699 SNP loci.The result of the PCA based on the SNP genotyping data is shown (Fig. 5).The first and second principal components explained 51.6 % and 11.5 % of the variances.The SNP genotyping data clearly separated the three species in the principal coordinate space; the plots of S. nishidae were separated from S. yuparensis and S. rebunshirensis along the first principal coordinate, and those of S. yuparensis and S. rebunshirensis were distinguished along the second principal coordinate.Note that S. yuparensis was not placed in-between of S. rebunshirensis and S. nishidae.

Chromosome number
The chromosome number of S. yuparensis was first determined as 2n=48 (Fig 6-A & D).In S. rebunshirensis the chromosome number was 2n=48 in all the three individuals examined (Fig 6 -B & E) and this is the same as reported by Sakai (1935) but different from the number reported by Funamoto & Nakamura (1990) (2n = 50).The chromosome number of S. nishidae was 2n=26 in all the four individuals (Fig 6 -C & F), and that was the same as reported by Funamoto & Nakamura (1993).

Discussion
In the molecular phylogenetic analyses based on ITS and cpDNA, S. yuparensis was not monophyletic and recovered in the same clade with S. rebunshirensis plus several other species (S. funstonii and S. caulescens in ITS; 11 species in trnL-F).On the other hand, S. nishidae formed a well-supported monophyletic clade in both ITS and cpDNA, and this clade was distinct from the clades including S. yuparensis and S. rebunshirensis.If S. yuparensis were a hybrid derivative of the two species, two distinct sequences from each of the S. yuparensis samples would be recovered in both the clades of S. rebunshirensis and S. nishidae in the ITS phylogeny, although it cannot be excluded that the ITS allele of one parental species (here S. rebunshirensis) became dominant by concerted evolution via unequal crossingover and gene conversion through time (given that S. yuparensis is not an immediate but old hybrid).However, the PCA based on the SNP genotyping data clearly separated the three species and S. yuparensis samples were not placed in-between of S. rebunshirensis and S. nishidae.This result, being based on multiple independent nuclear loci, also negate the hybrid hypothesis of S. yuparensis.The result also likely negates the scenario that ITS allele of S. rebunshirensis became dominant by concerted evolution in S. yuparensis.Chromosome observation revealed that S. yuparensis is tetraploid.If it were an immediate hybrid derived from diploid-tetraploid crossing, it would be triploid.Thus, the chromosome data negate, at least, the immediate hybridity.Considered in a comprehensive way, all the results based on ITS and cpDNA sequences, genome-wide SNPs, and chromosome count negate the hypothesis that S. yuparensis is derived from hybridization between S. rebunshirensis and S. nishidae.Concerning the infertility of S. yuparensis, Shimizu (1983) reported that the seeds of S. yuparensis did not germinate, and we have not observed mature fruits in the wild population.The infertility of S. yuparensis is probably caused by self-incompatibility, not by triploidy.Self-incompatibility in Saxifraga was reported for species of the same section Bronchiales (Goertzen 1996).Currently there are two patches of S. yuparensis at Mt. Yubari but these patches were originally one patch, that were subsequently fragmented by rock debris (authors' observation), and the two patches likely contain very closely related plants such as those of the same sibling that are incompatible.
Concerning the taxonomic treatment of S. yuparensis, the SNP analysis suggests that S. yuparensis is a lineage distinct from S. rebunshirensis.Note, however, that S. rebunshirensis itself has taxonomic problems and is sometimes treated as a subspecies or variety of S. cherlerioides, S. bronchialis L. (1753: 400), or S. funstonii (Engler et al. 1919, Hara 1937, Hara 1952, Toyokuni 1975, Takahashi 2015b).The present phylogenetic analyses are not conclusive about the species relationships because monophyly was not supported for each of these species and because the trees lacked necessary resolution.The SNP analysis did not include the other species treated in the phylogenetic analyses, specifically S. funstonii and S. caulescens clustered with S. yuparensis and S. rebunshirensis in the ITS tree (and in the cpDNA tree with the other nine species, although it showed lower resolution).Full revision of their taxonomy needs more intensive analysis of the allied species using genome-wide SNPs, but it is beyond the scope of the present study (i.e., test on the hybrid hypothesis).When we follow the idea to treat S. rebunshirensis as an independent species (Siplivinsky 1971;DeChaine 2014), what is the diagnostic character to separate it from S. yuparensis?Slightly tricuspidate leaf is a character based on which S. yuparensis has been distinguished from S. rebunshirensis.In the morphological observations of specimens, it was revealed that S. rebunshirensis from Hokkaido, Honshu, and Sakhalin sometimes had slightly tricuspidate leaves.The leaf apex shape of S. yuparensis fell within the variation of S. rebunshirensis and thereby is not a diagnostic character.In addition, Barkalov (2009) conducted morphological comparison between S. yuparensis from Mt. Yubari and S. rebunshirensis samples from the southern Kuril Islands, where S. yuparensis is not reported and which should be treated as S. rebunshirensis (Charkevics 1989; Takahashi 2015b); he recognized no morphological difference.According to literature, another diagnostic character of S. yuparensis and S. rebunshirensis is a conspicuous claw of petal base in the former species, that is lacked in the latter species (Nosaka 1980;Shimizu et al. 2014), while Barkalov (2009) did not mention the character.Although it was difficult to examine the character with dried herbarium specimens, by observing living plants in wild and botanic gardens, we also recognized the difference in the two species (authors' observations, Appendix).
In conclusion, S. yuparensis is not an immediate hybrid of S. rebunshirensis and S. nishidae but a distinct lineage, as was revealed by the DNA sequence, genome-wide SNP, and cytological evidences.This is an extremely narrow endemic species, that is morphologically distinguished from the congeners in Hokkaido by having a claw of petal base, and deserves for intensive conservation.

FIGURE 2 .
FIGURE 2. Maps of East Asia (A) showing the species range of Saxifraga rebunshirensis and of Hokkaido (B) showing the range of S. yuparensis (red; square), S. rebunshirensis (blue; circle), and S. nishidae (green; triangle).(the map image from Google Earth).The figure is available in color online.

FIGURE 3 .
FIGURE 3. Bayesian majority-rule consensus tree of Saxifraga sect.Bronchiales based on ITS.The numerals on branches are Bayesian posterior probabilities (PP: upper) and bootstrap percentages (BP: lower) in the MP analysis.The figure is available in color online.

FIGURE 4 .
FIGURE 4. Bayesian majority-rule consensus tree of Saxifraga sect.Bronchiales based on trnL-F.The numerals on branches are Bayesian posterior probabilities (PP: upper) and bootstrap percentages (BP: lower) in the MP analysis.The figure is available in color online.

FIGURE 6 .
FIGURE 6. Photomicrographs of somatic metaphase chromosomes of Saxifraga yuparensis (A), S. rebunshirensis (B) and S. nishidae (C).D, E and F are explanatory drawings of A, B and C, respectively.Scale bar represents 5 μm.

FIGURE 7 .
FIGURE 7. Specimens of Saxifraga yuparensis (A & D), S. rebunshirensis (B, E & F), and S. nishidae (C & G).Scale bars are 5 mm for D-G.The figure is available in color online.

TABLE 1 .
Samples used for ITS and cpDNA sequencing analyses.N is the number of sequenced individuals.