Phenotype evaluation
The tillering ability of YD was stronger than that of BN at several developmental stages or under different growth environments. In the field, the mean tiller number of YD and BN was 34.17 (range 27–39) and 23.38 (range 15–29) at the seedling stage (Feekes’ scale 29). In the greenhouse, the mean tiller number of YD and BN seedlings were 5.7 (range 2–7) and 2.46 (range 1–3), respectively (Fig. 1).
The average tiller number of RILs in three experiments were 5.24, 3.90 and 3.48, where ranges from 2.75-8.50, 2.40-7.25 and 1.40-7.20, respectively (Table 1). The distribution of tiller number of the RILs in the three trials exhibited an asymmetric normal distribution (Fig. 2), with most lines producing a low tiller number. The pairwise comparisons of tiller numbers of the RILs among the three trials were significantly correlated and exhibited a moderate correlation coefficient, ranging from 0.222 (T1 vs. T3) to 0.365 (T2 vs. T3) (p < 0.05).
QTL mapping
A total of six QTLs on chromosomes 2B, 3B, 4A, 6D, and 7B (with two QTLs) were detected in the three trials and from the mean value, and were designated qTN-2B, qTN-3B, qTN-4A, qTN-6D, qTN-7B.1, and qTN-7B.2 (Table 2).
Among these QTLs, the positive allele of “qTN-7B.2” was from BN, whereas the other five positive alleles of QTLs were from YD. Among these QTLs, qTN-7B.1 was stably detected in T1, T2, T3, and Mean value, and explained 12.91%, 10.21%, 18.89%, and 11.87% of the phenotypic variations, respectively. The QTL peak was located between the genetic markers Xgwm333 and wsnp_BE443010B_Ta_2_1 in T1, T3, and Mean value, which corresponded to the physical region of 475646000–518919079 bp in the IWGSC RefSeq 1.0 genome assembly. In addition, qTN-7B.1 identified in T2 was a subtle shift, which was located between 2ABD7ABD_wsnp_BE443010B_Ta_2_1 and 7ABD_wsnp_be518436B_Ta_2_1 and corresponded to the physical region of 518919079–525133262 bp. The other five QTLs (qTN-3B and qTN-7B.2, qTN-6B, and qTN-2B and qTN-4A) were detected in T1, T2, and T3, respectively.
Enriched potential pathways associated with tillering ability
The DEG analysis identified a set of 2806 and 2853 up-regulated genes for BN and YD, respectively. To evaluate the potential functions of genes associated with tillering ability, gene ontology (GO) enrichment analyses were conducted. A total of 431, 2469, and 1068 GO terms were enriched among 3514, 3537, and 3161 DEGs for the cellular component, molecular function, and biological process categories, respectively. The top five enriched categories were nutrient reservoir activity, manganese ion binding, ubiquitin-protein ligase binding, ubiquitin-like protein ligase binding, and nuclear chromatin (Fig. S1).
The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to identify the metabolic pathways that may be fully or partially represented by annotated DEGs. A total of 1856 DEGs were assigned to 137 KEGG pathways (Fig. S2). Among these pathways, metabolic pathways, biosynthesis of secondary metabolites, and plant–pathogen interaction were the three most highly enriched KEGG pathways.
Among metabolism pathways, the phenylpropanoid biosynthesis pathway was the most highly enriched among the DEGs, followed by starch and sucrose metabolism, and carbon metabolism.
Sequence variations of TaMoc1-7B between parents
According to annotation of CS RefSeq v2.0, the TaMoc1-7B were predicted to have only one exon. A total of 4055bp and 4053bp genomic sequences were obtained from BN and YD, respectively. The length of coding region in both two sequences were 1290bp, encoding 429 amino acids. The lengths of the 5’upstream and 3’downstream sequences of cloned TaMoc-7B in YD were 2052-bp and 711-bp, where 2054-bp and 711-bp were obtained in BN.
Alignment of sequences between two parents, 21 variations were found, 4 SNPs were at 3’ downstream region, 14 SNPs and two indels were at 5’ downstream region, and 1 SNPs were at the coding region (Fig S3). In particular, the one SNP at coding region could not affect amino acid change (CTC→CTG).
DEGs in qTN-7B.1
The confidence interval cross 133Mb (between 435789094-bp to 566605694-bp) of qTN-7B.1 was used to search candidate genes. A differential expression analysis of the RNA-sequencing data from the tiller bud for both parents was performed (Table S2). We first evaluated whether potential expression level differences of TaMoc1-7B between parents affect gene function. As the result, extremely low reads of TaMoc1-7B transcripts were detected in both parents (0, 3, 0 and 4, 5, 8 in BN and YD, respectively). The average FPKM value of TaMoc1-7B in BN and YD were 0.01 and 0.08. Next, thirty-eight showed significantly differential expression between BN and YD. Among these DEGs, 9 and 29 DEGs were up-regulated in YD and BN, respectively.
Next, sequence variation among the DEGs was surveyed to identify the most likely candidate gene. Moreover, 305 nonsynonymous SNV were identified in qTN-7B.1 region, where 44 were distributed on 20 DEGs (Table S3).
Development of genetic marker for TraesCS7B02G282100
TraesCS7B02G282100 encoding cinnamyl-alcohol dehydrogenase, which was related to Lignin biosynthesis. Based on RNA-seq data, the transcript level of BN was 3.179-fold up-regulated to YD. Besides, one SNP (T406G:p.F136V) was found at coding region that changed amino acid, where this SNP was located at conserved domain of protein (Fig 3). Searching conserved domain revealed that this SNP was located on the PLN02586 super family domain (probable cinnamyl alcohol dehydrogenase, from 13-bp to 1080-bp of CDS, searched in https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Based on this SNP, a KASP marker designate at KASP_282100 was developed. Through anchoring the genotypes of KASP_282100 in all RILs, a new genetic frame map was constructed. The KASP_282100 was mapped between Xgwm333 and 2ABD7ABD_wsnp_BE443010B_Ta_2_2, which was the qTN-7B.1 region (Table S4).