One hundred strains of Cyclocybe chaxingu were collected and are maintained at the Jiangxi Agricultural University Culture Collection Center (JAUCC). These strains came from three different sources: private mushroom farms, culture collections and the wild. Geographically, these 100 strains came from 12 provinces and two municipalities in China. These strains were selected for this study to represent the broad geographic and ecological origins of this species distributions. The details of these strains are presented in Suppl. material 1.

Our mushroom fruiting protocol followed that described in (Gao et al. 2026). Briefly, to fruit these strains, their mycelia were first grown on PDA (per litre: potato infusion from 200 g potato, 20 g glucose, 20 g agar; natural pH) and incubated at 25 °C in the dark for 7 d. Fully colonised agar blocks were then transferred to sterilised solid substrate packed in polypropylene filter bags. The cultivation substrate (% w/w, wet basis) consisted of 15% cottonseed hulls, 20% coarse sawdust, 20% fine sawdust, 20% bagasse, 15% wheat bran, 5% soybean meal, 3% cornmeal, 1% gypsum and 1% lime, adjusted to ~ 65% moisture. During incubation, environmental conditions were 20–25 °C, 60–70% RH, with adequate fresh-air exchange and darkness. After ~ 50 d, when bags were fully colonised, they were opened to induce primordia and fruiting at 22 °C, ~ 95% RH and ~ 600 lx diffuse light.

To analyse strain genotypes, their mycelia were ground in liquid nitrogen and genomic DNA from each strain was extracted by the CTAB method (Porebski et al. 1997). Universal fungal primer pairs were used to amplify regions of the large subunit of the ribosomal RNA gene (LSU) and the translation elongation factor 1-α (TEF1-α). These two markers were chosen for this study for several reasons. The first was their readily available universal primers that could be used to obtain sequences from all of our strains. The second was that these two gene fragments have been used for fungal phylogenetics and the identification of putative cryptic species, including as secondary fungal DNA barcodes (Xu 2016). The third is that, at least for TEF1- α, intraspecific variations have been reported within many fungal species (Buyck and Hofstetter 2011). PCR products were verified by 1% agarose gel electrophoresis and Sanger sequenced by Sangon Biotech (Shanghai) Co., Ltd. (Yang and Qin 2023). Primers and cycling profiles are listed in Table 1.

Raw reads were inspected and assembled to consensus chromatograms; sequences were queried by NCBI BLAST for preliminary identification. Locus-wise alignments were generated with MAFFT (default/normal mode) (Katoh and Standley 2013), trimmed and concatenated in MEGA (Sudhir et al. 2018). Phylogenetic analyses were run in PhyloSuite (Zhang et al. 2020) with ModelFinder to select the best-fit partition models (Kalyaanamoorthy et al. 2017). MrBayes v.3.2.7a was used under a partitioned scheme (two parallel runs, 2,000,000 generations, sampling every 1,000; first 25% discarded as burn-in). Convergence was checked by the average standard deviation of split frequencies (< 0.01)

Ten agronomic traits were quantified on market-mature fruiting bodies (definitions in Table 2; measurement scheme in Fig. 1). Specifically, for each strain, six bags were prepared and inoculated with agar plugs of the same dimensions. All bags were cultivated under identical conditions for 55 days, after which they were opened and randomly arranged to initiate fruiting. Data for DHF, YLD, AFW and FBN were collected individually for each bag. For traits SL, SD, CD, CT and GW, ten representative fruiting bodies with standard shape and a range of sizes were chosen from each bag and measured (Table 2). Subsequently, the mean, standard deviation and other essential statistics were computed for every trait, as described below.

All analyses were performed in R. Descriptive statistics (means, standard deviations, distributions) were computed with dplyr. Pairwise associations amongst the ten agronomic traits were estimated using Pearson’s r and visualised as correlation heatmaps with corrplot and pheatmap (Taiyun and Viliam 2021). For all analyses that involved multiple comparisons — specifically those performed with the corrplot, pheatmap, dendextend, FactoMineR and factoextra packages — P-values were adjusted to minimise false discovery rate (FDR).

Before multivariate analyses, trait values were z-standardised. Phenotypic dissimilarity between strains was defined as 1 – r (Pearson correlation across traits). Strains were grouped by hierarchical clustering (UPGMA); dendrogram construction and visualisation were performed with dendextend (Galili 2015) and UPGMA routines (Moulton et al. 2018). The resulting tree depicts phenotypic similarity patterns amongst the 100 strains.

Global structure was summarised using principal component analysis (PCA) with FactoMineR and the scores and loadings were visualised with factoextra (Iezzoni and Pritts 1991; Lê et al. 2008; Kassambara and Mundt 2016), providing a compact view of trait covariation and strain grouping in C. chaxingu.

All strains produced fruiting bodies under the standardised cultivation conditions and ten traits were recorded for each of the 100 strains (Table 3). All 10 traits showed significant variation amongst strains (one-way ANOVA, p < 0.01). Coefficients of variation (CV) were highest for FBN (45.6%), followed by AFW (35.0%), MI (22.7%), GW (19.7%) and CT (19.7%), indicating pronounced phenotypic diversity in fruiting body number- and size-related traits within our strain collection. Moderate variations were also observed for the remaining five traits: SD (CV = 16.8%), YLD (14.6%), CD (14.2%) and SL (13.1%), while DFH varied the least (7.6%). Collectively, these metrics reveal substantial phenotypic diversity within C. chaxingu, providing abundant germplasm for downstream strain selection and breeding.

Correlation heatmaps (hierarchical ordering) revealed both positive and negative associations amongst traits (Fig. 2A). MI showed positive correlations with CD and CT and a negative correlation with SL, indicating that larger, thicker caps combined with shorter stipes yield higher MI values. SD correlated positively with AFW, CD and CT, but negatively with SL, suggesting that thicker stipes accompany heavier, larger fruiting bodies with shorter stipes. For yield-proximal patterns, YLD correlated positively with SD, AFW, CD, CT, GW (and weakly with DFH), but negatively with FBN, consistent with a trade-off between producing fewer, heavier versus many, lighter fruiting bodies. AFW was strongly and positively associated with SD, CD, CT and GW and negatively with FBN, reinforcing this pattern. The abbreviations in Fig. 2B and Fig. 2C are listed in Table 4.

Reporting note: In Fig. 2C, post-hoc contrasts are summarised as Tukey’s Honest Significant Difference (HSD) test letters (FDR-adjusted), which are shown above each box to replace multiple asterisk bars.

UPGMA clustering on standardised trait data grouped strains into four major phenotypic branches (Fig. 3). Basidiome colours largely tracked these branches and served as a practical discriminator; therefore, strains were partitioned into five colour-based groups: white (WH), black (BK), dark brown (DB), brown (BN) and yellow (YE) (Fig. 4). Between-group comparisons showed pervasive differences across traits (overall p < 0.001; Fig. 2C). Briefly, DFH was shortest in BN and longest in YE; YLD was highest in BK and lowest in YE; AFW tended to be higher in WH and lower in YE; FBN was highest in YE and lowest in WH. Significant differences were also observed in SL, SD, CD, CT, GW and MI, indicating multivariate differentiation amongst the five groups.

The first three principal components (PCs) accounted for 83.98% of the total variance (PC1 = 57.74%, PC2 = 15.87%, PC3 = 10.37%; Table 5). PC1 loaded strongly on CT, AFW, SD, CD, GW, MI (positive) and FBN (negative), representing a size/mass axis contrasting a few large fruiting bodies with many small ones. PC2 was dominated by SL (positive) and, to a lesser extent, MI (negative), capturing a stipe-elongation gradient. PC3 loaded mainly on DFH (positive) and YLD (positive), describing earliness/yield timing. The PC1–PC2 scatter (Fig. 1B) showed partially separated clouds for the five colour groups: YE tended to the lower-left, WH to the lower-right, BN to the upper-left and BK near the upper-right; DB lay intermediate between BN and BK, suggesting that colour phenotypes are associated with multivariate agronomic patterns.

To investigate the potential relationships between genetic and phenotypic differences amongst strains, we juxtaposed a Bayesian phylogeny from concatenated LSU + TEF1-α sequences with a phenotypic dendrogram derived from the ten agronomic traits (Fig. 5). The phylogeny resolved four principal genetic clusters, whereas phenotypic clustering defined five colour-coded groups (WH, BK, DB, BN, Y). WH and BK frequently co-clustered genetically and a subset of BN strains were mixed with DB/WH/BK. Overall, several phenotypic groupings were partially supported by the phylogeny, suggesting that some agronomic patterns reflect underlying genetic structure, while others may arise from phenotypic plasticity or polygenic variation beyond the resolution, based on the two sequenced loci.

Trait standardisation and DUS (Distinctness, Uniformity, Stability) guidelines are essential for variety protection and orderly industry growth. As fruiting and fruiting body features of edible mushrooms are highly sensitive to environmental conditions, DUS trials should minimsze environmental noise by testing candidate and reference varieties under identical conditions with adequate plot size and sample number to identify heritable variations; where feasible, use centralised test sites to improve accuracy and reproducibility (Chu et al. 2023). Under current Chinese practice, distinctness can be established when a test variety shows a repeatable, significant difference in ≥ 1 trait versus the reference. Based on our variance structure and between-strain contrasts, stipe diameter (SD), stipe length (SL) and gill width (GW) emerged as practical indicator traits for DUS testing in C. chaxingu.

Since the 1980s, C. chaxingu cultivation has expanded significantly in China (Liu et al. 2020; Gao et al. 2026). In our panel of 100 strains, all 10 traits showed significant variation amongst strains (ANOVA, p < 0.01), with the greatest trait value dispersion in FBN (CV 45.6%) and the second highest in AFW (CV 35.0%). This breadth of phenotypic diversity indicates rich germplasm diversity for selection and provides a foundation for breeding and variety protection (Li et al. 2024).

Correlation patterns support two contrasting fruiting body yield strategies: (i) fewer, heavier fruiting bodies (higher AFW; larger CD/CT/SD) and (ii) many, lighter fruiting bodies (higher FBN). YLD correlated positively with AFW, SD, CD, CT, GW and negatively with FBN, suggesting that increasing cap/stipe size while moderating FBN is a plausible route to higher yield. DFH is associated with CT and GW, so these can serve as auxiliary indicators when selecting for shorter cycles. As correlations do not imply causation, these directions should be validated in multi-environment fruiting trials and through selection and response-to-selection trials.

Our analyses identified strains with superior performances in one to multiple economic traits. These strains are candidate parental strains for future breeding programmes that target the specific sets of traits for improvements. For example, we identified strains with either short cultivation periods, high overall productivity. These strains could serve as parental stocks for breeding new strains with three combined high-quality traits. In addition, SD, SL and GW may also be employed as practical metrics for C. chaxingu in DUS testing, offering valuable reference points for the outcomes of cultivation experiments conducted under diverse conditions.

UPGMA clustering resolved four major phenotypic branches; inspection of basidiome colour supported five colour groups (WH, BK, DB, BN, YE) with pervasive differences across traits (overall p < 0.001). PCA indicated three interpretable axes: PC1 (57.74%) = size/mass (CT, AFW, SD, CD, GW vs. FBN); PC2 (15.87%) = stipe elongation (SL); PC3 (10.37%) = earliness/yield timing (DFH, YLD). Operationally, BK tended to have higher YLD, BN had shorter DFH and YE had many. but lighter fruiting bodies. Thus, fruiting body colour is a useful field cue, but because it can be influenced by environmental conditions, other quantitative trait descriptors should be included for more accurate identification of strain groupings.

Concatenated LSU + TEF1-α sequences could intuitively reflect the five groups defined by phenotypic clustering, particularly groups BN, DB and YE. Despite overall genetic co-clustering between WH and BK, most admixed individuals clustered separately on distinct branches. This partial congruence suggests that some agronomic differentiation reflects genetic structure, while other patterns may arise from phenotypic plasticity or polygenic backgrounds not captured by the two sequenced loci. Future work with genetic crosses and denser markers (e.g. multilocus sequence typing, amplified fragment length polymorphisms and whole-genome sequencing) and formal congruence statistics (Mantel, Procrustes, adjusted Rand index) will refine these links.

Our study has several limitations. First, we analysed the fruiting pattern at only one site under a single cultivation condition. Multi-site, multi-fruiting condition and multi-season trials are needed to quantify G×E and trait repeatability (ICC). DUS panels will benefit from BLUP-based means and effect sizes (η² or Cliff’s δ) and from the inclusion of biological efficiency (BE) and quality traits (texture, aroma) alongside yield. Second, only two loci were sequenced and both loci are highly conserved within a species. Analysing more loci will allow greater differentiations amongst strains, which may help identify potential linkages between genomic regions and phenotypic traits. The tetrapolar heterothallic system of C. chaxingu (Cai et al. 2007; Zheng et al. 2007) is ideal for investigating such relationships and can be harnessed in breeding schemes.