Development of Polymorphic Simple Sequence Repeat Markers using High-Throughput Sequencing in Button Mushroom (Agaricus bisporus).

The white button mushroom (Agaricus bisporus) is one of the most widely cultivated species of edible mushroom. Despite its economic importance, relatively little is known about the genetic diversity of this species. Illumina paired-end sequencing produced 43,871,558 clean reads and 69,174 contigs were generated from five offspring. These contigs were subsequently assembled into 57,594 unigenes. The unigenes were annotated with reference genome in which 6,559 unigenes were associated with clusters, indicating orthologous genes. Gene ontology classification assigned many unigenes. Based on genome data of the five offspring, 44 polymorphic simple sequence repeat (SSR) markers were developed. The major allele frequency ranged from 0.42 to 0.92. The number of genotypes and the number of alleles ranged from 1 to 4, and from 2 to 4, respectively. The observed heterozygosity and the expected heterozygosity ranged from 0.00 to 1.00, and from 0.15 to 0.64, respectively. The polymorphic information content value ranged from 0.14 to 0.57. The genetic distances and UPGMA clustering discriminated offspring strains. The SSR markers developed in this study can be applied in polymorphism analyses of button mushroom and for cultivar discrimination.

Introduction

The white button mushroom (Agaricus bisporus) belongs to genus Agaricus, family Agaricaeae, order Agaricales [1]. It is a plant biomass degrading fungus with a wide geographical distribution [2]. Apart from its ecological role, A. bisporus is considered as a healthy functional food due to the high contents of polyphenols, vitamins, minerals, polysaccharides, and proteins [3]. The button mushroom is one of the most widely cultivated edible mushrooms in the world and is the basis of a multibillion dollar industry [4,5]. Its cultivation began in the early 17th century in France [6], however, since 1950s, the production of this mushroom has been continued to increase worldwide and now represents 30% of the global mushroom production [7]. Currently, North America, Europe, India, and China are major mushroom cultivating regions [8], and cultivation efforts have recently increased in China and Korea [9]. Developing new cultivars is essential, as the mushroom production is limited owing to their complex life cycle patterns.

A. bisporus has a bipolar mating system [10], which can vary depending on varieties [4]. The A. bisporus var. burnettii and A. bisporus var. eurotetrasporus varieties produce spores with single nucleus. This enables breeding of new cultivars through conventional mating. The wild varieties have been used to generate commercial lines. However, A. bisporus var. bisporus generally shows better quality than A. bisporus var. burnettii in the commercial lines [11]. The cultivation of the button mushroom is dominated by a single hybrid strain produced from cross breeding of the Horst U1 and U3 strains developed early 1980s in the Netherlands [12]. Currently, most hybrid varieties are genetically similar to Horst U1, which has resulted in very low diversity. This lack of diversity makes it difficult to discriminate between cultivars.

Many countries have signed the treaty with the International Union for the Protection of New Varieties of Plants and have implemented laws related to breeders’ rights [13]. Moreover, the demand for royalties on the use of cultivated varieties is also expected to increase [14]. Therefore, there is an urgent need to develop a precise discrimination method for cultivars for the protection of breeder’s rights.

The morphological and biochemical characteristics of plants are the main basis of cultivar discrimination [15]. Similarly, in mushrooms, morphological characteristics that include color, shape, size, attachment, ornamentation, and other developmental process of fruiting bodies have been used [16]. However, these methods have limitations in their ability to discriminate closely-related cultivars, especially concerning environmentally sensitive traits [17]. The use of DNA-based molecular markers could overcome this problem, which makes it possible to evaluate genetic characters more easily and accurately than using morphological characters [18]. In this regard, molecular techniques that include restriction fragment length polymorphism (RFLP), randomly amplified polymorphism DNA (RAPD), amplified fragment length polymorphism (AFLP), single nucleotide polymorphism (SNP), and simple sequence repeat (SSR) have been widely used in the field of breeding for the evaluation of genetic diversity, introduction and selection of genetic resources, and discrimination of varieties [19-23].

Among the molecular markers, SSR microsatellite markers can be used to identify polymorphisms due to differences in the frequency of repetition in a genome. Moreover, SSR markers are highly polymorphic and thus have been widely used to evaluate genetic diversity and relationships [24]. SSR markers can be developed from genomic or expressed sequence tag (EST) libraries [25]. However, the genomic SSRs distributed throughout the genome have a great advantage over EST SSRs, which are only distributed in the transcribed region [26].

The present study aimed to develop genomic SSR markers to discriminate parental strains of button mushroom cultivars, which may provide a sound basis for button mushroom breeding programs.

Materials and methods

Strains and sample collection

Six strains of A. bisporus, including S1038-211, S1346-15, S1346-17, S1346-20, S1346-26, and S1346-110, were obtained from the National Institute of Horticultural and Herbal Science (http://www.nihhs.go.kr/). These strains originated from two parental strains, ASI 1038 and ASI 1346, which have been frequently used for the development of new cultivars in Korea [27]. The mycelia of each strain were cultured on CDA medium in the dark at 25 °C for a month.

DNA extraction and sequencing

The extraction of genomic DNA from the cultured mycelia was performed using a GenEX Plant Kit (GeneAll Biotechnology Co. Ltd., Seoul, Korea) following the manufacturer’s instructions. The extracted genomic DNA was subjected to high-throughput methods of re-sequencing. Illumina paired-end DNA library (average insert size of 500 bp) was constructed using the Illumina TruSeq library preparation kit following the manufacturer’s instructions. The libraries were sequenced with 2 × 300 bp on the Illumina MiSeq platform at LabGenomics (http://www.Lab.genomics.com/kor/).

Genome assembly

The genomic reads generated by the Illumina MiSeq platform were quality-trimmed with Trimmomatic ver. 0.33 to obtain reads with high quality (Phred score >20). These were de novo assembled into contigs using Platanus ver. 1.2.1 with default parameters [28]. Scaffolds were generated from contigs using Platanus ver. 1.2.1 with default parameters. Those exceeding 500 bp were selected for further study. De novo assemblies of five samples were performed independently to generate contigs and scaffolds sets for each sample.

Gene prediction and functional annotation

Genes in scaffold sequences were predicted based on a self-training algorithm using GeneMark-ES ver 4.32 [29]. Deduced protein sequences of five samples were compared with each other based on the similarity with protein sequences from a reference genome, A. bisporus (GenBank Accession No: GCF_000300575.1). BLASTP searches with cut-off E-value of 1e−4 were employed to investigate similarity. Protein sequences specific or common to each of five samples were searched and depicted by a Venn diagram. Five gene sets of five samples and the reference genome were compared using BLAST ver. 2.3.31 and then specific genes were selected for gene ontology (GO) analysis using Blast2GO [30]. GO terms assigned to three categories—Molecular function, Biological process, and Cellular component—were investigated and the frequency of GO terms was compared among the five gene sets.

SSR identification and primer design

MIcroSAtellite software was used as the identification tool to detect the SSR markers present in the cp genome (http://pgrc.ipk-gatersleben.de/misa/). This software allows the localization and identification of both perfect and compound microsatellites with 1 to 6 nucleotides in the basic repeat unit. We selected the 1,858 SSR motif, which showed the largest diversity among all SSR motifs, and selected a total of 91 SSR markers. Primer design parameters were set as follows: length range, 18–23 nucleotides with 21 as the optimum; PCR product size range, 200–300 bp; optimum annealing temperature (Ta), 55 °C; and GC content 50–60%, with 51% as the optimum.

Validation of polymorphic SSR markers

The selected SSR markers were used for the identification of genetic diversity in the six A. bisporus strains. The DNA extracted from these strains for PCR templates was quantified with a model K5600 micro-spectrophotometer (Shanghai Biotechnology Co. Ltd., Shanghai, China) and diluted to 20 ng/μL. PCR was performed with 20 ng of template DNA using 2× PCR Master Mix Solution (i-Taq, iNtRON Biotechnology Inc., Seongnam, Korea) according to the manufacturer’s recommendations. Final volume and concentration of primers in all reactions were fixed to 20 µL and 10 µM, respectively. PCR reactions were performed as follows: 95 °C for 3 min; 35 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s; and finally 72 °C for 20 min. The size of the PCR product was confirmed using a fragment analyzer (Advanced Analytical Technologies Inc., Ankeny, IA, USA). Allele size was scored with PRO Size® 2.0 software (Advanced Analytical Technologies Inc.). To calculate the genetic diversity of each SSR locus, the major allele frequency (MAF), number of genotypes (NG), number of alleles (NA), observed heterozygosity (HO), expected heterozygosity (HE), and polymorphic information content (PIC) values were calculated with Power Marker V3.25. The genetic distance between the strains were calculated by the Shared Allele method and a dendrogram was constructed using the unweighted pair-group mean algorithm (UPGMA) method.

Results and discussion

Five different Agaricus offspring were sequenced with the Illumina MiSeq system, which produced 6,358,118–9,337,246 paired-end raw reads. A general summary of the statistical data on sequencing and assembly is provided in Table 1. The genomic reads were quality-trimmed with Trimmomatic ver. 0.33. The total number of valid reads after quality trimming (Phred score >20) was 5,980,422–8,764,486. Then the quality-trimmed reads were de novo assembled separately using Platanus ver. 1.2.1, which produced 10,676–11,841 completely assembled contigs. The overall sequence assembly produced 2,497–3,614 scaffolds ensuring the rules of minimum length requirement (500 bp).

Table 1.

Summary of Illumina paired-end sequencing and assembly for A. bisporus used in this study.

Sample	Raw Reads	Cleaned Reads	Contigs	N50	Scaffolds	N50	Predicted Gene	Avg Gene Length
S1038-211	7,338,790	7,080,774	11,414	5,019	3,226	30,120	9875	1394
S1346-110	9,337,246	8,764,486	11,268	5,306	2,681	42,867	9627	1448
S1346-15	6,440,746	6,089,334	11,841	4,112	3,316	25,673	9518	1395
S1346-17	6,358,118	5,980,422	12,374	3,739	3,614	20,952	9323	1395
S1346-20	8,107,872	7,620,128	10,676	5,413	2,497	49,431	9631	1446
S1346-26	8,909,410	8,336,414	11,601	4,811	2,803	38,055	9620	1433

In total, 9,323–9,875 gene models were predicted, with average size of 1,394–1,448 bp among the five samples with GeneMark-ES ver 4.32 software program (Table 1). The deduced protein sequences were compared with the A. bisporus reference genome. The BLASTP searches with a cut-off E-value of 1e−4 indicated that the annotated genes were highly homologous with A. bisporus in the database. The Venn diagram (Figure 1) illustrated the interrelation of Blast hits for predicted unigenes against the A. bisporus reference genome database. The results revealed that 6,559 unigenes could be considered as orthologous genes across all five samples, as shown in the Venn diagram.

Figure 1.

Venn diagram of Blast hits for unigenes against the reference genome, Agaricus bisporus.

To obtain better understanding on the functions of the predicted unigenes, GO classification was applied (Figure 2). The GO classifications with the reference genome assigned many unigenes to a wide range of GO categories. Different numbers of annotated unigenes were distributed under different GO categories, with 1,554 unigenes for biological processes, 722 unigenes for cellular components, and 454 unigenes for molecular functions, which indicated the high accuracy of the annotation. However, GO classification with five samples revealed only 221–552 unigenes related to different GO classes, which was markedly less compared to the reference genome. The GO annotation analyses revealed 90–302 unigenes for biological processes, 54–174 unigenes for cellular components and 20–76 unigenes for molecular functions, which indicated the reduced number of unigene annotations.

Figure 2.

Gontology classification of predicted unigenes.

The frequency distribution of the SSRs depending on the repeat type is shown in Figure 3. According to the motif type, the SSR motifs included 3 bp with 898 (48%), 2 bp with 527 (28%), 4 bp with 221 (12%), 6 bp with 123 (7%) and 5 bp with 89 (5%). Among the 91 selected SSR markers, 43 were polymorphic and the major allele frequency (MAF) ranged from 0.42–0.92, with an average of 0.711. The number of genotypes (NG) ranged from 1 - 4 with an average of 2.3, whereas the number of alleles (NA) ranged from 2–4 with an average of 2.5. The observed heterozygosity (HO), expected heterozygosity (HE), and polymorphic information content (PIC) averaged 0.291, 0.395, and 0.340, respectively (Table 2). The UPGMA clustering showed that the six button mushroom strains could be divided into two clusters when analyzed with the 43 SSR markers. All offspring could be discriminated from the parental cultivars (Figure 4).

Figure 3.

Percentage distribution of different repeat types of SSRs markers in the button mushroom strains based on resequencing data.

Table 2.

Diversity statistics of the 43 polymorphic simple sequence repeats (SSR) markers in A. bisporus used in this study.

Locus	SSR motif	Primer Sequence (5’-3’)	MAF	NG	NA	HO	HE	PIC
AB-gSSR-0113	(AG)6	F-TACTCAACAACATGATGCCATT	0.58	3	3	0.5	0.57	0.5
		R-CTTATATTCTCCTCGCTGTTGG
AB-gSSR-0126	(CT)6	F-ATCATCGATTAGTCCTTCAACG	0.83	2	2	0	0.28	0.24
		R-TATCAGTGACTGGTTGCAGAAG
AB-gSSR-0139	(AG)6	F-GAATACGTAAGCTGACCTCCTG	0.67	3	3	0	0.5	0.45
		R-CAATTCCTGCTTTCTCTCTACG
AB-gSSR-0150	(CT)6	F-CTCTAACAGCGTCTAACAGCCT	0.83	2	2	0	0.28	0.24
		R-TATTCAATTTCATGCTCAGTGG
AB-gSSR-0179	(CT)6	F-CTATTACTTCCTCGCTCACCAC	0.5	3	4	0.83	0.63	0.56
		R-CTTACGAAATAGAACGAGGGTG
AB-gSSR-0182	(CT)6	F-CATCGAAGCGATATAAGAAAGG	0.83	2	2	0	0.28	0.24
		R-TAGTGTTGTTGTATTTCGGCTG
AB-gSSR-0193	(CT)7	F-ATGTACTCTCACCTGCCTCACT	0.58	3	2	0.5	0.49	0.37
		R-AGGACGACTTCGATGAGTTAAA
AB-gSSR-0199	(AG)9	F-TCAAATTCCATGGTGTTGAGTA	0.42	2	4	1	0.64	0.57
		R-GGACATGCGTATAGACCTCATT
AB-gSSR-0218	(AC)6	F-GAACGTATGGGATGTCAAGACT	0.67	2	2	0	0.44	0.35
		R-AATAGAAGATTCCTCTCCTCCG
AB-gSSR-0222	(AC)6	F-TTGAAGGAGGCTTAGAGTTCTG	0.83	2	2	0	0.28	0.24
		R-ATTGTTGATCACAAGGAAGAGG
AB-gSSR-0270	(TA)6	F-TCGAAATCGCTTAACTTGGTAT	0.83	2	2	0	0.28	0.24
		R-TATGTCTGGACGTTTCCTCTTT
AB-gSSR-0296	(TA)6	F-AGATTATGGAGATTCATGGCAC	0.92	2	2	0.17	0.15	0.14
		R-AAGAACCAGTTTAAGAATGGCA
AB-gSSR-0298	(TA)6	F-TGGATACGCATCTTTCATGTTA	0.67	3	3	0	0.5	0.45
		R-TCTGAAATCGAGTGCTGATATG
AB-gSSR-0305	(TA)6	F-GAGCGATGGGAACAGAATATAA	0.83	2	2	0	0.28	0.24
		R-ACAAGGTTATACACACCGAAGG
AB-gSSR-0310	(TA)6	F-CTGATAATCTCCAATCGTCCAT	0.83	2	2	0	0.28	0.24
		R-TTACTCGGAGTTGTTTGTTGAA
AB-gSSR-0391	(TC)6	F-AGGTGGTATGATAGTCCTGCAC	0.83	2	2	0	0.28	0.24
		R-CACTAACAATCGTGAGCAGTTC
AB-gSSR-0396	(TC)6	F-AGCAACTGCTTACCTCCATTAG	0.83	2	2	0	0.28	0.24
		R-TAGCCAATGGGTTCACCTATAC
AB-gSSR-0428	(TC)6	F-CCCTCAACTGTATTCAACCATT	0.83	2	2	0	0.28	0.24
		R-TCAACAAGTCCAACAACATCAT
AB-gSSR-0432	(TC)6	F-TTGCTAGTAAAGCTCAGCAACA	0.5	2	3	1	0.61	0.54
		R-GTCTTTCGAATGAAGACCTTTG
AB-gSSR-0444	(TC)6	F-CGGGTCTCAATAAGTTCAAGAG	0.83	2	2	0	0.28	0.24
		R-GGTGGTTGCATTATTATTCGAT
AB-gSSR-0465	(TC)6	F-TGTCTTTGGCTCTAGGACGTAT	0.83	2	2	0	0.28	0.24
		R-CATTATTAACTCTAGGCGTCGG
AB-gSSR-0474	(GA)6	F-GACGATGGCACCATATAAATCT	0.5	3	4	1	0.63	0.56
		R-TACCGAAGTAGGAGTTGGAGAA
AB-gSSR-0488	(TC)7	F-AATTCATGAGAAAGAAACGAGC	0.5	1	2	1	0.5	0.38
		R-ACTCTTTAAATTGTGGCCTTCA
AB-gSSR-0489	(GA)6	F-CAATACTCAACGGATTCGACTT	0.5	2	3	1	0.57	0.48
		R-CTCATCCTCATCAACAGCATAA
AB-gSSR-0512	(TC)6	F-TGTTGGATTAGAAGGAGCAAGT	0.75	2	2	0.5	0.38	0.3
		R-GGTTAACCAATATTACGACGGA
AB-gSSR-0532	(ATT)5	F-TTTCATGTCGCAAATTGAGTAG	0.83	2	2	0	0.28	0.24
		R-ACCCTTAGGCTAGGTACGAATC
AB-gSSR-0562	(CTT)5	F-ATCCCTTCTTGACATGAACATC	0.83	2	2	0	0.28	0.24
		R-TAGGTGTCGACGTTAGACCTTT
AB-gSSR-0564	(CTT)6	F-GAGGAAAGTTGCTGGTATGAAG	0.92	2	2	0.17	0.15	0.14
		R-GAAGAAGACGACGAAGAAGAAA
AB-gSSR-0574	(AAG)6	F-ACACGTATCCGAGTAGAGGAGA	0.82	2	2	0.17	0.15	0.14
		R-AATTGCACAACCTCCAAATAAC
AB-gSSR-0580	(GTT)5	F-AAGGGTGGACTTAGCTGATAGG	0.75	3	3	0.17	0.4	0.36
		R-ACTCTCTCACCCAGCTCAAGAT
AB-gSSR-0582	(GTT)5	F-TATTCCCAACTGTAGTCGCTTT	0.83	2	2	0	0.28	0.24
		R-GCTGAACTTGGTATCGAATCTC
AB-gSSR-0584	(GTT)6	F-CTCGTCCTCGAAATAATACTCG	0.58	2	2	0.83	0.49	0.37
		R-ATTTGAGAGGGAAAGTTGAACA
AB-gSSR-0603	(GTT)5	F-GCGAAACACATCTAGGAAACAT	0.67	3	4	0.17	0.51	0.48
		R-ACGATGATGAACACGATAATGA
AB-gSSR-0608	(TAT)5	F-AAGAGACTCACGATTTCTTTGG	0.83	2	2	0	0.28	0.24
		R-CGAAGTTTAGCAAGGGAGTATG
AB-gSSR-0611	(ATG)5	F-CTACCCACGCTAGTTACCTTTG	0.92	2	2	0.17	0.15	0.14
		R-TAATTCAACCCTGTCGGTATTC
AB-gSSR-0630	(CAT)6	F-ACTCGTATTCAACTGGGAGCTA	0.5	3	3	0	0.61	0.54
		R-TCCAACATATGAGAGTGATGCT
AB-gSSR-0687	(GAT)5	F-ATTTGATGCTTCACGTTTCTTT	0.58	2	2	0.83	0.49	0.37
		R-ATCTTCATCGACGTTGTCTTCT
AB-gSSR-0709	(TCT)5	F-GGATCGCCCTAATACACAATTA	0.42	2	4	1	0.64	0.57
		R-CAGTGGTAGTTGATCGCAGTAA
AB-gSSR-0713	(AGA)7	F-AATCATGGTTATTCAGGACTCG	0.5	2	3	1	0.57	0.48
		R-CTGGAAGACTGTTATCAGAGGG
AB-gSSR-0730	(TCT)6	F-AGAGATGAGGGAGGAATAGAGG	0.5	4	3	0.33	0.61	0.54
		R-AACAGGACCTCATACAACAAGC
AB-gSSR-0736	(AGT)6	F-CATTAGCAGGTTTAGGAGATGC	0.67	3	3	0	0.5	0.45
		R-ACCACTACTACTCCACGACACC
AB-gSSR-0784	(GCT)5	F-ACCTCTACAAGATCGAACAGGA	0.67	3	3	0	0.5	0.45
		R-CACGATACGTAGGTTGTCCAC
AB-gSSR-0788	(AGC)5	F-CTGGAGAATAGGAGGAGGAGTT	0.92	2	2	0.17	0.15	0.14
		R-TCCAACAATGCTAGTGTCCATA
		Mean	0.711	2.3	2.5	0.291	0.395	0.34

MAF: major allele frequency; NG: number of genotypes; NA: number of alleles; Ho: observed heterozygosity; He: expected heterozygosity; PIC: Mean polymorphic information content.

Figure 4.

UPGMA dendrogram of six button mushroom strains constructed by using genetic distance analysis based on 43 SSR markers.

A. bisporus is one of the most cultivated edible mushrooms globally and remains an important component of the human diet. Different trials have been performed to breed the cultivars and have succeeded in the release of the first hybrid strains [31]. However, genetic manipulation by mating A. bisporus strains was limited due to the formation of secondary homothallic mushroom species during bisporic production of basidiospores [32]. Therefore, it is necessary to explore the correlation between the genotypes to effectively operate breeding programs.

In this study, we chose five different diploid strains for genome sequencing because their genetic relationships with other strains are very limited. The genome sequences of diploid spores were analyzed using the Illumina Miseq platform. The sequence assembly produced 2,497–3,614 scaffolds (>500 bp) from the 6,358,118–9,337,246 paired-end raw reads, which revealed only partial genome assembly. In general, the high percentage of repetitive sequences and high expansion of retro-transposon gene families are considered to be major obstacles in genome assembly of mushroom species [33]. However, the next generation sequencing technologies have efficiently solved this high repeat ratio problem. The GeneMark-ES program based on a self-training algorithm predicted 9,323–9,875 gene models with average size of 1,394–1,448 bp among the five samples. In addition, 6,559 unigenes were considered as orthologous genes across all five samples, as shown by Venn diagramming. Different number of unigenes in the Leucocalocybe mongolica genome have been against three different databases [34]. Similarly, in the GO classification, only 221–552 unigenes were classified to predicted functions.

Identification and development of genomic SSR loci using whole genomes have been successfully applied in several plants [35,36] and the mushroom [37]. In general, the wild type has higher genetic diversity than cultivars [38-41]. Therefore, the SSR markers developed using offspring in this study will reveal more diversity when applied to the wild types. The SSR markers successfully discriminated offspring from common parent and thus could be used in cultivar discrimination, molecular breeding, and studies of the genetic structure of the button mushroom.