Flammulina filiformis strain Cha01 (Serial number: D01, D represents dikaryon) and its two constituent monokaryotic strains CHA-Y-4 (M01-1, M represents monokaryon) and CHA-Y-25 (M01-2); 00117 (D02) and its constituent monokaryon 00117-Y-1 (M02-1); Fv093 (D03) and its constituent monokaryon Fv093-Y-B (M03-1); 03878 (D04) and its constituent monokaryon 03878-Y-1 (M04-1); 03890 (D05) and its constituent monokaryon 03890-Y-1 (M05-1); Chuan6 (D06) and its constituent monokaryon C6-2-3 (M06-1); Huang1 (D07) and its constituent monokaryon Huang-1-12 (M07-1); Su6 (D08) and its constituent monokaryon Su6-1-2 (M08-1); F2927 (D09); JIN19 (D10); FHJ-12 (D11); JHH (D12) and its constituent monokaryon JIN-1-16 (M12-1); BJP2 (D13) and its constituent monokaryon BJP2-Y-1 (M13-1); F007 (D14) and its constituent monokaryon F007-Y-A (M14-1); F015 (D15) and its constituent monokaryon F015-Y-A (M15-1); WL1073 (D16) and its constituent monokaryon WL1073-Y-3 (M16-1); WS2154 (D17) and its constituent monokaryon WS2154-Y-3 (M17-1); 01922 (D18); JIN4 (D19); FHB01 (D20); FHB07 (D21); FHB021 (D22); F004 and its constituent monokaryon F004-Y-A (M23-1); F006 and its constituent monokaryon F006-Y-A (M24-1) were stored in the Institute of Cash Crops, Hebei Academy of Agriculture and Forestry Sciences. Strain NJ6 and its constituent monokaryon Fv6-3 (As1, as represents the strain was previously released genome assemblies; Liu et al., 2014; Huang et al., 2019); strain 1,123 and its constituent monokaryon W23 (As2) and L11 (As3; Wang et al., 2015); strain C01 and its constituent monokaryon Fv01 (As4) were stored in Fujian Edible Fungi Germplasm Resource Collection Center of China (Yao et al., 2019). The information of strain KACC43778 (constituent monokaryons: KACC42780 and KACC43777) from Konkuk University (Park et al., 2014) and TR19 from Kinki University (Kurata et al., 2016) were referenced. The genome of KACC42780 and TR19 was also used for analysis in this study. The constituent monokaryotic strains of each dikaryon strain in Figure 1 were obtained via protoplast-regenerating of dikaryon and isolating monokaryons (Yoo et al., 2004; Huang et al., 2009; Meng et al., 2021). The detailed information including origin, type and color of all the above strains was described in Figure 1, and all the F. filiformis strains were maintained at 25°C on potato dextrose agar medium (PDA; 200 g/L potato; 20 g/L glucose; 20 g/L agar).

The cultivation of the fruiting bodies of F. filiformis dikaryotic strains was performed per the method described by Tao et al. (2019) with some modifications. The cultivation medium was composed of 52.5% cottonseed shells, 15% sawdust, 25% wheat bran, 5% corn powder, 2% calcium sulfate dehydrate, and 0.5% pulverized lime in 60% water. The cultivation bottles containing 200 g sterile medium were inoculated with PDA blocks with mycelia and placed in an incubator at 23°C for mycelia growth (20 days). After the mycelia filled the bottles, the bottles were moved to an incubator at 10°C and 90% humidity to promote primordia formation and fruiting body development.

Total genomic DNA of F. filiformis was extracted by the cetyl trimethyl ammonium bromide (CTAB) method and sequenced using two platforms including Illumina NovaSeq and BGISEQ500 with a paired-end read length of 150 bp. The raw sequencing reads were processed using Fastp (Chen et al., 2018) with default parameters to clean reads with low quality and adapters.

The clean short reads of all samples were aligned to the reference genome, KACC42780 with 11 complete chromosomes of F. filiformis (BioProject ID: PRJNA191921; Park et al., 2014) using BWA-MEM (Li and Durbin, 2010) with the default parameters. The alignment hits of mapped reads were sorted using SAMtools (Li et al., 2009). The statistics and summary of mapping results were performed using the SAMtools stat command. Variant detection and joint genotyping of multiple samples were performed with FreeBayes using the default parameters (Garrison and Marth, 2012).

To obtain a set of genetic polymorphic markers, variants were filtered to include: (1) only bi-allelic SNPs, (2) minor allele frequency > 0.05, which removed rare variants, (3) at least two samples with homozygous but different genotypes (i.e., both AA and aa homozygous genotypes were required), (4) a missing rate < 0.2, and (5) no variant sites with heterozygous genotypes in monokaryotic samples. This variant filtering procedure was performed using our custom Perl scripts.

The reference genomes of the Fv6-3 (NCBI accession no. PRJNA603211), W23 (NCBI accession no. PRJNA191864), L11 (NCBI accession no. PRJNA191865), Fv01 (NCBI accession no. PRJNA769814), KACC42780 (NCBI accession no. PRJNA191921), and TR19 (NCBI accession no. PRJDB4587) strains were downloaded from the GenBank database of NCBI (Park et al., 2014; Wang et al., 2015; Kurata et al., 2016). Aside from the reference genome KACC42780 used in this study, five other assemblies were chopped into many 5 kb segments with a step of 200 bp when chopping. These chopped segments (simulated long reads) were aligned to the reference genome using Minimap2 (Li, 2018). Variant calls were implemented using BCFtools (Danecek et al., 2021) with the: bcftools mpileup--gvcf, bcftools call, and bcftools merge commands. The allele information for these five assemblies was added to the genotype matrix of the marker set, which contained three types of allele information: missing, 0 (reference allele), or 1 (alternative allele). The haplotype of the reference genome, KACC42780, included all of allele 0 and was also added to the genotype matrix of the marker set.

This study included 13 combinations harboring one sequenced dikaryon and one sequenced monokaryon (Figure 1). Genotype phasing of dikaryon and monokaryon haplotype inference in the 13 combinations were performed. The following two procedures were proposed and implemented using our custom Perl scripts: (1) for all SNP sites, the heterozygous genotypes of the dikaryotic strain were phased using the alleles of its sequenced monokaryotic strain, and (2) based on the phased genotypes of the dikaryon, the un-sequenced monokaryotic haplotype was generated as an individual sample.

The variant call format (VCF) file of the marker set was converted into PLINK files using PLINK 2.0 (Chang et al., 2015). Based on PLINK files, population structure analysis was performed using ADMIXTURE (Alexander et al., 2009) with a parameter K that ranged from 2 to 6. The PCA analysis was performed using PLINK 1.9 (Purcell et al., 2007) with the parameter—pca 20, and the Identity-By-State (IBS) calculation was also performed using PLINK 1.9 with the parameter—distance square ibs.

All the 13 F. filiformis wild strains (D12-D22, F004, F006) collected from China are yellow strains, and 11 cultivated strains collected from China and Japan contained one brown strain (D01), three white strains (D03–D05), and seven yellow strains (D02, D06–D11; Figure 1). The detailed information including origin, type and color of fruiting body of all the above strains was showed in Figure 1. Of them, dikaryotic strains D01-D22 were performed the whole genome resequencing except F004 and F006 (Resequencing of the genomic DNA of F004 and F006 were failed accidentally), and some available monokaryotic strains of potential value (marked with blue triangle in Figure 1) were also selected to perform the whole genome resequencing. This study describes the whole genome sequencing (WGS) of 22 dikaryotic strains and 17 monokaryotic strains of F. filiformis (Figure 1).

A total of 450 Gb sequencing data was generated (Supplementary Table S1), and all strains were sequenced with high depth (above 100-fold, D05 and D13 even reached 1,000-fold). The 39 sets of WGS data involved 22 dikaryotic strains belonging to two categories: wild (n = 11) and cultivated (n = 11). The dikaryotic strain, Cha01 (D01), had three WGS datasets and the 13 dikaryotic strains (D02–D08 and D12–D17) had two WGS datasets (one dikaryon and one monokaryon, respectively). The remaining eight dikaryotic strains (D09–D11 and D18–D22) and two monokaryotic strains (M23-1 and M24-1) had only one set of WGS data each. WGS datasets of the 39 strains were aligned to the reference genome, KACC42780 (NCBI accession no. PRJNA191921). The overall alignment and coverage rates were 69–83% and 85–93%, respectively (Supplementary Table S2). The presence of unmapped reads may indicate that some genomic components are unique to particular strains, suggesting that one reference genome is not enough to represent the whole population of F. filiformis.

Through SNP detection and joint genotyping of all 39 sequenced strains (Figures 1, 2), a set of SNP markers with intrapopulation polymorphism, including 849,987 bi-allelic SNPs, was constructed (the strategy for variant filtering is detailed in the Methods, called this marker set as name of SNP850K). Those variants from a strain-specific genomic region not existing in the reference genome were not included in this marker set. However, this marker set SNP850K was basically constructed from shared genomic regions among sampled strains. These bi-allelic SNP markers basically covered all 11 chromosomes (Figure 2A) with a high distribution density of 24.16 SNP markers per kb. Most of the genomic regions (sliding windows of size 500 kb) had more than 10 SNP markers per kb. Sufficient markers with high density and high resolution ensure the validity of subsequent analyses and the representativeness of the study findings.

Using these high-resolution markers, a quality control (QC) analysis of the samples was conducted and the proportion of missing and heterozygous genotypes from each WGS sample was calculated and visualized (Figure 2B). No sample was an extreme outlier. The missing rate was generally low (<2%). The proportion of heterozygous genotypes showed three different levels: high, medium, and low. All monokaryotic strains exhibited a low level of heterozygous proportion (nearly zero), and two dikaryotic strains [D04 (03738) and D03 (Fv093)] also displayed a low level, indicating a high degree of homozygosity between two monokaryons. There were two dikaryotic strains [D05 (03890) and D14 (F007)] with a medium level of ~20%. The rest of the dikaryotic strains had a high level of homozygosity (>30%). We also estimated the minor allele frequency (MAF) of all markers. The histogram showed that there were many variant sites with low allele frequency (Figure 2C).

The previously released genome assemblies of Fv6-3 (one monokaryon of NJ6 from China), W23 and L11 (two monokaryons of dikaryotic strain 1,123 from China), Fv01 (one monokaryon of C01 which was largely cultivated in Chinese factory and purchased from Japan), KACC42780 (one monokaryon of KACC43778 from Korea), TR19 (dikaryon from Japan) were used for analysis with no extra processing (Figure 1; Supplementary Table S3). Combined with these well-studied strains, an expanded population could help to understand the population structure and breeding history. Each genome assembly was also aligned to the reference genome, KACC42780 (Supplementary Table S3), and the mapped rate was around 85% with a coverage rate between 75 and 84%. After variant detection and genotype identification based on the marker set SNP850K, apart from the reference strain, KACC42780, haplotypes of the other five monokaryotic strains were integrated into the monokaryotic population. The haplotype information (all sites are REF. allele) of the monokaryotic reference genome, KACC42780, was added. A population of monokaryons (n = 35) was constructed that contained 17 sequenced and 13 inferred monokaryons, and five monokaryotic reference genome assemblies.

A population containing 22 dikaryotic strains from the 39 WGS datasets was used for admixture and PCA on the marker set SNP850K. Relatedness analysis between strains was also conducted, using the Identity-by-Stat (IBS) method. The result of cross-validation of multiple K-values (K = 2 to 6) implied that there were most likely three ancestral genetic components in this population (Supplementary Figure S1). Using these results, the 22 dikaryotic strains were divided into three groups (Figure 3). As the respective representatives of the three groups, there were seven non-admixed strains that showed an almost 100% single ancestry component (Figure 3A), while the other strains were admixed. These non-admixed strains were positioned in the outermost area of the PCA plot (Figure 3B), displaying a great genetic distance between different groups. The admixed strains were located in the middle of the PCA plot.

To present breeding history more clearly, the cultivated dikaryotic (CD) strains and the wild dikaryotic (WD) strains were further classified. The CD strains were divided into three subgroups (n = 5, 5, and 1, respectively), and the WD strains were divided into two subgroups (n = 9 and 2, respectively). The WD subgroup 1 (n = 9) contained nine strains (apart from D14 and D20). The WD subgroup 2 (n = 2) contained two non-admixed strains (D14 and D20) that were distantly related to the other strains. The population structure and relationship of the wild strains had varied contributions of genetic components that were less related to each other (Figures 3B,C), suggesting that there was high genetic diversity in the wild population.

The three CD subgroups are described below. CD subgroup 1 (n = 5) included five strains (D01–D05). There were three non-admixed strains (Figure 3A): D03 (Fv093) and D04 (03878), which had high homozygosity and genetic similarity (Figure 3C), and D05 (03890), which was close but not identical to the other strains. The other two strains D01 (Cha01) and D02 (00117) were also related to them (Figure 3C). The CD subgroup 2 (n = 5) included five strains (D06– D10). The PCA result showed that a focused point was overlapped by all five strains (Figure 3B), indicating that they had a highly similar genetic background. High pairwise relatedness scores (IBS values) that were close to 1.0 also confirmed this conclusion (Figure 3C). The two CD subgroups were from the same group and showed strong relatedness (Figure 3C). CD subgroup 3 (n = 1) had only one strain, D11 (FHJ-12), which is a newly bred cultivar crossed between a cultivated monokaryon (CM) from the strain D07 (Huang1) and a wild monokaryon (WM) from the wild strain D12 (JHH) in the WD subgroup 1. Thus, D11 (FHJ-12) and D12 (JHH) were related and clustered (Figure 3C) and were close in the PCA plot (Figure 3B). The dikaryotic strain D11 (FHJ-12) was obtained by crossing of monokaryotic strains M07-1 (HUANG-1-12) and M12-1 (JIN-1-16). There was also relatedness between D11 (FHJ-12) and D07 (Huang1).

An admixture and relatedness for the 35 monokaryotic strains were conducted on the marker set SNP850K. The admixture analysis results clearly showed that there were three well-separated groups (Groups 1–3) that contained 13, 8, and 14 strains, respectively (Supplementary Figure S2; Figure 4A). Each group had a dominant ancestral component and many of the monokaryotic strains were non-admixed (Figure 4A), indicating a 100% contribution from a single ancestral component. There were 12, 4, and 6 (38%) non-admixed strains in the three groups. The cultivated monokaryon (CM) strains primarily came from groups 1 and 2. Moreover, all non-admixed CM strains were highly similar or identical to each other (Figure 4B). The high pairwise relatedness score (IBS values) indicated a narrow genetic background between two monokaryons. Group 3 was almost entirely populated by wild monokaryotic (WM) strains, and the non-admixed WM strains exhibited low pairwise relatedness (Figure 4B). Therefore, it can be concluded that the WM strains were all different from each other, showing a wealth of genetic diversity.

The 35 monokaryotic and 23 dikaryotic strains were combined and PCA analysis was performed on the marker set SNP850K. The results from this integrated PCA analysis (Figure 5) supported the grouping of the strains into three CD and two WD subgroups. The PCA plot (Figure 5A), combined with monokaryons, was consistent with the PCA result of the dikaryotic strains (Figure 3C), e.g., there was the same focused point of CD subgroup 2 [D07 (Huang1) et al]. However, some dikaryons had two monokaryons on either side (Figures 5B–E). In addition, the dikaryon was located at the middle position of its two monokaryons in the PCA plot. Furthermore, for each eigenvector of PC1, PC2, PC3, etc., the value of the dikaryon was equal to the mean of its two monokaryons (Supplementary Table S4). In this study, a total of 14 sets harboring one dikaryon and two monokaryons (Figure 1) had symmetrical monokaryons positions (Figure 5A). This provides evidence for inferring the history of cross-breeding. For example, by using only a dikaryotic population to perform the population structure analysis, it was found that the CD D11 (FHJ-12) strain was closer to the WD D12 (JHH) strain than the CD D07 (Huang1) strain but did not clearly demonstrate the nature of the relationship. In the PCA plot combined with monokaryons, the CD D11 (FHJ-12) strain was in the middle of M12-1 and M07-1, which were the monokaryons of D12 and D07, respectively.