Ecological speciation of Japanese hedgehog mushroom: Hydnum subalpinum sp. nov. is distinguished from its sister species H. repando-orientale by means of integrative taxonomy

Hydnum repando-orientale is an East Asian species closely related to H. boreorepandum and H. repandum; all three species produce edible mushrooms. We identified two ecological groups of H. repando-orientale in Japan: a temperate group occurring in Fagaceae-dominated forest at < 1200 m a.s.l. (ROF) and a subalpine group occurring in coniferous forest in highland at > 1900 m a.s.l. (ROC). We re-examined the taxonomy of the two ecological groups of H. repando-orientale using integrative approaches. Phylogenies of the two ecological groups and other related species were inferred from the internal transcribed spacer (ITS) and gene portions encoding the large subunit of nc rRNA (LSU), translation elongation factor-1 alpha (TEF1), RNA polymerase II largest subunit (RPB1), and RNA polymerase II second-largest subunit (RPB2). The concatenated phylogenetic tree separated the two ecological groups into well-supported sister clades. Also, species delimitations based on the topological congruence (GCPSR) and coalescent models (GMYC and BP&P) supported to separate the two ecological groups. Morphological analysis showed that ROC specimens had significantly larger basidiospores, compared with ROF specimens. Mon-mon mating tests using six ROF, three ROC, and three H. boreorepandum strains each showed independent incompatible groups, whereas one ROC strain showed compatibility with both ROC and ROF populations. Based on these results, we defined the ROC group as a new species, H. subalpinum. Because H. repando-orientale and H. subalpinum have smaller genetic divergence in nc rDNA and maintain slight sexual compatibility, they may have recently speciated in East Asia.

Hydnum repandum has been also reported from Japan (Kawamura 1913(Kawamura , 1929(Kawamura , 1954Yasuda 1913;Asahina 1939;Ito 1955;Imazeki and Hongo 1957;Kikuhara 1987;Yanaga 2015;Yanaga et al. 2015;Sugawara et al. 2019 (Yanaga 2015;Yanaga et al. 2015;Niskanen et al. 2018;Sugawara et al. 2022a)-of these, only H. repandoorientale has not been re-classified into another subgenus. Currently, only two species of subg. Hydnum (i.e., H. boreorepandum and H. repando-orientale) have been proven to be distributed in Japan (Niskanen et al. 2018;Sugawara et al. 2022a). Hydnum boreorepandum and H. repando-orientale were recognized as sister species of H. repandum in the phylogeny inferred from the nc rDNA internal transcribed spacer (ITS) sequences (Niskanen et al. 2018); these three species were distinguished by their geographical distribution patterns in Eurasia and by forest habitats. Hydnum repandum occurs in Europe, H. repandoorientale occurs in East Asia, and H. boreorepandum occurs in both areas; H. boreorepandum prefers coniferous forests in boreal climate, whereas H. repandum has a wider range of hosts and climate regions (Niskanen et al. 2018). We discovered two ecological groups of H. repando-orientale in Japan: a temperate group (ROF) which occurs in temperate Fagaceae-dominated forests at ≤ 1200 m a.s.l. and a subalpine group (ROC) which occurs in subalpine coniferous-dominated forests at ≥ 2000 m a.s.l. (Sugawara et al. 2022a). A previous study suggested small differences between ROC and ROF groups in terms of morphological characters (basidiospore size) and sequence data (ITS of nc rDNA operon and translation elongation factor 1-alpha (TEF1)). Because we have not found intermediate types in Japan, the two groups may have adapted to their different ecological niches. We tentatively concluded that the ROC and ROF groups are conspecific (= H. repando-orientale s. lat.) because of very small variations in ITS sequences between them (< 0.3%); these variations are within the range of common intraspecific variations in the genus Hydnum (i.e., 0-1% (or 1.5%) for intraspecific variations and > 1% (or > 1.5%) for species delimitation (Niskanen et al. 2018)). However, mating incompatibility among ROC and ROF groups could not be assessed because of the slow growth of their monospore isolates; thus, the taxonomic designations of the two ecological groups of H. repando-orientale require further investigation.
The Hydnum species classification emphasizes on the phylogenetic relationships inferred from ITS sequences (Niskanen et al. 2018;Swenie et al. 2018;Cao et al. 2021b;Sugawara et al. 2022a); however, a phylogeny of the non-coding ITS region alone sometimes leads to unreliable data because it depends on the evolutionary history of a single-DNA sequence, rather than a species. For this reason, the use of multiple molecular markers is strongly recommended in a taxonomic framework (Lücking et al. 2020;Aime et al. 2021;Cao et al. 2021a). Molecular data from multiple loci contribute to species delimitations based on topological congruence and the multispecies coalescent model. The genealogical concordance phylogenetic species recognition (GCPSR) approach proposed by Taylor et al. (2000) is used in Fungi to define phylogenetic species based on congruent clades from multiple genealogies. Coalescent-based species delimitation using multiple loci is a powerful method for estimating evolutionary lineages that involve ancestral polymorphism which lacks reciprocal monophyly among alleles (Fujita et al. 2012). Unfortunately, species delimitation approaches have not been used in the current classification systems in the genus Hydnum. However, these approaches enable the recognition of cryptic (or pseudocryptic) species as part of an integrative taxonomy, together with other species concepts such as morphology, ecophysiology, and hybrid incompatibility (Fujita et al. 2012;Looney et al. 2020;Cao et al. 2021a).
Here, we aimed to (1) reconsider species boundaries among H. repando-orientale and related species in Japan by integrative taxonomic approaches, and (2) present a taxonomic treatment for the ROC group of H. repando-orientale s. lat. under the nomenclature. Multi-locus molecular phylogeny was inferred from the sequences of the ITS and the large subunit (LSU) of nc rDNA, TEF1, RNA polymerase II largest subunit (RPB1), and RNA polymerase II second-largest subunit (RPB2). In addition, species delimitations were performed using the GCPSR approach and the coalescent model based on the generalized mixed Yule coalescent (GMYC) and Bayesian framework (BP&P). A mating test using monospore isolates was performed to evaluate mating compatibility among the ROC and ROF groups, as well as H. boreorepandum, a sister species of H. repando-orientale s. lat. Finally, we performed morphological characterization of basidiomata including a statistical analysis of mean basidiospore size. Based on findings indicating that the ROC and ROF groups are distinct species, we provide a detailed description of the ROC group as a new species.

Basidiomata specimens
We examined 33 basidiomata specimens of H. boreorepandum and H. repando-orientale s. lat. (Table 1): 23 were collected in our studies in 2016-2020 (Sugawara et al. 2019(Sugawara et al. , 2022a; 4 were collected in this study; 2, including holotype of H. repando-orientale (TUMH 60745), were loaned from the Tottori University Mycological Herbarium (TUMH), Fungus/ Mushroom Resource and Research Center, Faculty of Agriculture, Tottori University; 4 were loaned from the TNS herbarium of the National Museum of Nature and Sciences, Tokyo. These specimens were collected at 25 sites in Honshu, mainland Japan, at 10 to 2500 m a.s.l. (Fig. 1). We defined H. repando-orientale s. lat. as the two ecological groups (ROC and ROF) based on the recorded altitude and forest habitat (Table 1).

Morphological analysis
The morphological characterization was conducted in accordance with the method in our previous report (Sugawara et al. 2022a). We observed macroscopic characters of fresh basidiomata and used the color code of the Methuen Handbook of Colour (Kornerup and Wanscher 1978). Microscopic characters of dried specimens were observed in 3% KOH or Melzer's reagent using differential interference contrast microscopy. Thirty basidiospores were measured from each of 21 specimens, and previous measurements were reused for 19 specimens. We calculated length/width ratio of each basidiospore (Q), its mean of a specimen (Q m ), and statistics of basidiospore length and width (5% smallest value (a), mean-SD (b), mean + SD (c), 5% highest value (d)), in which we described each highest or lowest value among all specimens as "(a)b-c(d)." Length or width of other tissues were expressed as range.
Because the ROC group showed slightly larger basidiospores, compared with the ROF group, we statistically analyzed the mean basidiospore sizes (MBS) by unpaired two-sample Wilcox tests. We also included eight MBS values of H. repandum and H. boreorepandum measured by Grebenc et al. (2009), and Niskanen et al. (2018; overall, our analysis included measurements of 5 H. boreorepandum, 7 H. repandum, 11 ROCs, and 16 ROFs. Using the "stats" package in R v. 4.0.5 (R Core Team 2021), the length (L) and width (W) of MBS were compared among ecological/phylogenetic groups (ROC, ROF, H. boreorepandum, and H. repandum) by the function "pairwise.wilcox.test" and the Bonferroni method.

Monospore isolation
Fresh basidiomata materials were used for monospore isolation, in accordance with the method proposed by Sugawara et al. (2019). Basidiospores from each hymenophores were collected on an axenic plastic Petri dish, suspended in sterile distilled water, and inoculated onto modified Norkrans' C medium (MNC; Yamada and Katsuya 1995) solidified with 1.5% gellan-gum (MNC1.5G; Sugawara et al. 2019). Using a platinum loop, inoculated spores were streaked in zigzag to create a spore concentration gradient; they were then incubated at 15°C in the dark (MIR-254-PJ, Panasonic Healthcare, Japan, Tokyo). Infrequent (ca. < 1%) basidiospore germination was observed 1-2 months after incubation. When mycelial colonies reached approximately 2-mm diameter, each was transferred to MNC medium solidified with 1.5% agar (MNC1.5A). To conduct mating tests, we selected monospore isolates that exhibited better growth and lacked clamp connections on hyphal septa. In total, three strains from three ROC basidiomata, six strains from four ROF basidiomata, and three strains from three H. boreorepandum basidiomata were isolated and established from September 2020 to March 2021.

Mating tests
We used only newly obtained strains for mating tests. Monokaryotic strains of the ROC group showed poor growth on MNC1.5A medium (e.g., 0.5-cm diameter per 1-month incubation at 20°C); for this reason, we conducted mating tests three times on different dates and under different culture conditions. The first test involved nine monokaryotic strains that showed better growth on MNC1.5A medium from April 27 to August 4, 2021. Mycelial agar blocks (ca. 3 × 3 mm) pre-cultured for 1-2 months at 15°C were cut and placed near a pair-strain for 3-mm spacing on MNC1.5A medium. Next, each mating pair was incubated at 15°C for 2 months, then at 20°C for 2 months. For the second and third mating tests, each mycelium was pre-cultured in MNC liquid medium for 6 weeks and subsequently transferred to MNC medium plates. From this pre-culture, we could obtain enough mycelial biomass from all monokaryotic strains. Then, second and third tests were conducted using MNC1.5A and one-tenth Site code number indicates collection site in Fig. 1 -: personal number or herbarium number absent. +: specimen or its monospore isolate was examined in each experiment concentration of MNC medium solidified with 2.0% agar (1/10MNC2.0A), respectively. Each test was started at 6 August 2021 or 12 August 2021 and incubated for 2 months at 20°C. The formation of clamp connections was observed under a light microscope at 200× magnification with a long workingdistance objective lens and at 1000× magnification on slide glasses mounted with distilled water. Some strains showed possible crossing between ecological/phylogenetic groups and resulted in intermediate mating incompatibility groups; therefore, we also examined the nuclear phases of hyphae that formed clamp-like cells to determine whether dikaryotization occurs. A portion of mycelium was mounted in a mixture of 4',6-diamidino-2-phenylindole (DAPI; Fujifilm Wako Pure Chemical, Japan, Osaka) 0.04% (v/v) and Calcofluor White (Fujifilm Wako Pure Chemical) 0.02% (v/v), then observed under a fluorescence microscope (Eclipse Ei, Nikon Imaging, Japan, Tokyo) equipped with a mercury lamp (Intenslight C-HGFI, Nikon Imaging). To assess the viability of descendant mycelium in artificial medium, a clamped mycelium was inoculated onto a MNC1.5A plate and incubated at 20°C.

Phylogenetic analyses
Phylogenetic analyses were performed based on the (i) ITS dataset and (ii) concatenated dataset of five loci (ITS, LSU, TEF1, RPB1, and RPB2). The ITS dataset included sequences of Holarctic species of the genus Hydnum downloaded from the INSD/UNITE databases, which comprised > 60 phylogenetic species (Table S2; Olariaga et al. 2012;Yanaga et al. 2015;Niskanen et al. 2018;Swenie et al. 2018;Cao et al. 2021b;Sugawara et al. 2022a). As an outgroup, we selected mycorrhizal Sistotrema spp. Niskanen et al. 2018;Sugawara et al. 2022a); we excluded S. confluens Pers. And S. subconfluens L.W. Zhou because they showed extremely high genetic divergence in the ITS sequences (Niskanen et al. 2018;Sugawara et al. 2022b). The 175 ITS sequences including the outgroup were sampled and automatically aligned using MAFFT online v. 7 (Katoh et al. 2019). The alignment was manually refined, and the best substitution model for RAxML program was estimated using ModelTest-NG v. 0.2.0 (Flouri et al. 2015;Darriba et al. 2020). We identified the best maximum likelihood tree using the rapid bootstrap algorithm in RAxML v. 8.2.10 (Stamatakis 2014) on the raxmlGUI v. 2.0.5 platform (Edler et al. 2021). This analysis was computed under a HKYGAMMA substitution model with 1000 replications of bootstrap analyses (MLBS).
A more detailed phylogenetic analysis was performed using the sequences of five loci (ITS, LSU, TEF1, RPB1, and RPB2) obtained from 24 materials of H. boreorepandum, the ROC group, and the ROF group. Other taxa belonging to the subg. Hydnum were included in the dataset if ≥ 3 sequences of the five loci to be analyzed were available in the INSD. Furthermore, we used newly obtained 45 sequences of eight Hydnum species and four Sistotrema species (Table 2). Thus, 51 total specimens was included in the second analyses.
First, we annotated each locus and constructed independent phylogenies. The ITS and LSU sequences were contiguously connected and subsequently annotated using ITSx v.  (Matheny et al. 2007). TEF1 includes three coding (ca. 500 bp) and two intronic regions (ca. 100 bp); RPB1 comprises mostly a coding region (ca. 850 bp) with a small intronic region (ca. 15 bp); RPB2 has only a coding region (ca. 800 bp). The RAxML phylogenies were independently constructed using each alignment as described above. The following substitution models were selected as recommended by ModelTest-NG: GTRGAMMAI for ITS+ LSU, HKYGAMMAI for TEF1, GTRGAMMAI for RPB1, and GTRGAMMA for RPB2. For phylogenetic analysis of the concatenated dataset using MrBayes v. 3.2.7 (Ronquist et al. 2012), we determined the partitioning schemes and substitution models using PartitionFinder 2 v. 1.1 (Lanfear et al. 2017) under a "greedy" scheme search algorithm (Lanfear et al. 2012) and AICc criteria. By this analysis, we set nine subsets that included independent substitution models as shown in Table 3. Because the TEF1 alignment showed partition scheme trends that differed from the RPB1 and RPB2 alignments, we performed further model estimation and scheme search. Next, PartitionFinder 2 analysis was performed on only the TEF1 alignment under an "all" scheme search algorithm, in which all possible combinations of data blocks were analyzed. This analysis yielded the same partitioning scheme as greedy option; thus, we adopted the scheme in Table 3. The MrBayes analysis was computed with 2,000,000 generations of two iterations of four Markov chain Monte Carlo (MCMC) chains, where trees were sampled every 100 generations.
Chain convergence was confirmed by both visualization by Tracer v. 1.7.2 (Rambaut et al. 2018) and a small value (< 0.01) of average standard deviation of split frequencies (ASDSF). After burn-in of the first 25% generations, the consensus topology was constructed based on the 50% majority-rule of whole topologies. We also constructed concatenated gene trees based on the RAxML and maximum-parsimony methods. The RAxML phylogeny was computed under the GTRGAMMAI substitution model with 1000 replications of bootstraps. The maximum-parsimony trees were inferred by Subtree-Pruning-Regrafting (SPR) methods with 1000 replications of bootstraps using MEGA 7 v. 0.26 (Kumar et al. 2016), and a bootstrap consensus tree was constructed from the resulting 10 unrooted trees.
The alignments of ITS and concatenated dataset are provided in the Supplementary Data. Numerical information concerning the alignments (e.g., numbers of parsimonyinformative sites and distinct alignment patterns) is shown in Table S3.

Species delimitations
Because there was a little conflict among topologies constructed from single alignments, we performed species delimitation using GCPSR (Taylor et al. 2000), GMYC (Pons et al. 2006), and BP&P (Yang 2002;Rannala and Yang 2003) approaches. The GCPSR approach defines a congruence of multi-gene genealogies as a phylogenetic species. In the GCPSR protocol proposed by Dettman et al. (2003Dettman et al. ( , 2006, congruent clades were recognized as genealogical concordance and/ or genealogical non-discordance; they were then defined  Table 1. The map was downloaded from the Geospatial Information Authority of Japan as phylogenetic species via exhaustive subdivision, in which all individuals were required to be placed within a phylogenetic species without creating conflicts with other phylogenetic species. Here, genealogical concordance was defined as the same monophyletic clade recognized in most topologies (≥ 3 of 4 of ITS-LSU, TEF1, RPB1, and RPB2 phylogenies); genealogical non-discordance was defined as the recognition of supported-monophyly (MLBS ≥ 70) in ≥ 1 topology, the clustering of which is never contradicted with the same level of support in other topologies. Finally, specimens were defined as the smallest phylogenetic groups that did not create discordance.
The GMYC approach explores the switching of branching by speciation under the Yule model (interspecific) to neutral coalescent within a species (intraspecific) based on the topology of an ultrametric tree (Pons et al. 2006;Fujisawa and Barraclough 2013). The species delimitation based on GMYC was analyzed using "splits" package in R v. 4.0.5. We generated an ultrametric tree for each locus (ITS-LSU, TEF1, RPB1, and RPB2) under Bayesian inference implemented in BEAST 2 v. 6.6 (Bouckaert et al. 2019). We set a strict clock model for estimating branch lengths and tree priors under the Yule model. The MCMC analysis was performed for 10,000,000 generations and sampled every 1000 steps. The convergence of chain was confirmed by higher values (≥ 200) of effective sample size (ESS) for each parameter on Tracer, and a consensus topology was summarized after a 25% burn-in. For each locus, we removed an outgroup (i.e., Sistotrema spp.) from the topology using "ape" package (Paradis and Schliep 2019) and assigned the topology by GMYC analysis with a single threshold using "splits" package (Thomas et al. 2021).
The BP&P is a Bayesian MCMC program that infers species tree and species delimitation under the multispecies coalescent model using multiple-locus alignments (Yang 2002;Rannala and Yang 2003). We performed an unguided species delimitation of "A11" algorithm in BP&P v. 4.3., which performs joint species delimitation and species tree inference using the reversible-jump MCMC algorithm (Yang and Rannala 2014;Yang 2015). We assigned 26 individuals for four inferenced species as a prior (H. repandum, H. boreorepandum, ROC, and ROF); four alignments (ITS-LSU, TEF1, RPB1, RPB2) were assigned as multiple loci. Because we could not provide the corroborated values of the population size parameters (θs) and the divergence time at the root of the species tree (τ 0 ) for each inverse gamma distribution, we assigned four combinations of rate parameters between higher (0. The evolutionary divergences in each gene were analyzed to estimate overlapping of inter/intraspecific variations within/ between each group. This analysis is useful for determining the optimal DNA barcode (e.g., Harder et al. 2013;Li et al. 2017;Wang et al. 2018). Using MEGA 7, the mean evolutionary divergence "within groups" and "net between groups" were estimated for each locus (ITS, LSU, TEF1, RPB1, and RPB2). The maximum-composite-likelihood method with 1000 replications of bootstrap analysis was implemented to analyze pairwise distances. The patterns among lineages were set as gamma distribution rates, including site heterogeneity.

Phylogenetic analyses
The RAxML phylogram obtained from the ITS dataset showed topology similar to the phylogram in our previous study ( Fig. 2; Sugawara et al. 2022a). Sequences of ROC, ROF, H. repandum, and H. boreorepandum formed a monophyletic clade with high support (MLBS = 89) within the subg. Hydnum (MLBS = 75). Hydnum repandum sequences formed a paraphyletic group with H. boreorepandum and H. repandoorientale s. lat. Hydnum boreorepandum including European and East Asian specimens showed monophyly with strong support (MLBS = 91). Hydnum repando-orientale s. lat. (i.e., the assemblage of ROC and ROF groups) formed a monophyletic clade with strong support (MLBS = 100). Among them, all 18 sequences in the ROF group, including holotype of H. repando-orientale (TUMH 60745), were slightly separated from the sequences in the ROC group as a subclade with low support (MLBS = 61). Three specimens in the ROC group formed a further subclade with the other eight sequences in the ROC group with moderate support (MLBS = 81).

Species delimitations
A summary of each species delimitation is shown in Fig. 4. The GCPSR criterion supported both genealogical concordance and non-discordance of ROC, ROF, and H. boreorepandum, respectively (see "C/nDC" on branches in Fig. 4). A subclade of the ROF group containing three specimens (TNS-F-78326, TUMH 63125, and TUMH 64069) exhibited genealogical nondiscordance. Finally, the exhaustive subdivision process (Dettman et al. 2006) in the GCPSR approach recognized ROC, ROF, and H. boreorepandum as three phylogenetically distinct species and the subclade in ROF as an intraspecific variation.
The GMYC analyses of each topology rejected the null model (likelihood ratio test p-value < 0.001) and supported the assumption that all H. boreorepandum specimens belonged to a single group. The unity of a mixed group of ROC and ROF specimens was supported by the ITS-LSU topology (AICcsupported-value = 1.00) but not by the TEF1, RPB1, and RPB2 topologies (< 0.15). ROC and ROF specimens were clearly separated as two distinct groups based on the RPB1 and RPB2 topologies. In the TEF1 topology, most ROC and  Total 3924 sites * ITS, ITS1, and ITS2 regions. 5.8S and 28S, 5.8S and 28S of nc rDNA, respectively. C1, c2, and c3, first, second, and third positions of the coding region (TEF1, RPB1, and RPB2), respectively. Int, intronic regions of RPB1 and TEF1 ITS1 and ITS2 regions were input as the same scheme. Three loci of TEF1 intronic regions were input as the same scheme ROF specimens were assigned to distinct two groups; one ROF specimen (TUMH 63126) formed a separate group. The BP&P analysis under all combinations of θs and τ 0 priors supported a four-species model composed of H. boreorepandum, H. repandum, ROC, and ROF, with the highest posterior probabilities (0.98-1.00). These analyses supported the assumption that the ROC and ROF groups were two distinct species with high posterior probability support (≥ 0.98) in the all priors set; they did not support a single species model for ROC and ROF. Therefore, all species delimitation approaches supported the assumption that ROC and ROF are distinct species, rather than a single species.
For all loci (ITS, LSU, TEF1, RPB1, and RPB2), the ranges of evolutionary divergences within each ecological group estimated by the maximum-composite-likelihood method were 0.000-0.001 in ROC, 0.000-0.004 in ROF, and 0.000-0.001 in H. boreorepandum (Table 4). The variation of ITS sequences among H. boreorepandum, ROC, and ROF was considerably higher (0.012-0.014) than in H. boreorepandum alone (0.000-0.001); however, the variation between ROC and ROF (0.001) overlapped with the variation within each group (0.000-0.001). Compared with ITS, LSU sequences showed lower values of evolutionary divergences in all groups; thus, three groups could not be distinguished by ITS and LSU markers. However, each of the TEF1, RPB1, and RPB2 genes showed substantially higher variation among the three groups, compared to within the three groups: divergences of single group vs. all three groups were [≤ 0.004 vs. 0.008-0.019] in TEF1, [≤ 0.001 vs. 0.006-0.010] in RPB1, and [0.000 vs. 0.005-0.014] in RPB2 (Table 4).

Mating incompatibility tests
Because clamped hyphae were observed (although infrequently) compared to absent in the original simple-septate hyphae, mating compatibility was evaluated based on the presence of clamp connection at junctions between two confronting colonies. Mating tests showed mating compatibility exists in each ecological/phylogenetic group (Table 5). Although most strains did not show mating compatibility with a strain of a different group, two strains showed mating compatibility beyond their potential incompatibility group (Table 5): one ROC strain (SuR20200920-007 ST03: TUMH 60412) formed clamps with all ROF strains, and one H. boreorepandum strain (SuR20201011-301 ST03: TUMH 64408) formed clamps with several ROF strains (Fig. S5). The clamped hyphae between ROC and ROF strains showed dikaryotization of hyphal cells (Fig. S5e), and we successfully isolated it as a dikaryotic culture strain (SuR20201024-101 ST01 × SuR20200920-007 ST03). Because the clamped hyphae generated by crossing of H. boreorepandum and ROF strains were very rare and sparse, their nuclear phase could not be observed and they could not be isolated in culture; we thus presumed that their dikaryotization failed or the dikaryon lost viability on culture medium. In conclusion, ROF, ROC, and H. boreorepandum groups belong to different mating incompatibility groups but retain partial mating ability; one ROC strain potentially has the ability to form a hybrid with ROF; one H. boreorepandum strain contingently forms clamp-like structures by crossing with ROF but has the lost normal mating ability.
Hydnum subalpinum R. Sugaw. & N. Endo, sp. nov. (Fig. 6). MycoBank no.: 844782 Diagnosis: Hydnum subalpinum is a sister species of H. repando-orientale but differs in that its basidiospores are slightly larger (ave. 8.0-9.1 × 6.9-8.0 μm) and it occurs in subalpine forest habitats associated with Gymnosperm. This species is also related to H. boreorepandum: H. subalpinum and H. boreorepandum show similar morphologies and ecologies but differs in terms of phylogeny and biological isolation, as indicated by in vitro mating incompatibility.
Stipe robust, 20-60 × 9-16 mm, central or eccentric, equal to slightly enlarged at the base, solid, glabrous, whitish cream to cream, not turning color where scratched. Context flesh, whitish cream to cream, when young turning yellowish where scratched. Rhizomorphs abundant, white. Odor mild, strong.
Pileal color of H. boreorepandum and H. subalpinum was slightly whiter, compared with H. repandum and H. repando-orientale (Niskanen et al. 2018;Sugawara et al. 2022a); we presumed that this was affected by environmental condition. Indeed, older basidiomata of H. boreorepandum were tinged with orange hues (Sugawara et al. 2022a). This notation is supported by two TNS specimens of H. subalpinum that were labeled as "H. rufescens" in the subg. Rufescentia Niskanen & Liimat. because of their brownish-orange color. Yanaga et al. (2015) showed that H. repando-orientale (as "H. repandum var. repandum" and "H. repandum var. album") has color variations of the pileus, even at the same collection site.

Discussion
In this study, integrative taxonomic approaches verified that two ecological groups of Hydnum repando-orientale s. lat. constituted two independent species. Phylogeny inferred from five loci demonstrated the distinct divergence between temperate ROF (H. repando-orientale) and subalpine ROC (H. subalpinum) groups. While H. repando-orientale and H. subalpinum have not yet accumulated sufficient variations in ITS, LSU, and TEF1, all species delimitation analyses (GCPSR, GMYC, and BP&P) separated them into two distinct phylogenetic clades; this genetic divergence implies reproductive isolation derived from their different ecological niches. Contrary to our expectations, mating tests among monospore isolates in these two species did not show complete incompatibility between species-one strain showed intermediate incompatibility. However, the other strains exhibited mating incompatibility with other species. Morphological difference between these two species are scarce, but H. subalpinum specimens have significantly larger basidiospores, compared with specimens of H. repando-orientale. We suspect that these two species have lost gene flow because of local adaptation to different vegetation/climate conditions; genetic divergence accumulated independently, resulting in mating incompatibility and phenotypic differences.
One strain of H. subalpinum (SuR20200920-007 ST03: TUMH 64012) clearly showed in vitro dikaryotization with H. repando-orientale strains. This hybridization ability indicates recent speciation into H. repando-orientale and H. subalpinum. However, although the sister species have not established genetic hybrid incompatibility, they presumably have other isolating barriers such as immigrant inviability or ecological hybrid inviability (Nosil et al. 2005) because of their allopatric distribution and niche divergence. In many basidiomycetes, the crossing of different species is frequently observed in mating tests (Le Gac and Giraud 2008). A typical

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.
Competing interests The authors declare no competing interests.