Genetic Diversity and Genome-Wide Association Study of Pleurotus pulmonarius Germplasm

Abstract

Pleurotus pulmonarius is prized by consumers for its distinct flavor, strong aroma, and dense, crispy texture. Although China has extensive germplasm resources for P. pulmonarius, only a limited number of cultivars are commercially available. A comprehensive evaluation and detailed analysis of P. pulmonarius germplasm, alongside the exploration of superior germplasm resources, are essential for developing new varieties. In this study, we resequenced the genomes of 47 P. pulmonarius strains collected nationwide, identifying a total of 4,430,948 single nucleotide polymorphism (SNP) loci. After filtering based on minor allele frequency and data integrity, 181,731 high-quality SNP markers were retained. Phylogenetic analysis grouped the strains into six clusters, with strains from similar geographical regions clustering together. Most CBS strains formed a single cluster; cultivated varieties exhibited higher genetic similarity, whereas wild strains displayed greater diversity. Principal component analysis (PCA) and population structure analyses, using the same SNP markers, corroborated the phylogenetic findings. DNA fingerprinting, derived from 369 core SNPs, further underscored the genetic diversity among strains. Significant morphological variation was observed, with strains in groups ZP, CBS, and WHLJ exhibiting notably higher yields and cap widths compared to other groups. Correlation analysis revealed associations among various phenotypes, while genome-wide association study (GWAS) identified multiple SNP markers within candidate genes linked to agronomic traits, most of which were controlled by multiple genes. This research offers a molecular-level characterization and evaluation of P. pulmonarius germplasm resources, providing a scientific basis for enriching available germplasm and advancing breeding materials.

1. Introduction

Pleurotus pulmonarius, a member of the class Agaricomycetes and subclass Agaricomycetidae, is one of the main species within the oyster mushroom group. According to the China Edible Fungi Association, P. pulmonarius production in China reached 636,700 tons in 2022. China has abundant P. pulmonarius germplasm resources, with wild populations found in regions such as Jiangsu, Zhejiang, Guangdong, Fujian, Henan, and Guangxi. The successful cultivation of P. pulmonarius was first achieved in Luoyuan County, Fuzhou City, in 1998, leading to its rapid and widespread adoption thereafter [1]. Known for its high nutritional value, P. pulmonarius is rich in protein, amino acids, sugars, unsaturated fatty acids, vitamins, folic acid, and essential trace elements. Its distinct aroma, firm yet tender texture, and mild flavor make it highly desirable to consumers.

Studying and assessing the genetic diversity of P. pulmonarius germplasm is essential for breeding applications, providing a foundation for identifying superior germplasm and selecting parental lines with optimal genetic backgrounds. Mating type analysis of 30 P. pulmonarius cultivars revealed only four A factors and three B factors, indicating high genetic similarity among cultivars. The A factor controls nuclear binding and binucleation, while the B factor regulates hyphal fusion and nuclear migration [2]. Clonal propagation through tissue culture is commonly practiced, occasionally leading to the misidentification of distinct strains under the same name. Therefore, accurate variety identification and genetic relationship analysis between cultivated and wild strains are crucial for ensuring high-quality production and advancing scientific research. With the rapid development of molecular biology, molecular markers have become invaluable for analyzing interspecies genetic differences. In addition to traditional morphological identification and antagonistic reactions, molecular markers such as RFLP, RAPD, ISSR, AFLP, and SNPs are widely used to analyze the genetic diversity and phylogenetic relationships of P. pulmonarius germplasm [3]. However, these methods often face limitations in marker quantity, reproducibility, and stability. In contrast, single-nucleotide polymorphisms (SNPs) are abundant, stable, and suitable for automated, high-throughput applications, making them the most widely used genetic markers today. As a primary genotyping marker, SNPs are extensively applied in identifying edible fungi and analyzing genetic diversity across strains, including species like Lentinus edodes [4,5,6], Flammulina filiformis [7,8], and Agaricus bisporus [9]. SNPs are also critical in analyzing P. pulmonarius germplasm diversity, identifying hybrids, and constructing genetic linkage maps. Genome resequencing enables the acquisition of SNP markers evenly distributed across the genome [10]. For example, SNP markers have been used to identify laccase genes in Lentinus edodes [6]. Okuda et al. constructed a genetic linkage map for P. pulmonarius using 150 single-spore hybrid populations and 300 AFLP markers along with two mating-type factors and sporeless traits. This map comprises 12 linkage groups, covering a total length of 971 cM [11]. Vidal-Diez et al. recently assembled and annotated high-quality genomes of P. pulmonarius strains ss2 and ss5. Each genome consists of 23 scaffolds, with sizes ranging from 5.55 Mb to 11 kb (ss2) and 5.06 Mb to 21 kb (ss5), and N50 values of 3.2 Mb and 3.4 Mb, respectively, providing valuable reference genomic data for developing genome-wide SNP markers [12].

Genome-wide association study (GWAS) is commonly employed to identify genes associated with agronomic traits across the entire genome [13,14,15,16]. However, applications of GWASs in edible fungi remain limited. Li et al. used 297 genome-wide markers to genotype 89 cultivated varieties of L. edodes. Among these, 43 markers were significantly associated with four phenotypic traits, leading to the identification of 97 nearby genes and establishing a foundation for molecular breeding in L. edodes [17]. Corner et al. classified the genus Pleurotus into three categories based on hyphal types [18]. Morphological characteristics in edible fungi are highly susceptible to environmental influences [19]. Using a GWAS, Yu identified that DNA damage repair and protein translation processes may play roles in L. edodes responses to cadmium (Cd) stress, providing valuable genomic and population data for functional genomics and artificial breeding research [20]. Zhang applied bulked-segregant analysis (BSA) and an extreme-phenotype GWAS (XP-GWAS) to map cap color traits in Pleurotus cornucopiae and identified the tyrosinase-encoding gene PcTYR through comparative transcriptome analysis [21]. Traditional QTL analysis primarily relies on genetic diversity among parent populations, and its detection efficiency varies between populations. As it often involves multiple genes, accurately pinpointing candidate genes is challenging. In contrast, GWAS enhances allele diversity and improves mapping accuracy. Advances in sequencing technology and mutation analysis tools have significantly increased marker density, while the incorporation of mixed linear models with random and fixed effects has substantially reduced false-positive results.

In this study, the genomes of 47 P. pulmonarius strains collected from various regions across the country were resequenced to detect genome-wide genetic polymorphisms. SNP markers were screened to construct an SNP fingerprint for the strains. Using these genome-wide SNP markers, a phylogenetic tree was developed to analyze the genetic relationships among the strains. Additional population genetic analyses, including population structure, principal component analysis, and linkage disequilibrium, were performed. Fruiting body-related traits of the strains were measured, and the phenotypic data were statistically analyzed. Finally, GWAS between SNP markers and phenotypic traits was conducted to identify SNP markers with significant influence on the traits. This study provides a theoretical foundation for future screening of superior germplasms, gene mining for trait variation, and targeted breeding efforts.

2. Materials and Methods

2.1. Experimental Materials

The tested strains were all stored in the China Center for Mushroom Spawn Standards and Control (CCMSSC) (

Table 1. Tested strains and their origins.

Sample Name	Group	Sample Name for Analysis	Name Interpretation
CCMSSC00313	ZP 1	ZP1	Cultivation
CCMSSC00493	ZP	ZP2	Cultivation
CCMSSC00494	ZP	ZP3	Cultivation
CCMSSC00498	ZP	ZP4	Cultivation
CCMSSC00499	ZP	ZP5	Cultivation
CCMSSC03815	ZP	ZP6	Cultivation
CCMSSC03836	ZP	ZP7	Cultivation
CCMSSC03886	ZP	ZP8	Cultivation
CCMSSC03897	ZP	ZP9	Cultivation
CCMSSC04423	ZP	ZP10	Cultivation
CCMSSC04537	ZP	ZP11	Cultivation
R08026	ZP	ZP12	Cultivation
R08027	ZP	ZP13	Cultivation
CCMSSC01106	WSC	WSC1	Sichuan wild
CCMSSC01108	WSC	WSC2	Sichuan wild
CCMSSC01109	WSC	WSC3	Sichuan wild
CCMSSC01110	WSC	WSC4	Sichuan wild
CCMSSC01111	WSC	WSC5	Sichuan wild
CCMSSC01114	WSC	WSC6	Sichuan wild
CCMSSC01116	WSC	WSC7	Sichuan wild
CCMSSC01118	WSC	WSC8	Sichuan wild
CCMSSC01119	WSC	WSC9	Sichuan wild
CCMSSC01120	WSC	WSC10	Sichuan wild
CCMSSC01123	WSC	WSC11	Sichuan wild
CCMSSC04199	WHLJ	WHLJ	Heilongjiang wild
CCMSSC04458	CBS	CBS1	CBS
CCMSSC04459	CBS	CBS2	CBS
CCMSSC04460	CBS	CBS3	CBS
CCMSSC04461	CBS	CBS4	CBS
CCMSSC04462	CBS	CBS5	CBS
CCMSSC01301	WYN	WYN1	Yunnan wild
CCMSSC04584	WYN	WYN2	Yunnan wild
CCMSSC04585	WYN	WYN3	Yunnan wild
CCMSSC04586	WYN	WYN4	Yunnan wild
CCMSSC04587	WYN	WYN5	Yunnan wild
CCMSSC04588	WYN	WYN6	Yunnan wild
CCMSSC04592	WYN	WYN7	Yunnan wild
CCMSSC04593	WYN	WYN8	Yunnan wild
CCMSSC04591	WSX	WSX	Shaanxi wild
CCMSSC04594	WLN	WLN1	Liaoning wild
CCMSSC04595	WLN	WLN2	Liaoning wild
CCMSSC04597	WJL	WJL1	Jilin wild
CCMSSC04598	WJL	WJL2	Jilin wild
CCMSSC04599	WHB	WHB1	Hubei wild
CCMSSC04600	WHB	WHB2	Hubei wild
CCMSSC04601	WHB	WHB3	Hubei wild
CCMSSC04602	WHB	WHB4	Hubei wild

¹ The cultivated strains are labeled with the format ‘ZP’ followed by a numerical code, while the wild strains are labeled with ‘W’ followed by the abbreviation of the collection region and a numerical code.

).

2.2. DNA Extraction and Sequencing

Genomic DNA from the tested strains was extracted from mycelia grown on PDA (Difco™ Potato Dextrose Agar) plates using the cetyltrimethylammonium bromide (CTAB) method [22]. DNA concentration and purity were assessed with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Sequencing libraries of 350 bp fragments were constructed and sequenced on the Illumina/BGI platform (Beijing BioMarker Technology Co., Ltd., Beijing, China), generating 150 bp paired-end reads.

2.3. Identification and Screening of SNP

Single nucleotide polymorphisms (SNPs) were primarily detected using the GATK software (v4.0) toolkit [23]. Redundant reads were filtered with samtools (v1.9) based on localization results of clean reads aligned to the reference genome to ensure accuracy [24]. GATK’s HaplotypeCaller algorithm was employed to detect SNP and InDel mutations, generating gVCF files for each sample, followed by joint genotyping across groups. Rigorous filtering was applied to ensure reliability. The main filtering parameters included using vcfutils.pl in bcftools to exclude SNPs within 5 bp of an InDel and adjacent InDels within 10 bp; limiting variations to a maximum of two in any 5 bp window; and excluding variants with a Phred quality score (QUAL) below 30, a quality-to-depth ratio (QD) below 2.0, a mapping quality (MQ) below 40, and a Fisher Strand value (FS) above 60. Default GATK parameters were used for additional filtering. These steps ensured the accurate and reliable detection of SNPs and InDels.

2.4. Population Genetic Analysis and Construction of Fingerprints

In the analysis, all selected high-quality SNP markers were used for population genetic analysis. MEGA X software [25] was applied to construct phylogenetic trees for each sample using the neighbor-joining method with the p-distance calculation model and 1000 bootstrap replications. Population genetic indicators were calculated using VCFtools (v0.1.15) [26]. Admixture software assessed population structure [27], while EIGENSOFT [28] performed principal component analysis to cluster the samples. Genetic diversity analysis was conducted using PLINK software [29]. A DNA fingerprint was created to identify markers suitable for genetic typing based on specific criteria: markers were evenly distributed across the genome, exhibited 100% completeness (no missing data), had a minor allele frequency (MAF) above 20%, a polymorphic information content (PIC) above 0.35, conformed to Hardy–Weinberg equilibrium (p-value > 0.01), and showed no nearby mutations within 100 bp of the marker [30].

2.5. Cultivation Test and Trait Screening

The 47 tested strains were cultivated in polycarbonate plastic bottles (diameter: 70 mm, height: 90 mm, volume: 280 mL) filled with 180 g (wet weight) of substrate. The substrate consisted of cottonseed shell (94%), wheat bran (5%), and gypsum (1%) with approximately 65% water content. The bottles were sterilized for 2 h at 0.15 MPa and 126 °C. Six replicate bottles were inoculated for each strain. Once fully colonized, the bottles were transferred to a fruiting room maintained at 16 °C to 18 °C, with 80% to 90% relative humidity and a 12 h photoperiod (300 to 350 lx). CO₂ concentration was kept below 1000 ppm through ventilation.

Phenotypic data were collected for various agronomic traits, including the period from inoculation to primordia generation (PPG) and the period from inoculation to harvesting (PIH) for each bottle. Measurements of mature fruiting bodies included cap length (CL), cap width (CW), stipe length (SL), and stipe diameter (SD). Yield of fruiting bodies per bottle (YFB) was also recorded. All strains were cultivated under the same climate-controlled conditions. Broad-sense heritability (H²) was calculated across tester lines as $H^2={\sigma }_G^2/[{\sigma }_G^2+\left({\displaystyle \frac{{\sigma }_e^2}{n}}\right)]$, where ${\sigma }_G^2$ represents the genotypic variance, ${\sigma }_e^2$ represents the residual variance, and n is the number of testers (indicate heterokaryon set for which this was the case) (n = 42). The phenotypic data were recorded individually and analyzed statistically. One-way analysis of variance and correlation analysis were performed using SPSS software (v27.0) [31]. The best linear unbiased prediction (BLUP) [32] method was used to estimate trait values for association analysis. Descriptive statistics were calculated for each phenotype, including the number of samples after removing missing values (SamNum), mean (Mean), standard deviation (SD), median (Median), minimum (Min), maximum (Max), range (Range), and coefficient of variation (CV). This comprehensive approach was designed to assess the growth and phenotypic characteristics of P. pulmonarius strains under controlled cultivation conditions, enabling a detailed analysis of mushroom production traits.

2.6. Genome-Wide Association Study

SNP data were filtered based on a minor allele frequency (MAF ≥ 0.05) and data integrity (INT ≥ 0.8). Phenotypic data were used to construct frequency distribution histograms. GEMMA software [33] employing the linear mixed model (LMM) was used for GWAS, represented by the model:

$$y=W\alpha +x\beta +\mu +e$$

The analysis generated p-values, which were used to create Manhattan and QQ plots. A threshold was set by dividing the number of effective markers by quality control criteria, with regions meeting −log₁₀(p) ≥ 5 designated as candidate intervals. Functional annotation was applied to genes identified within these candidate regions to clarify their potential biological roles. All experimental steps were completed in 2023.

3. Results

3.1. Quality Statistics of Sequencing Data and Discovery of SNPs

Whole-genome resequencing of P. pulmonarius strains yielded clean reads ranging from 6,119,174 to 11,448,795, with an average of 7,998,114 reads per strain. The Q30 quality score ranged from 88.4% to 95.81%, with an average of 92.66%. GC content varied from 43.58% to 50.35%, averaging 47.33%. In total, 112.15 Gbp of clean data was generated across the 47 strains, with an alignment rate to the reference genome of 75.02%, an average coverage depth of 41×, and genome coverage of 89.84% (at least one base covered). These metrics satisfied the requirements for sequence analysis.

Using the GATK toolkit, 4,430,948 SNPs were detected between the samples and the reference genome (NCBI: ASM1298053v1). The largest number of SNPs was observed between cultivated varieties ZP1 and ZP13 (234,642 SNPs), while the smallest number was observed between ZP2 and ZP12 (2151 SNPs). Analysis of SNP distribution across the genome showed that 65.51% of SNPs were located within gene-coding regions, 13.97% were within 5 kb upstream of genes, and 5.95% were in intronic regions. Synonymous coding mutations accounted for 87.13% of the total, whereas nonsynonymous coding mutations constituted 12.69%. Other types of mutations, such as SNPs within 5 kb downstream of genes, accounted for 10.59% (Figure 1).

3.2. Population Genetic Analysis Based on Whole-Genome Wide SNP Markers

3.2.1. Phylogenetic Relationship

After filtering, a total of 181,731 high-quality SNP markers were obtained. Phylogenetic analysis revealed distinct clusters among the tested P. pulmonarius strains. Cultivated varieties (ZP8, ZP11, ZP10, ZP7, ZP3, ZP6, ZP4, ZP12, ZP2, ZP9, and CBS3) grouped closely within Cluster I, indicating high genetic similarity among cultivars (Figure 2). In contrast, two other cultivated varieties, ZP5 and ZP13, clustered with a wild strain from Yunnan in a separate branch within Cluster I, suggesting that these cultivars may have been domesticated from wild strains originating in Yunnan. Most strains collected from CBS clustered in Cluster II, except for wild strain WSC5 from Sichuan and cultivated strain ZP1, which also fell into this cluster but on different branches.

The wild strains exhibited relatively high genetic diversity, with strains from similar geographical regions clustering together. For instance, wild strains from Liaoning, Jilin, and Heilongjiang grouped in Cluster III, while those from Sichuan were predominantly found in Clusters IV and V, and those from Hubei in Cluster VI.

PCA analysis illustrated the clustering pattern of the samples, confirming the results of the phylogenetic analysis (Figure 3, Table S1). The two-dimensional plot of PC1 and PC3, as well as the three-dimensional plot of PC1, PC2, and PC3, depicted the spatial distribution of each strain based on genetic variation. PC1 accounted for 36.58% of the variation, PC2 explained 12%, and PC3 accounted for 6.52%. The analysis revealed four distinct clusters among the strains. With the exception of ZP1 and ZP5, most strains and CBS3 formed a cohesive group characterized by PC1 values ranging from 0.2125 to 0.2122, closely matching the branches observed in the phylogenetic tree. Strains such as WYN8, WYN7, WYN3, WYN4, and WYN2 exhibited PC1 values between 0.1507 and 0.1368, while ZP5, ZP13, WYN6, ZP1, WSC5, WYN5, WYN1, and WSX ranged from 0.0505 to 0.0491, aligning closely with branches in Clusters I and II in the phylogenetic tree.

Wild strains, represented by Clusters III, IV, V, and VI from the phylogenetic tree, showed a broader distribution range of PC1, spanning from −0.1236 to −0.1385. The PCA results closely aligned with the phylogenetic tree, reflecting high genetic similarity among cultivated varieties and significant variability and divergence among wild strains.

3.2.2. Genetic Diversity Analysis

Polymorphic markers among strains within each group exhibited average polymorphic information content (PIC) ranging from 0.226 to 0.327 (

Table 2. Analysis table of population genetic diversity of P. pulmonarius.

Group	Average MAF	Nei Diversity Index	Observed Heterozygous Number	Polymorphism Information Content
CBS	0.229	0.200–0.667 (0.361) 1	0.200–1.000 (0.243)	0.164–0.375 (0.265)
WHB	0.255	0.250–0.667 (0.399)	0.250–1.000 (0.409)	0.195–0.375 (0.281)
WJL	0.326	0.500–0.667 (0.550)	0.500–1.000 (0.606)	0.305–0.375 (0.326)
WLN	0.328	0.500–0.667 (0.552)	0.500–1.000 (0.598)	0.305–0.375 (0.327)
WSC	0.202	0.091–0.556 (0.303)	0.091–1.000 (0.293)	0.083–0.375 (0.238)
WYN	0.231	0.125–0.667 (0.343)	0.125–1.000 (0.296)	0.110–0.375 (0.260)
ZP	0.207	0.077–0.533 (0.287)	0.077–1.000 (0.435)	0.071–0.375 (0.226)

¹ The number in parentheses represents the average value, and the number range outside parentheses represents the minimum and maximum.

). Among domestic wild strains, the average minor allele frequency (MAF) of Sichuan wild strains was notably low at 0.202, with predominant alleles more widely represented. The average Nei diversity index for this group was only 0.303. In contrast, wild strains from Liaoning and Jilin exhibited higher average MAF values of 0.328 and 0.326, respectively, with Nei diversity index averages of 0.552 and 0.550, respectively. Further comparison showed that average heterozygosity values were notably high in Jilin wild strains (0.606) and Liaoning wild strains (0.598), along with their respective PIC values of 0.326 and 0.327, indicating significantly greater genetic diversity in Northeast China compared to Sichuan. Wild strains from Yunnan and Hubei displayed lower average PIC values of 0.260 and 0.281, respectively. The genetic characteristics of CBS strains were similar to those of wild strains from Yunnan, with relatively lower diversity among the tested strains. Cultivated strains exhibited an average MAF of 0.207, an average Nei diversity index of 0.287, observed heterozygosity of 0.435, and an average PIC of 0.226. These values indicated lower genetic diversity in cultivated strains compared to wild strains.

3.2.3. Population Structure

Based on the 181,731 high-quality SNP markers, Admixture software was used to analyze the population structure of the samples. The analysis considered K values ranging from 1 to 10 to determine optimal clustering. At K = 2, the cross-validation error rate was minimized at 0.514, categorizing the 47 tested strains into two main clusters, with 17 strains showing ambiguous classification (Figure 4, Table S2). The first cluster included 11 cultivated strains previously identified along with CBS3, while the second cluster predominantly consisted of 10 wild strains from Sichuan, 5 from Northeast China, and 4 from Hubei. This analysis underscored the genetic differentiation between cultivated varieties and specific wild strains, highlighting substantial genetic divergence.

Additionally, at K = 4, the cross-validation error rate was lower than at K = 3, with a rate of 0.548, and the strains were distinctly divided into four clusters. Compared to clustering at K = 2, strains CBS1, CBS2, CBS5, and CBS4 formed a separate cluster, while ZP1 and WSC were also categorized into their own cluster.

A total of 369 candidate core SNP markers were identified, facilitating the construction of a DNA fingerprint map to differentiate the 47 P. pulmonarius strains (Figure 5, Table S3). The fingerprint map clearly depicted the close genetic relationships and divergence among the strains. The minor allele frequency was reported as 0.416 (Table S4). These core SNPs exhibited high quality, even distribution across the genome, and strong discriminatory power among strains. The map illustrated that cultivated strains had high genetic similarity with minimal SNP variation, whereas wild strains exhibited greater diversity with numerous variations, indicating substantial genetic differences among wild strains.

3.3. Evaluation of Agronomic Traits in P. pulmonarius

A total of 42 P. pulmonarius strains were capable of producing fruiting bodies (Table S5). The frequency distribution indicated continuous variation for nearly all seven agronomic traits, suggesting that these traits are quantitative and controlled by multiple genes (Figure 6). Since the period from inoculation to primordia generation (PPG) represents the earliest time point for primordia formation, it was excluded from the ANOVAs. The ANOVA results revealed significant variation among strains for the remaining six traits (α = 0.05) (Table S6).

The boxplot clearly illustrates the variation for each trait (Figure 7). The most pronounced differences among strains were observed for the trait YFB, while PPG exhibited the least variation. However, it is worth noting that there were four significant outliers for PPG, indicating that although most strains showed a similar period from inoculation to primordia formation, a few strains had notably different times.

The average period from inoculation to primordium generation (PPG) was approximately 35 days (

Table 3. P. pulmonarius phenotypic descriptive statistical table.

Trait	SamNum 1	Mean	Sd 2	Median	Min 3	Max 4	Range	CV 5
CL	42	30.75	5.06	31.36	20.32	47.85	27.53	16.45%
CW	42	33.81	6.33	33.29	23.29	50.99	27.70	18.72%
SD	42	13.80	4.01	13.33	5.50	25.56	20.07	29.05%
SL	42	27.47	7.47	27.03	10.73	45.77	35.04	27.18%
PPG	42	34.86	6.54	33	31	64	33	18.76%
PIH	42	42.40	10.37	38	36	79	43	24.45%
YFB	42	15.30	11.06	10.41	2.15	42.76	40.61	72.32%

¹ SamNum: number of samples after removing missing values; ² Sd: standard deviation; ³ Min: minimum; ⁴ Max: maximum; ⁵ CV: coefficient of variation.

). Notably, strains ZP3, ZP4, WHLJ, and CBS5 exhibited faster primordium formation, occurring in as few as 31 days, while ZP2 required a longer period of 64 days. The broad-sense heritability for stipe diameter (SD) and yield of fruiting bodies (YFB) was relatively high, indicating that the variations observed for these traits were primarily influenced by genetic factors. The average period from inoculation to harvest (PIH) was around 42 days; however, strains ZP4, WHLJ, and CBS3 required only 36 days, while WSC5 took 79 days to reach harvest.

Cap size varied significantly among the 42 strains. For example, strain CBS3 had a cap length (CL) of 47.85 mm and a cap width (CW) of 50.99 mm, while strain WJL1 had a CL of only 20.32 mm and a CW of 23.29 mm. Stipe length (SL) also varied greatly, with WSC7 having the shortest length at 10.73 mm and CBS5 the longest at 45.77 mm. Stipe diameter (SD) showed notable variation as well, with WSC5 having the smallest diameter at 5.5 mm, while ZP9 had the largest diameter at 25.56 mm. In terms of YFB, strain ZP3 achieved the highest fruiting body weight at 42.76 g per bottle. Cultivated strains generally produced higher fruiting body weights with more considerable variation, as indicated by a coefficient of variation of 72.32%.

Significant differences were detected for CW, SD, and YFB among the groups. Strains in groups ZP, CBS, and WHLJ exhibited significantly larger yields and cap widths, while strains in groups WHB, WJL, and WLN showed significantly larger stipe diameters compared to other groups (Figure 8). No significant differences were observed for the other four traits among the groups.

Correlation analysis revealed a significant positive correlation between PPG and PIH (correlation coefficient = 0.77). Cap length and width were also significantly positively correlated (correlation coefficient = 0.91). Additionally, YFB showed a positive correlation with CL, CW, and SL, decreasing in strength in that order. Conversely, YFB exhibited a negative correlation with SD (Figure 9, Table S7). The negative correlation between YFB and SD suggests that a thicker stipe may adversely affect overall quality. Furthermore, harvesting time was negatively correlated with cap length and width, indicating that later fruiting resulted in smaller cap sizes.

3.4. GWAS Analysis Results

In this GWAS, significant SNPs across various scaffolds were identified for key traits (Figure 10,

Table 4. Proportion of variance explained (PVE%) and minor allele frequency (MAF) per significant SNP detected in whole population sampling.

Trait	Scaffold	p-Value	−log10p	Allele	MAF 1	SNP ID	PVE 2
CL 3	SJKF01000008.1	6.03 × 10−6	5.22	A/G	0.39	rs_474319	0.02%
CW	SJKF01000016.1	3.92 × 10−6	5.41	A/C	0.41	rs_338157	40.45%
	SJKF01000008.1	3.36 × 10−6	5.47	A/G	0.39	rs_474319	0.38%
	SJKF01000011.1	6.11 × 10−6	5.21	T/C	0.22	rs_949638	0.30%
	SJKF01000010.1	2.40 × 10−6	5.62	C/T	0.41	rs_1934108	0.43%
	SJKF01000006.1	3.44 × 10−7	6.46	T/A	0.17	rs_3071735	47.30%
	SJKF01000003.1	8.07 × 10−6	5.09	T/C	0.2	rs_3643242	48.76%
SD	SJKF01000005.1	2.34 × 10−6	5.63	A/G	0.21	rs_3145491	4.29%
SL	SJKF01000013.1	7.60 × 10−6	5.12	G/T	0.28	rs_459476	19.54%
	SJKF01000007.1	3.04 × 10−6	5.52	G/A	0.13	rs_499031	28.29%
	SJKF01000004.1	1.14 × 10−6	5.94	C/T	0.12	rs_1310326	20.78%
	SJKF01000004.1	7.99 × 10−6	5.1	T/C	0.43	rs_2224699	4.55%
	SJKF01000005.1	9.00 × 10−6	5.05	G/A	0.27	rs_2285326	32.72%
YFB	SJKF01000008.1	7.24 × 10−6	5.14	T/G	0.17	rs_885485	0.64%
	SJKF01000002.1	5.17 × 10−6	5.29	A/C	0.27	rs_4651725	65.72%

¹ MAF: minor allele frequency; ² PVE: phenotypic variance explained; ³ cap length (CL); cap width (CW); stipe diameter (SD); stipe length (SL); yield of fruiting body per bottle yield per bottle (YFB).

). On scaffold SJKF01000008.1, SNP rs_474319 explained only 0.02% of the variance in cap length (CL), indicating a minimal effect. Conversely, SNP rs_338157 on scaffold SJKF01000016.1 had a strong influence on cap width (CW) with a proportion of variance explained (PVE) of 40.45%. Scaffolds SJKF01000006.1 and SJKF01000003.1 also showed high PVEs for CL with SNPs rs_3071735 (47.30%) and rs_3643242 (48.76%). For stipe diameter (SD), SNP rs_3145491 on scaffold SJKF01000005.1 showed a smaller effect (PVE 4.29%). In stipe length (SL), SNPs rs_499031 and rs_2285326 on scaffolds SJKF01000007.1 and SJKF01000005.1, respectively, made significant contributions with PVEs of 28.29% and 32.72%. The highest impact was observed on scaffold SJKF01000002.1, where SNP rs_4651725 influenced yield with a PVE of 65.72%.

To avoid false positives, the long PPG periods of four strains (ZP1, ZP2, WCS5, and WYN1) and the PIH periods of five strains (WSC5, ZP2, WYN1, WJL1, WLN1) were excluded as outliers in the GWASs. After their exclusion, no significant SNP loci were identified for traits PPG and PIH.

4. Discussion

4.1. Genetic Diversity Analysis of P. pulmonarius Based on SNP

Traditional germplasm identification methods are characterized by long identification cycles, high costs, and susceptibility to environmental factors. Therefore, DNA molecular markers are crucial for germplasm identification at the molecular level. Compared to other molecular markers, SNP-based diversity analysis provides a more accurate and detailed interpretation of phenotypic differences at a lower cost, making it the most widely applied molecular marker [34,35]. SNPs, which represent the most common types of nucleotide sequence polymorphisms—including insertions and deletions—are abundant throughout the genome. These co-dominant markers are highly reproducible and suitable for genetic diversity and population structure mapping [36,37,38,39,40]. Previous studies have shown that SNPs are closely associated with individual traits, where different SNP positions may affect gene expression and phenotype through various mechanisms. The functional impact depends on whether the nucleotide variation affects protein structure or function.

Despite its relatively short cultivation history in China, P. pulmonarius holds promising economic potential and market prospects. Assessing genetic diversity in germplasm resources is essential for resource identification, new gene discovery, and breeding of new varieties. Current challenges in P. pulmonarius breeding include issues such as misidentification of varieties, insufficient research on breeding practices, and limited depth in genetic studies. Additionally, some P. pulmonarius varieties suffer from degeneration and poor disease resistance after several years of introduction. Understanding the genetic relationships and population structure of P. pulmonarius strains from different regions of China at the genomic level is essential for addressing these challenges.

To advance the breeding and application of P. pulmonarius germplasm resources, this study conducted whole-genome resequencing of 47 strains collected nationwide. A total of 181,731 high-quality SNP loci were identified after filtering, indicating substantial genetic variation among the strains. The significant number and proportion of SNPs make these strains well suited for subsequent genetic evolution andGWAS. These SNPs were used for DNA identity verification, phylogenetic relationship analysis, and genetic diversity assessment. Phylogenetic analysis revealed that most currently cultivated P. pulmonarius germplasm originates from the same source, exhibiting close genetic relationships and potential synonymy issues. In contrast, wild germplasm showed greater genetic diversity. Strains from the same geographic origin were grouped into the same cluster, a pattern also observed in other mushroom species [41]. The results of principal component analysis and population structure analysis were consistent with the phylogenetic tree, underscoring the reliability of SNP markers in assessing genetic relatedness in P. pulmonarius strains.

DNA fingerprinting, based on molecular markers or specific sequences, was used to assess genetic diversity and effectively differentiate strains. DNA fingerprinting has been widely applied in genetic diversity analyses of plants and microbes [42], providing valuable tools and technical support for germplasm identification and SNP-based fingerprinting methods [43]. In this study, a fingerprint atlas was constructed using 369 core SNP loci, offering unique genetic identification for each strain. These results provide a foundation for studies on genetic diversity, quality control, and molecular breeding.

4.2. Evaluation of Agronomic Traits

P. pulmonarius, one of the most widely cultivated edible fungi in China, faces challenges due to chaotic germplasm and limited genetic diversity in cultivars. These limitations stem from tissue separation, monosporal and polysporal hybridization, long-term artificial selection, and strain introductions across cultivation areas. These factors complicate the effective utilization and innovation of P. pulmonarius germplasm resources, restricting available breeding materials. Furthermore, unlike green plants, edible fungi exhibit lower morphological differentiation, and oyster mushroom fruiting body morphology is highly influenced by environmental conditions. Therefore, optimizing multiple traits is essential to maximize productivity when selecting superior varieties, necessitating the evaluation of germplasm based on multiple traits rather than a single trait [44].

In this study, we conducted comprehensive statistical analyses on seven key traits related to mushroom production in 47 P. pulmonarius strains. According to the ANOVA results, all seven traits were identified as quantitative traits primarily determined by genetic variation. The high broad-sense heritability values indicated that genetic factors are the primary contributors to the observed variation. Correlation analyses revealed significant associations, such as between the period of primordia generation and the period from inoculation to harvesting, as well as between cap length and cap width. Strains from Groups ZP, CBS, and WHLJ generally exhibited superior agronomic traits compared to others, suggesting their potential as breeding materials. The outliers observed for the traits PPG and PIH may be attributed to the low mycelial activity of certain strains, underscoring the need for improved strain preservation to ensure a robust analysis and evaluation system for P. pulmonarius germplasm. These findings provide valuable insights for optimizing breeding and cultivation strategies for P. pulmonarius.

4.3. Genome-Wide Association Study

GWASs rely on a comprehensive set of variants covering the entire genome. Initially, these studies used low-throughput markers like RFLP, AFLP, and SSR. However, with advancements in technology, the focus has shifted to higher-throughput and more stable SNPs, which now dominate the field. A GWAS identifies specific alleles contributing to desirable traits, facilitating genetic improvement efforts. While GWASs have been extensively applied in plants and animals to uncover genes associated with agronomic traits, their use in edible fungi remains limited. In this study, 47 P. pulmonarius strains from diverse geographic origins, showing significant phenotypic variation in agronomic traits, were used for a GWAS. All seven traits were identified as polygenic, controlled by multiple genes. Effective GWAS models, such as the linear mixed model (LMM) provided by GEMMA software, helped mitigate false positives and identified numerous key SNPs and candidate genes associated with these traits.

For instance, SNP rs_474319 on scaffold SJKF01000008.1, though statistically significant, had a low proportion of variance explained (PVE) at 0.02% for cap length, indicating a minimal impact on this trait. In contrast, SNPs like rs_338157 on SJKF01000016.1, rs_3071735 on SJKF01000006.1, and rs_3643242 on SJKF01000003.1 exhibited high PVE values (ranging from 40.45% to 48.76%), underscoring their significant roles in the genetic control of cap width and cap length. Similarly, SNPs associated with stipe length, such as rs_499031 on scaffold SJKF01000007.1 and rs_2285326 on SJKF01000005.1, showed high PVEs (28.29% and 32.72%, respectively), highlighting their substantial influence on this trait. These findings suggest that the genetic architecture of traits like cap length, cap width, and stipe length is shaped by SNPs distributed across multiple scaffolds, reflecting the complex genetic regulatory networks underlying these phenotypic characteristics.

Additionally, scaffold SJKF01000005.1 exhibited pleiotropic effects on both stipe diameter (SD) and stipe length (SL), while scaffold SJKF01000008.1 showed pleiotropy on cap length (CL), cap width (CW), and yield of fruiting bodies (YFB). The observed phenotypic correlations were effectively explained by genetic correlations, providing valuable insights for future breeding programs and molecular improvement strategies. Notably, SNP rs_4651725, associated with the YFB trait, exhibited a high PVE of 65.72% despite a relatively modest -log10P value of 5.29, making it less prominent in the Manhattan plot. This discrepancy may arise because the p-value primarily reflects statistical significance, which can be influenced by limited sample size, discontinuous phenotypic expression from the micro-cultivation system, and lack of strong natural selection, potentially causing the SNP to appear less significant. Additionally, the large PVE value might result from the high coefficient of variation of YFB among strains. The PVE quantifies the proportion of phenotypic variance explained by the SNP, highlighting its contribution to the traits. Thus, an SNP that marginally exceeds the significance threshold may still demonstrate a substantial PVE.

After removing outliers, no significant SNP associations were found for the period from inoculation to harvest or the primordium formation period. To detect significant SNPs associated with these traits, future experiments should increase sample size and conduct repeated trials to address outliers.

To date, GWASs have been used in various species, including staple crops like rice [45], model organisms such as Arabidopsis thaliana [46], and other plants like sorghum [47]. Their combination with genomics, transcriptomics, and metabolomics has provided a deeper understanding of species biology and genetics. This multidisciplinary approach helps clarify the genetic basis of complex traits. Looking ahead, GWAS technology is expected to continue evolving, with promising applications in the field of edible fungi. Future developments are likely to enhance its effectiveness in understanding genetic structures, optimizing breeding strategies, and contributing to the development of improved fungal varieties.

5. Conclusions

In this study, sequencing of 47 P. pulmonarius strains identified 4,430,948 SNPs. Analysis using 181,731 SNP markers revealed six distinct subpopulations, and a DNA fingerprint was developed using 369 core markers. Genetic diversity analysis showed minimal variation within cultivated strains from the same area but significant diversity among wild strains from different regions. GWASs identified key SNPs and candidate genes associated with five phenotypic traits, providing valuable insights for targeted breeding of improved edible fungus strains.