A global phylogenomic analysis of the shiitake genus Lentinula

Abstract

Lentinula is a broadly distributed group of fungi that contains the cultivated shiitake mushroom, L. edodes. We sequenced 24 genomes representing eight described species and several unnamed lineages of Lentinula from 15 countries on four continents. Lentinula comprises four major clades that arose in the Oligocene, three in the Americas and one in Asia-Australasia. To expand sampling of shiitake mushrooms, we assembled 60 genomes of L. edodes from China that were previously published as raw Illumina reads and added them to our dataset. Lentinula edodes sensu lato (s. lat.) contains three lineages that may warrant recognition as species, one including a single isolate from Nepal that is the sister group to the rest of L. edodes s. lat., a second with 20 cultivars and 12 wild isolates from China, Japan, Korea, and the Russian Far East, and a third with 28 wild isolates from China, Thailand, and Vietnam. Two additional lineages in China have arisen by hybridization among the second and third groups. Genes encoding cysteine sulfoxide lyase (lecsl) and γ-glutamyl transpeptidase (leggt), which are implicated in biosynthesis of the organosulfur flavor compound lenthionine, have diversified in Lentinula. Paralogs of both genes that are unique to Lentinula (lecsl 3 and leggt 5b) are coordinately up-regulated in fruiting bodies of L. edodes. The pangenome of L. edodes s. lat. contains 20,308 groups of orthologous genes, but only 6,438 orthogroups (32%) are shared among all strains, whereas 3,444 orthogroups (17%) are found only in wild populations, which should be targeted for conservation.

Keywords: domestication; evolution; fungi; mushrooms; population genomics.

Fig. 1.

Genome characteristics of Lentinula. (A) Nonrepeat regions and number of genes. The vertical lines indicate the average sizes of the groups. (B) Repeat contents of the eight PacBio assemblies. Repeat regions were estimated using the JGI annotation pipeline. LTR: long terminal repeats.

Fig. 2.

Phylogenetic relationships, clade ages, genome sizes and contents, and classification of Lentinula. Chronogram was produced using BEAST. Node bars indicate 95% h.p.d. of age estimates. All branches receive 100% bootstrap support unless otherwise indicated. KO = Korea, JP = Japan, CN = China, RU = Russia, TH = Thailand, VN = Vietnam, NP = Nepal, AU = Australia, NZ = New Zealand, PNG = Papua New Guinea, CR = Costa Rica, US = United States, BR = Brazil, and MX = Mexico. Numbers following isolate and country codes indicate ITS groups. Core genes are those conserved in all genomes, while ITS-group specific genes include the genes conserved in all members of a group but lacking in others.

Fig. 3.

Species or minimal cluster trees inferred with STACEY with matrices showing posterior probability densities for species clusters obtained at collapse height (ε) settings of 0.0001 (Left), 0.0002 (Middle), and 0.0003 (Right).

Fig. 4.

Population structure in L. edodes s. lat. based on 3 × 10⁶ SNPs obtained from assemblies. (A) SplitsTree graph of SNP data, with Admixture pie graphs at K = 2. (B) Admixture results at K = 2, 3, and 4. (A and B) Labels for isolates with group 5 ITS type are shaded in purple (others are group 1). Labels for cultivars are shaded in gray (others are wild isolates). (C) PCA of SNPs.

Fig. 5.

The best-supported phylogenetic SNaQ network inferred from pseudolikelihood inference with incomplete lineage sorting based on 196 gene trees from an 86-genome dataset and enforcing M1 and M2 clades as hybrids, implemented in PhyloNet. The number of genomes per group is given in parentheses. Both M1 and M2 were resolved in a clade with L. edodes W with different amounts of introgression from L. edodes C + W. M1 as sister to the M2, and the W clade shows more introgression with the L. edodes C + W clade. This pattern may be indicative of cultivated strains being reintroduced to hybridize with natural populations of the W clade. Inset: Comparison of model likelihoods for different maximum number of putative hybridizations (hmax) shows the best model at hmax = 3.

Fig. 6.

Rarefaction curve showing pangenome (blue line) and core genome (red line) in L. edodes s. lat. X axis indicates number of isolates sampled; Y axis indicates number of orthogroups in each category (an orthogroup may contain one or more gene copies). Core are orthogroups that are present in all isolates up to that point, whereas the pangenome contains orthogroups present in at least one isolate prior to that point. Each defined range of isolates (cultivars: Left; wild: Right) is produced using 1,000 iterations with random genome addition sequences. Violin plots indicate the distribution of values for the different permutations.