Research Article |
Corresponding author: Patricia Velez ( pvelez@ib.unam.mx ) Corresponding author: Tsuyoshi Hosoya ( hosoya@kahaku.go.jp ) Academic editor: Francesco Dal Grande
© 2020 Jaime Gasca-Pineda, Patricia Velez, Tsuyoshi Hosoya.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Gasca-Pineda J, Velez P, Hosoya T (2020) Phylogeography of post-Pleistocene population expansion in Dasyscyphella longistipitata (Leotiomycetes, Helotiales), an endemic fungal symbiont of Fagus crenata in Japan. MycoKeys 65: 1-24. https://doi.org/10.3897/mycokeys.65.48409
|
During the Last Glacial Maximum (LGM), drastic environmental changes modified the topology of the Japanese Archipelago, impacting species distributions. An example is Fagus crenata, which has a present continuous distribution throughout Japan. However, by the end of the LGM it was restricted to southern refugia. Similarly, Dasyscyphella longistipitata (Leotiomycetes, Helotiales, Lachnaceae) occurs strictly on cupules of F. crenata, sharing currently an identical distribution. As the effects of the LGM remain poorly understood for saprobiotic microfungal species, herein we identified past structuring forces that shaped the current genetic diversity within D. longistipitata in relation to its host using a phylogeographic approach. We inferred present and past potential distributions through species distribution modeling, identifying environmental suitability areas in mid-southern Japan from which subsequent colonizations occurred. Our findings suggest that current high genetic diversity and lack of genetic structure within D. longistipitata are the result of recent multiple re-colonization events after the LGM.
Divergence, gene flow, geographic information systems, haplotype network, host distribution, intraspecific variation, species distribution modeling
Dramatic cyclical glacial advances and retreats characterize the Pleistocene (
Strong evidence suggests that, during the LGM, Japan was extensively covered by coniferous forests of Pinus, Picea, Abies, and Tsuga. Small populations of Fagus, however, were limited to the Pacific coasts of Kyushu and Shikoku (
The Japanese endemic canopy tree Fagus crenata has been extensively studied in this sense, representing an acknowledged model. Independent lines of evidence suggest that by the end of the LMG, the populations of F. crenata were restricted to environmentally stable southern areas (refugia) along the shoreline in southern areas of Japan, and then expanded northwards via two migration routes through the western and eastern shore of Japan (
Although the LGM has long been a focal point for research in a wide range of plant and animal taxa, information on the microbial component remains lacking. Fungal phylogeographic analyses are recently increasing (e.g., to identify glacial refugia, re-colonization routes, and interglacial population expansions;
Dasyscyphella longistipitata Hosoya (Leotiomycetes, Helotiales, Lachnaceae) overcomes the above-mentioned impediments. This culturable ascomycete is highly specific to F. crenata, sharing an identical wide geographical distribution throughout Japan (
In this study, we reconstructed the historical demography of D. longistipitata based on genetic variation data, in addition to information on the geographical distribution and ecological attributes for D. longistipitata as well as for its host F. crenata. The correspondence between both species in their response to major past climate changes was assessed in order to attain a better understanding of the impacts of glacial cycling on the Japanese biota and neglected fungal taxa.
A total of 270 fungal collections were sampled during 2010–2018 in 14 localities covering the entire geographical distribution of the host species F. crenata from Jogakura (Locality 1) in the northeast to Mt. Shibi (Locality 14) in the southwest (Fig.
Geographical distribution of the sampling localities for Dasyscyphella longistipitata associated with cupules of Fagus crenata in Japan. Red dots and numbers correspond to D. longistipitata; whereas blue dots represent F. crenata study sites from
Sampling localities, decimal geographic coordinates, and summary statistics for Dasyscyphella longistipitata across Japan based on the concatenated ITS and beta-tubulin markers. Significance at p < 0.05 in indicated by *. Latitude and Longitude are indicated in decimal following WGS84 datum; number of isolates analyzed (N), segregating sites (S), haplotype number (h), haplotype diversity (Hd), nucleotidic diversity (Pi), Tajima’s D neutrality test parameter (D).
Number | Locality | Voucher specimen | Collection year | Latitude / Longitude | N | S | h | Hd | Pi | D |
---|---|---|---|---|---|---|---|---|---|---|
1 | Jogakura | TNS-F-32166 | 2010 | 40.65806, 140.8378 | 14 | 56 | 14 | 1 | 0.0106 | -1.6175 |
2 | Sado Isl. | TNS-F-38439 | 2013 | 38.06667, 138.2928 | 21 | 60 | 21 | 1 | 0.0099 | -1.3972* |
3 | Tsuchiyu | TNS-F-46835 | 2012 | 37.66302, 140.2851 | 24 | 103 | 23 | 0.996 | 0.0090 | -1.5928* |
4 | Mt. Atema | TNS-F-35002 | 2010 | 37.01183, 138.7792 | 17 | 49 | 17 | 1 | 0.0086 | -1.5025 |
5 | Obora | TNS-F-32140 | 2010 | 36.50295, 138.3289 | 16 | 55 | 16 | 1 | 0.0105 | -1.4934 |
6 | Suganuma | TNS-F-32150 | 2010 | 36.401, 136.8891 | 21 | 70 | 21 | 1 | 0.0101 | -1.7464* |
7 | Mt. Tsukuba | TNS-F-31635 | 2010 | 36.22789, 140.1017 | 20 | 51 | 18 | 0.984 | 0.0080 | -1.5806* |
8 | Mt. Ougiyama | TNS-F-46830 | 2012 | 35.46278, 134.3214 | 20 | 63 | 20 | 1 | 0.0098 | -1.6550* |
9 | Mikunitouge | TNS-F-35001 | 2010 | 35.40472, 138.9164 | 20 | 59 | 20 | 1 | 0.0102 | -1.3855 |
10 | Ashiu | TNS-F-46827 | 2012 | 35.35045, 135.7326 | 12 | 46 | 12 | 1 | 0.0096 | -1.5862* |
11 | Yamakita | TNS-F-32167 | 2010 | 35.22131, 138.9831 | 20 | 50 | 20 | 1 | 0.0092 | -1.1969 |
12 | Mt. Takanawa | TNS-F-32168 | 2010 | 33.94556, 132.85 | 23 | 77 | 23 | 1 | 0.0111 | -1.7589* |
13 | Mt. Shiragadake | TNS-F-54834 | 2018 | 32.15147, 130.9429 | 21 | 74 | 21 | 1 | 0.0097 | -2.4506* |
14 | Mt. Shibi | TNS-F-54822 | 2018 | 31.98015, 130.3676 | 21 | 59 | 21 | 1 | 0.0093 | -1.7216 |
Overall | 270 | 270 | 255 | 0.999 | 0.0096 | -2.4491* |
Culturing, DNA extraction and sequencing methods were conducted following
The following population genetics summary statistics were calculated to estimate the amount of genetic diversity: Number of polymorphic (segregating) sites (S), number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (Pi). The Tajima’s D neutrality test (
In addition, to compare our results with further fungal population-level genetic variation estimates, we analyzed ITS and beta-tubulin sequences from NCBI GenBank (by
To estimate the genetic structure among the 14 sampling localities of D. longistipitata, we calculated the AMOVA-based pairwise PhiST (
We evaluated effective population size changes over time using the Bayesian Gaussian Markov Random Field (GMRF) Skyride Plots (
We implemented the analysis of species distribution modeling (SDM) to evaluate the geographical areas of conserved environmental suitability using the BioClim environmental layers (BIO1–BIO19; http://worldclim.org/bioclim;
The present and past distribution models were generated using the R package biomod2 v 3.1 (Thuiller et al. 2016) by means of ensemble models (
Areas of conserved environmental suitability across time were delimited as follows: the ensemble layers were transformed to presence/absence considering a threshold of 90% of the distribution of suitability values, and the resulting layers were merged. The area of the resulting layer was classified as highly conserved (overlapping of four layers), medium (at least three overlapping layers), and low (two overlapping layers).
Finally, to estimate the ancestral locations and the spatial dynamics of the two species, we performed a continuous-space phylogeographic analysis using BEAST software. Runs were performed using the molecular clock rates, and substitution models above mentioned. The statistics for spatial traits (i.e., localities) were generated using the Cauchy RRW model. Also, to avoid noise due to duplicate location traits, a jitter window size of 0.01 was used. Posterior to convergence checking, the consensus discrete trait-annotated trees were obtained by means of treeAnnotator using the “Common Ancestor Heights” option and discarding the first 10% trees. Then, this information was used as an input in SPREAD3 v1.0.7 (
We analyzed a total of 270 isolates (12–24 per site). Overall, 85 ITS-based haplotypes were identified. The majority was represented by H12 (40.7%), followed by H28 (14.4 %), and H10 (10.7 %), and the 27.8 % were singletons. Both H12 and H28 were found in all the 14 sites, while H10 was found from 12 sites.
In contrast, the beta-tubulin sequences were remarkably diverse, with 224 recognized haplotype patterns. The most frequently observed haplotypes were B4 (3.7 %), B27 (2.2 %), and B75 (1.9 %), and the majority (74.8 %) was occupied by singletons. The ITS and beta-tubulin concatenated sequences ranged from 980–1022 bp, resulting in 255 haplotypes with 270 segregating sites (Table
The paired PhiST yielded low values, ranging from 0 to 0.0358 (Suppl. material
Multivariate analyses of the genetic diversity inferred from ITS and beta-tubulin concatenated sequences of Dasyscyphella longistipitata A Principal Component Analysis (PCA) of genetic diversity at the individual level B principal Correspondence Analysis (PCoA) of genetic diversity using the localities as grouping factor. Colors represent the locality of origin arranged as a latitudinal gradient where red represents the further north site.
The haplotype network revealed intricate relationships among D. longistipitata haplotypes (Fig.
Haplotype network based on the concatenated sequences of ITS and beta-tubulin of Dasyscyphella longistipitata. The size of the circles represents the haplotype frequency; white dots represent mutational steps between haplotypes (note that the branches lengths do not correspond to genetic distances). Colors represent the locality of origin arranged as a latitudinal gradient where red represents the further north site.
The Skyride plots (Fig.
Bayesian Skyride Plot for A Dasyscyphella longistipitata using the concatenated ITS and beta-tubulin, and B Fagus crenata, using the reported sequences in
The SDM of D. longistipitata resembles the current distribution of F. crenata (
Our results on the ancestral localities and dispersion trends, suggested that D. longistipitata originated from Mt. Ougiyama (site 8), Obora (site 5), Ashiu (site 10), and Suganuma regions (site 6; Fig.
Saprophytic fungi (those obtaining nutrients from dead organic matter) are critical to the dynamics and resilience of ecosystems (
One of the most remarkable results of this study was the unexpected high genetic diversity in D. longistipitata (pi = 0.0094; Table
The estimated nucleotide diversity in D. longistipitata was lower compared to C. purpurea (pi = 0.0101), which is reasonable, as C. purpurea comprises three divergent lineages (
Several factors influence the genetic structure in fungal populations, due in part to the diverse ecological roles of these osmotrophs. For instance,
Even though the sampled localities are distributed across a broad geographic range covering ~1700 km, our genetic structure estimates were practically negligible (Fig.
According to the negative values of the neutrality tests and the star-shaped connections in the haplotype network, the Bayesian Skyride plots revealed that both D. longistipitata and F. crenata underwent important recent demographic growth events. However, the effective-size of D. longistipitata was several orders of magnitude larger than F. crenata. For instance, assuming a fixed mutation rate of 1×10^-10 for both species, we would expect a current effective-size of ~36×10^9 for D. longistipitata, and 35,000 for F. crenata. Moreover, considering a generation time of 10 years for F. crenata along with a mutation rate range of 1×10^-11 to 1×10^-12 per site/year (which is feasible for the matk gene given the reported mutation rate and the slow evolution rates reported for Fagaceae family;
Aside from the lack of mutation rates for the genetic markers, and precise reports on the number of reproductive events per year (i.e., the generation time), our rough calculations on the estimates for population growth agree with the hypothesis of post-glacial demographic expansion. Furthermore, the demographic increase trend in D. longistipitata may be related to a higher mutation rate and shorter generation times. So, growth trends for D. longistipitata in the Skyride plots could correspond to a fraction of the time scale in F. crenata. In this sense, the mast seeding (large, synchronic seed production) in F. crenata (
The SDM analyses revealed several restricted areas of climatic stability distributed in mid-southern Japan (Fig.
Modifications in species distribution ranges result from the interaction between ecological and evolutionary processes (
Dasyscyphella longistipitata has been recognized as a fine woody debris decomposer, in particular of cupule litter, an essential component of litterfall in beech forests (
This project was funded by Grant-in-Aid for scientific research (C) 17570088. We thank Dr. A. Yabe for his valuable discussion about past distribution of Fagus crenata in Japan. We also appreciate sampling assistance by Kazuaki Tanaka, Harue Abe, Makoto Hashiya, Etsuko Kurokawa, Kentaro Hosaka, Takashi Shirouzu, Hideji Ushijima, Konomi Yanaga, Yoko Kimura, Daisuke Sakuma, Tomio Okino, Miyuki Kobayashi. This paper was written during a research stay of PV at the National Museum of Nature and Science, Tsukuba, with support of FY2018 JSPS Invitational Fellowship for Research in Japan (ID No. S18062). We acknowledge anonymous reviewers for constructive comments on earlier versions of the manuscript.
Figure S1. Species distribution modeling (SDM) for Dasyscyphella longistipitata.
Data type: occurrence
Explanation note: A Projection for the current major environmental conditions. Color bar represents environmental suitability values, where warmer colors indicate higher suitability, B projections for the four time periods analyzed in this study. LGM = Last Glacial Maximum, LIG = Last Interglacial period. The areas represent the upper 90% of the distribution of the suitability values for each period of time.
Table S1
Data type: statistical data
Explanation note: Summary statistics for ITS and beta-tubulin markers, showing the sample size (N), number of polymorphic sites (S), number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (Pi).
Table S2
Data type: statistical data
Explanation note: Paired PhiST values between D. longistipitata localities (1–14). Bold numbers denote significance at p < 0.05. For nomenclature on sampling localities please refer to Table
Table S3. List of isolates and their haplotype properties
Data type: (measurement/occurrence/multimedia/etc.)
Explanation note: All the isolates obtained were numbered with a suffix “DL-”, and their ITS-5.8S and beta-tubulin regions were sequenced. Only isolates that were successful sequenced in both regions were used. In total, 270 isolates were used.
Table S4. List of accession numbers for ITS-5.8S and beta-tubulin haplotypes
Data type: (measurement/occurrence/multimedia/etc.)
Explanation note: When a haplotype was obtained from multiple isolate, one isolate was randomly chosen and its sequence was registered. SeqID is for the authors’ private use, and indicated for convenience.