Identification of novel putative causative genes and genetic marker for male sterility in Japanese cedar (Cryptomeria japonica D.Don)

EST collection, sequencing and de novo assembly

Two sequencing platforms, the Roche 454 and Illumina HiSeq system, were used to sequence 19 EST libraries constructed from multiple organs (Additional file 1); the sequencing and assembly results are summarized in Additional file 2. The sequence data used for the following analysis, including already reported data (cambium during the active season (DRA000525) [34] and shoots (DRA001261) [35]), were from approximately 3 million reads with an average length of 437.8 bp, amounting to 1.7 Gbp. The variation in obtained isotig number is presumably related to the sequencing read number for each library.

All sequences generated from the nine libraries generated from the Roche 454 platform were assembled into 34,731 isotigs and these assembled isotigs were regarded as reference sequences. The mode and arithmetic average of isotig length were 1000 bp and 1724. 57 bp, respectively (Additional files 2 and 3). More than 99% of the isotigs were longer than 100 bp, and the median (N₅₀) was 2075 bp. The depth of reference sequences ranged from 1 to 364.99, with an average of 25.18 reads/site in each isotig (Additional file 2). The size distribution of the reference sequences was compared with several model plants such as Arabidopsis thaliana, Oryza sativa, Zea mays, Selaginella moellendorffii, Physcomitrella patens, Populus trichocarpa, P. abies, and P. taeda (Additional file 3). The length of the reference sequences of Japanese cedar was at least as long as the model plants, suggesting that most of the isotigs obtained include nearly full-length gene sequence information. Our reference sequences generated from these nine libraries were BLASTN searched against EST sequences in the Japanese cedar database (ForestGEN), resulting in 15,815 sequences with significant matches in ForestGEN (Fig. 1); more than half of the sequences (18,916 isotigs) were regarded as newly collected in Japanese cedar. Thus, our effort nearly doubled the number of EST resources in Japanese cedar.

All the reference sequences generated from the nine libraries on the Roche 454 platform were compared with isotigs of each library and classified into four categories of organs (male strobili, shoots, wood and roots) (Additional file 4). The number of homologous isotigs thought to be expressed in each organ was 8234 in male strobili, 19,882 in shoots, 19,942 in wood and 16,302 in roots (Additional file 4). Though 5486 isotigs (16.3%) were common across organs, 477 isotigs (1.4%) in strobili, 3277 isotigs (9.4%) in shoots, 3573 isotigs (10.3%) in wood, and 2458 isotigs (7.1%) in roots were organ specific (Additional file 4). These results influenced the number of reads obtained; however, these distinct isotigs presumably reflect organ-specific or common functions of genes expressed in each organ.

SNP discovery

In the process of SNP discovery, from a total of ~ 540.1 billion cleaned reads, a total of 573,795 SNPs were identified (Additional file 5) and 73,274 SNPs considered as robust for genotyping were selected as Axiom_Cj70K_ver. 1 (Additional file 1); from secondary SNP discovery, 89,993 SNPs were identified (Additional files 1 and 5).

In the Axiom genotyping assay for the set of reads called Axiom_Cj70K_ver. 1, genotyping showed that 53,378 out of 73,274 SNPs obtained a reliable genotype (Hiraoka et al. unpublished) for 386 plus trees; therefore, from the genotyped SNPs and 89,993 SNPs used for secondary SNP discovery, 73,638 new SNPs obtained following removal of duplicate and undesignable SNPs and addition of 2 SNPs (cj_gSNP01452 and cj_gSNP0438) reportedly near the ms-1 locus [17] were selected as Axiom_Cj70K_ver. 2 (Fig. 2, Additional file 1). The 73,638 SNPs based on isotigs covered 26,689 isotigs with an average number of SNPs per isotig of 2.76. Nearly 40% had one SNP per isotig and nearly 60% had two or more (Fig. 3).

Phenotyping

In this study, the fertile and male-sterile progeny segregated at 148:42 (approximately a 3:1 ratio; chi-square test P > 0.01), suggesting that inheritance of the male-sterile trait is controlled by a single recessive Mendelian locus.

Developing male-sterile ‘Sosyun’ pollen has aberrant morphology in the transient phase from the tetrad to microspore stage and thereafter [16]. This phenotype resembles that of sterile trees controlled by the recessive ms-1 locus, which is also manifested at the tetrad stage [4, 5, 17], suggesting that ‘Sosyun’ male sterility may be controlled by ms-1.

Genotyping

A total of 73,638 SNP markers were used for genotyping 190 progeny from the F₂ and their three contributing parents (Fig. 4). Out of 73,638 SNP markers, 7168 were polymorphic and categorized as poly-high resolution by the Axiom quality control criteria, segregating in F₂ progeny in a way concordant with the genotype of their parents and grandparents, with no missing genotypes. Out of the 7168 markers and 31 SSRs, 7199 markers in total, 3205 were used to construct a genetic linkage map after excluding markers with genotypic segregation significantly different from Mendelian segregation (when markers with identical segregation pattern were included, the number rose to 6629 markers; Fig. 2, Additional file 6).

Construction of genetic linkage map and QTL analysis

The high-density linkage map constructed from the 3205 markers allowed assignment of 11 linkage groups, covering 1492.8 cM (mean distance between adjacent markers of 0.47 cM/marker) and consisting of 4649 expected genes (Fig. 5, Additional file 6). The total genetic distance was close to the length of the reference map (1405 cM) for Japanese cedar, and the order of anchor markers (SSRs) in the new map also corresponded to that in reference map [17]. The average interval between markers in the new map was similar to that in reported reference maps for this species (1.1 cM/marker [17], 0.49 cM/marker [20]) and in map for other conifers (0.58-0.62 cM/marker in P. taeda, [36, 37], 0.6 cM/marker in a consensus map for P. taeda and Pinus elliottii [38], 0.93 cM/marker in Pinus pinaster [39], 0.92 cM/marker in a consensus map for Picea glauca and Pinus marina [40]).

We identified one notable QTL peak for male sterility. An LOD score of the strongest QTL peak was observed for the trait between 33.4 and 34.5 cM on linkage group 9 (LG9; Fig. 6), which overlapped with the reported QTL region for ms-1 [17, 18]. Eleven marker loci located to five isotigs, listed in Table 1, were found in the region. These loci were able to discriminate male-sterile progeny in the F₂ mapping population with an accuracy of 100%, higher than markers reported previously (96-98.5%) [17, 18], suggesting that these loci may be within the causative gene or tightly linked to the male sterility gene of ‘Sosyun’.

Locus (marker name)	Linkage group	Position (cM)	LOD score	Phenotypic variation explained	EST name	Sequence position	SNP	Segregation type	Description	Length of EST (bp)	E-value	Homologous Arabidopsis locus
AX-115676114	9	33.414	99.99	100.0	MF_BC_CLC_contig_46230	309	A/G	lmxll	SNARE associated Golgi protein family	530	4.29	AT4G14950.2
AX-115695236	9	33.414	99.99	100.0	reCj11611	1639	A/G	lmxll	Glycosyl hydrolase family 47 protein	2558	0.00	AT5G43710.1
AX-115678020	9	33.414	99.99	100.0	reCj19250	1927	A/G	hkxhk	PRH75 DEAD box RNA helicase (PRH75)	2841	0.00	AT5G62190.1
AX-115675027	9	33.414	99.99	100.0	reCj19250	2335	A/G	hkxhk	PRH75 DEAD box RNA helicase (PRH75)	–	0.00	AT5G62190.1
AX-115677965	9	34.123	947.48	100.0	MF_BC_CLC_contig_27478	273	A/G	nnxnp	Pre-mRNA-splicing factor 3	2139	0.00	AT1G28060.1
AX-115681104	9	34.475	99.99	100.0	reCj19314	126	T/G	hkxhk	Galactosyltransferase family protein	2081	0.00	AT5G62620.1
AX-115669347	9	34.475	99.99	100.0	reCj19314	166	T/C	nnxnp	Galactosyltransferase family protein	–	0.00	AT5G62620.1
AX-115676564	9	34.475	99.99	100.0	reCj19314	432	A/G	nnxnp	Galactosyltransferase family protein	–	0.00	AT5G62620.1
AX-115664995	9	34.475	99.99	100.0	reCj19314	859	T/G	hkxhk	Galactosyltransferase family protein	–	0.00	AT5G62620.1
AX-115679030	9	34.475	99.99	100.0	reCj19314	1623	A/G	nnxnp	Galactosyltransferase family protein	–	0.00	AT5G62620.1
AX-115665901	9	34.475	99.99	100.0	reCj19314	2587	T/G	nnxnp	Galactosyltransferase family protein	–	0.00	AT5G62620.1

Putative causative genes

In the region where an LOD score ≥ 99.9 was detected on LG9, the 11 loci derived from the five isotigs near the significant QTL explained 100% of the phenotypic variation. The A. thaliana homologs of these isotigs in the TAIR10 database are shown in Table 1. These putative causative genes for the male-sterile trait have been reported as causative genes of male sterility in other species [41–46].

As one example, pre-mRNA splicing plays a role in removing introns from nuclear genes to generate functional mRNA in most of eukaryotes, regulating important developmental mechanisms. In rice, temperature-sensitive splicing (controlled by a cis splicing cite, small nuclear RNA, trans pre-mRNA splicing protein and SR protein via a complex cell signaling pathway) is thought to be caused by photoperiod and temperature-sensitive genetic male sterility [41]. In A. thaliana, pre-mRNA-splicing factor 3 (AT1G28060.1; 33.4 cM; MF_BC_CLC_contig_27478) is a component of U4/U6 small nuclear ribonuclear protein particles, and encodes RDM16, which may indirectly regulate DNA methylation and other aspects of gene transcription [47]. Mutation of RDM16 and insertion of the corresponding DNA fragment into the At1g28060 promoter of an rdm16 ros1 mutant led to morphological defects of leaves and siliques [47]. DEAD-box proteins encode RNA-dependent ATPases or ATP-dependent RNA helicases, which mediate conformational changes with respect to spliceosome assembly and disassembly through all of the processes dealing with RNA including its synthesis, modification, cleavage and degradation [41, 48]. In male sterility of rice, abnormal programmed cell death during the degradation of tapetal cells is caused by disrupted apoptosis inhibitor5 (API5), a nuclear protein that interacts with two DEAD-box ATP-dependent RNA helicases, API5-INTERACTING PROTEIN1 (AIP1) and AIP2, which regulate the expression of the cysteine protease gene CP1 [42]. These observations support the genes related to splicing that were associated with the significant QTL as candidates for the causative gene for male sterility in our Japanese cedar population.

Genes related to cell wall components were also associated with the significant QTL. Glycosyl hydrolases (GHs) are involved in the metabolism of various carbohydrate containing compounds present in plant tissues, with a major function in metabolism of most cell wall polysaccharides [49]. GH family protein 17, which is involved in pollen wall development, is a candidate gene for rice cytoplasmic male sterility based on F₂ fine-mapping [43]. In Arabidopsis, excess microsporocytes1 (ems1)/extrasporogenous cells (exs) mutants, which produce excess microsporocytes at the expense of the tapetum, and a dysfunctional tapetum1 (dyt1) mutant show decreased expression of GH17 compared with wild type [50]. The dyt1 mutant is suggested to act downstream of ems1/exs [51]. Similarly, in Chinese cabbage, the expression of most GH family proteins is downregulated in male-sterile individuals [44]. The candidate genes in Table 1 also include a galactosyltransferase family protein gene sequence. In Arabidopsis, this protein (AT5G62620: hydroxyproline-O-galactosyltransferase) is involved in arabinogalactan protein function, associated with important roles in plant growth and development [52]. In maize, the loss of function of the causative gene Male sterile 8 (MS8), putatively encoding β-1,3-galactosyltransferase, leads to an abnormal phenotype in epidermal and tapetal cells [45], and a functional connection between arabinogalactan protein and MS8 is suggested in early anther development [46].

These observations thus indicate each of these genes as candidates for the causative gene for male sterility in Japanese cedar. We will be comparing these genes in male-sterile and male-fertile individuals.

Marker development and validation

Since conventional tree breeding has been constrained by delays for evaluation of most economically important traits, which are expressed only at the adult stage, and hence by long breeding cycles, improvement of trees using MAS for selection of markers tightly associated with QTLs of important traits is extremely attractive [53–55]. However, in coniferous species, the identified QTLs for most economically important traits such as growth and wood properties generally explain only a small portion of phenotypic variation, so a MAS approach is not likely to be a successful strategy due to the complex genetic backgrounds for these traits [56–60], with the exception of some pest resistance genes [61, 62]. Because Japanese cedar pollinosis is a serious social problem in Japan, the male-sterile trait is mandated for breeding these trees. MAS is considered an effective and attractive solution, because the male-sterile trait is controlled by a recessive major gene [17, 18, 42, 46], and hence a definitive identification is expected if tightly linked markers are available. We identified one notable region for male sterility, and around it 11 marker loci derived from 5 isotigs. The relationship of the 11 marker locus genotypes with male sterility in the F₂ population is shown in Fig. 7. Out of the 11 loci, all the male sterile individuals were homozygotes for one allele at four loci (AX-115678020, AX-115675027, AX-115664995, and AX-115681104), and the segregation ratio of the four loci was statistically consistent with Mendelian inheritance (1:2:1; chi-square test: P > 0.01). This finding indicated that these markers may also be able to discriminate heterozygotes (cryptic carriers). As male-sterile trees cannot be used as pollen donors, cryptic carriers with superior characteristics are essential for improving sterile trees and for establishing seed orchards, where sterile trees are pollinated by cryptic carriers, and sterile seedlings can be screened from the fertile cryptic carriers in the nursery using GA3 treatment to artificially promote induction of male strobili. To examine the potential of markers for precisely discriminating cryptic carriers, we developed four SNP markers (two markers for SNPs on reCj19250 (PRH75 DEAD-box RNA helicase) and two markers for SNPs on reCj19314 (galactosyltransferase family protein); Fig. 8, Table 2) that can be analyzed by a SNapShot assay, and genotyped 1082 plus trees and breeding materials selected from the Kanto breeding region. The allele and genotypic frequencies of the four markers are given in Table 3. The allele frequencies of nucleotide A in reCj19250_1927 and nucleotide T in reCj19250_2335, which are homozygous in ‘Sosyun’ were both 0.004, whereas the allele frequencies of nucleotide C in reCj19314_126 and reCj19314_859 in ‘Sosyun’ were 0.143 and 0.136, respectively (Table 3). The SNPs retained by ‘Sosyun’ at the two markers of reCj19250 can be regarded as very rare (≤0.01), and only seven plus trees were heterozygous at these two markers (Table 3). To confirm the screening potential of the two markers of reCj19250 as MAS markers outside of the mapping population, three of the seven heterozygous plus trees were crossed with ‘Sosyun’ or another heterozygous plus tree. The 63 progeny from the three families were genotyped using the SNaPshot assay and phenotyped for male sterility (Table 4, Additional file 7). Sterile progeny from all three families segregated, illustrating the effectiveness of the markers for screening cryptic carriers. The genotype arrays of the two markers for reCj19250 perfectly matched the male-sterile phenotype in all three families, whereas the markers for reCj19314 failed to discriminate sterile individuals precisely in one family. Thus, the reCj19250 markers are most suitable for MAS of male sterility.

Locus position	Forward/Reverse	Primer sequences	Amplicon size (bp)	Forward/Reverse	Extension primers	No.of nucleotides (bp)	Alleles (dyes) detected
				(Sense/Anti-sense)
reCj19250_1927	Forward	5’-CTGGAAAGGTCATGTACTCTGC-3’	190	Forward	5′-1 T-AACATATGCGTTCACTACCCTGC-3’	24	G (blue), A (green)
	Reverse	5’-TCCCATCTACAGTAAGAGACATCC-3’		(Sense)
reCj19250_2335	Forward	5’-GTTCTAGACGTGGCGGCAG-3’	180	Reverse	5’-TGCTCAATCCTGTACACAGTACCTGCAC-3’	28	C (yellow), T (Red)
	Reverse	5’-ACCTCGGAATGAACATTATGTGTTGAG-3’		(Anti-sense)
reCj19314_126	Forward	5’-CGCTCTTCGTGTTGGCTTAG-3’	133	Reverse	5′-7 T-TCTCTGTGTGCCGACGAAACCGAAGT-3’	33	A (green), C (yellow)
	Reverse	5’-GGCTCAGAGATTACCACCTCG-3’		(Anti-sense)
reCj19314_859	Forward	5’-GGTGGGCCTGTTTCCTATTCCG-3’	137	Reverse	5′-11 T-GGCCTCCATTGCAGACTTCTCAAACTG-3’	38	A (green), C (yellow)
	Reverse	5’-GCAGTTTCATCATCCATTACGGC-3’		(Anti-sense)

Locus	Genotype	Observed number	Frequency	Allele	Observed number	Frequency
reCj19250_1927	GG	1074	0.993	G	2155	0.996
	GA	7	0.006	A	9	0.004
	AA	1	0.001
reCj19250_2335	CC	1074	0.993	C	2155	0.996
	CT	7	0.006	T	9	0.004
	TT	1	0.001
reCj19314_126	AA	803	0.742	A	1854	0.857
	AC	248	0.229	C	310	0.143
	CC	31	0.029
reCj19314_859	AA	814	0.752	A	1869	0.864
	AC	241	0.223	C	295	0.136
	CC	27	0.025

Phenotype and developed marker	Number of identified type		Crossed families
			Plus-tree A (Ms1/ms1) × Plus-tree B (Ms1/ms1)	'Sosyun' (ms1/ms1) × Plus-tree B (Ms1/ms1)	'Sosyun' (ms1/ms1) × Plus-tree C (Ms1/ms1)
		Total number of Progeny	22	18	23
Phenotype	Fertile progeny		19	11	7
	Sterile progeny		3	7	16
reCj19250_1927	GG (MS1/MS1)		10	-	-
	GA (Ms1/ms1)		9	11	7
	AA (ms1/ms1)		3	7	16
reCj19250_2335	CC (MS1/MS1)		10	-	-
	CT (Ms1/ms1)		9	11	7
	TT (ms1/ms1)		3	7	16
reCj19314_126	AA (MS1/MS1)		10	-	-
	AC (Ms1/ms1)		9	11	23
	CC (ms1/ms1)		3	7	0
reCj19314_859	AA (MS1/MS1)		10	-	-
	AC (Ms1/ms1)		9	11	23
	CC (ms1/ms1)		3	7	0

This validation demonstrated that two of the markers could be used to i) screen for cryptic carriers outside of the mapping population and ii) precisely discriminate male-sterile individuals in the progeny. Thus, the markers we developed can realize the MAS approach for the first time in coniferous tree breeding.

RNA extraction and sequencing

For construction of EST libraries, we used several organs: the cambium region in the dormant season, sapwood and heartwood, apical shoots, gibberellin (GA)-treated shoots, male strobili from male-sterile individuals, male strobili from male-fertile individuals, and roots of seedlings (Additional file 2). Total RNA was isolated from these organs using an RNeasy Plant Mini kit (QIAGEN, Gaithersburg, MD, USA) and RNA Isolation Reagent (Thermo Fisher Scientific, Waltham, MA USA). The quality of total RNA was assessed via an Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). cDNA from each organ was synthesized from a mixture of RNA samples by nebulization, adaptor ligation, emulsion PCR and sequencing on a Roche 454 Genome Sequencer platform (Roche/454 Life Sciences, Branford, CT, USA) using FLX or Titanium technology, which was done at Hokkaido System Science Co., Ltd. (Sapporo, Hokkaido, Japan).

Assembly of EST sequences from several organs and construction of reference sequences of C. japonica

The ESTs sequenced by the Roche 454 Genome Sequencer were trimmed of adapter sequences and poly(A/T) sequences by the cutadapt tool [63]. Then, low-quality sequences and short sequences (< 50 bp) were removed. Japanese cedar ESTs collected from the cambium region during the active season (DRA000525) [34] and shoots during the annual season (DRA001261) [35] were added following assembly. First, all of the reads from each library were assembled using GS De Novo Assembler version 2.8 software (Roche, Indianapolis, IN, USA) with default settings. Next, all of the reads from all libraries were assembled in the same way, and these assembled isotigs were regarded as reference sequences (Additional file 1).

To classify reference sequences according to organ, we compared reference sequences with isotigs of each library. Because reads were mapped, for example using the BWA software package [64], when only portions of ESTs showed high similarity, we selected a broader comparison by the BLASTN algorithm with default parameters, using reference sequences as queries and isotigs of each library as subjects. We regarded a reference sequence as expressed in an organ if the sequence showed high similarity (> 95%) to isotigs of the library from that organ and that region covered more than 95% of the length of the subject or query. Similarly, using BLASTN, reference sequences were compared with sequences from ForestGEN, an EST database of C. japonica.

SNP discovery and SNP selection for axiom genotyping

For SNP discovery, resequencing to the reference sequences was performed for cambium of eight individuals (DRA00550 Experiment ID: DRX081263-70) and bulk samples of apical shoots of eight individuals (DRA00550 Experiment ID: DRX081272), shoots of eight individuals (DRA00550 Experiment ID: DRX081271), and male strobili of two male-sterile individuals (DRA00550 Experiment ID: DRX081274) and four male-fertile individuals (DRA00550 Experiment ID: DRX081273) using the Illumina HiSeq 2000 platform (Illumina, Branford, CT, USA) at Hokkaido System Science Co., Ltd. (Additional files 1 and 2). RNA extraction and the quality checks followed the same method used for EST library construction. Using a TruSeq RNA Sample Prep kit (Illumina), cDNA synthesis from an RNA sample from each organ, nebulization, adaptor ligation (including index tagging for individual recognition), bridge PCR and paired-end sequencing were performed on the Illumina HiSeq 2000 platform.

Reads sequenced on the Illumina HiSeq system were also trimmed of adapter sequences and poly(A/T) by cutadapt. Then, reads of each library were mapped to reference sequences by BWA [64] and SNPs were identified using SAMtools software [65] with default settings.

Because there are fewer available reference sequences for male strobili, we added more sequences from known ESTs of C. japonica and from de novo assembly for SNP detection in male strobili. The known ESTs were sequences from ForestGEN that did not show similarity to the reference sequences based on a BLASTN search. Illumina HiSeq reads from male strobili were de novo assembled by the CLC Genomics Workbench software package (QIAGEN) into more than 100,000 contigs (Mass Submission System ID: IABV01000001-01109392). About 20,000 contigs > 100 bp with average coverage > 10 reads/site were selected as references for SNP detection.

The 73,274 SNPs considered robust for genotyping were compiled into a set called Axiom_Cj70K_ver. 1 (Gene Expression Omnibus (GEO) Dataset GSE95616); we used SNPs obtained from Illumina HiSeq reads from a cambium and male strobilus library (Additional file 1). In Axiom genotyping, it is difficult to genotype SNPs adjacent to other SNPs. Hence, SNPs located within 20-50 bp of other SNPs were removed. SNPs with < 125 bp of sequence separating them from other SNPs were eliminated using a custom Perl script.

Next, we added more SNPs to Axiom_Cj70K_ver. 2 (GEO GSE95618), from Illumina HiSeq sequence reads of cambium of eight individuals, male strobili of two male-sterile individuals, male strobili of four male-fertile individuals, shoots of eight individuals, bulk samples of apical shoots of eight individuals and Roche 454 reads from GA-treated shoots (DRA00550 Experiment ID: DRX081260), sapwood and heartwood (DRA00550 Experiment ID: DRX081261) and roots (DRA00550 Experiment ID: DRX081262), which were also assembled into reference sequences, mapped to the reference sequence and used to identify SNPs (Additional file 5). As done for Axiom_Cj70K_ver. 1, SNPs located within 20-50 bp of other SNPs were removed. The surrounding sequence of the 53,378 SNPs was confirmed in Axiom_Cj70K_ver. 1 and 27,166 of these SNPs were compared to eliminate redundancy.

Genotyping

Genomic DNA was extracted from progeny, parents and grandparents of the F₂ population using a DNeasy Plant Mini Kit (QIAGEN). We used SNP markers obtained by resequencing to the reference sequence, which consisted of isotig sequences from next-generation sequencing data collected from various organs at several developmental stages and in different seasons using the Roche GS-FLX system. We also used two markers (cj_gSNP01452, cj_gSNP0438) reportedly in the vicinity of ms-1 [17]. In total, 73,638 SNPs were used for genotyping using the Affymetrix GeneTitan system as described in the manufacturer’s guide. The SNP genotyping data were divided into six categories according to their clustering performance with respect to Axiom quality control criteria: (i) polymorphic high resolution, where the SNPs were clustered at a higher resolution, including at least two minor alleles; (ii) monomorphic high resolution, where the SNPs were clustered into one group except for the presence of a minor allele in two or more samples; (iii) no minor homozygote, where the SNPs were clustered into two groups except for the presence of a minor homoallele in some samples; (iv) call rate below threshold, where the genotype call rate was < 97%; (v) off-target variant, where atypical cluster properties arose from variants in the SNP flanking region; and (vi) other, where the SNP did not pass quality control. We selected the genotyping categories in (i) polymorphic high resolution and (iii) no minor homozygotes except a monomorphic allele in both parents of the progeny.

SSR markers covering 11 linkage groups of Japanese cedar have been reported [19], and 40 SSR markers were used as anchors for the genetic linkage map in this study. The genotype of each DNA sample was scored using 31 pairs of existing microsatellite primers [9–11, 19, 66]. Multiplex PCR with three or four SSR primer pairs was performed using a Multiplex PCR Kit (QIAGEN, Hilden, Germany), with 2× QIAGEN multiplex PCR master mix, 0.25 μM each primer pair, and 40 ng genomic DNA in a total volume of 10 μl. Amplification was performed in a Veriti thermal cycler (Thermo Fisher Scientific) using an initial denaturation step at 95 °C for 15 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 1.5 min, and extension at 72 °C for 1 min, with a final extension at 60 °C for 30 min. PCR product (1 μl) was mixed with 0.2 μl GeneScan 500 LIZ size standard (Thermo Fisher Scientific) and 9.8 μl of Hi-Di formamide (Thermo Fisher Scientific) prior to electrophoresis. The length of the amplified fragments was analyzed with an ABI 3130xl sequencer (Thermo Fisher Scientific) and alleles were scored with GeneMapper v5.0 software (Thermo Fisher Scientific). In total, 73,678 markers were used for genotyping (Fig. 2).