Screening of Candidate Leaf Morphology Genes by Integration of QTL Mapping and RNA Sequencing Technologies in Oilseed Rape (Brassica napus L.)

Leaf size and shape play important roles in agronomic traits, such as yield, quality and stress responses. Wide variations in leaf morphological traits exist in cultivated varieties of many plant species. By now, the genetics of leaf shape and size have not been characterized in Brassica napus. In this study, a population of 172 recombinant inbred lines (RILs) was used for quantitative trait locus (QTL) analysis of leaf morphology traits. Furthermore, fresh young leaves of extreme lines with more leaf lobes (referred to as ‘A’) and extreme lines with fewer lobes (referred to as ‘B’) selected from the RIL population and leaves of dissected lines (referred to as ‘P’) were used for transcriptional analysis. A total of 31 QTLs for the leaf morphological traits tested in this study were identified on 12 chromosomes, explaining 5.32–39.34% of the phenotypic variation. There were 8, 6, 2, 5, 8, and 2 QTLs for PL (petiole length), PN (lobe number), LW (lamina width), LL (Lamina length), LL/LTL (the lamina size ratio) and LTL (leaf total length), respectively. In addition, 74, 1,166 and 1,272 differentially expressed genes (DEGs) were identified in ‘A vs B’, ‘A vs P’ and ‘B vs P’ comparisons, respectively. The Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to predict the functions of these DEGs. Gene regulators of leaf shape and size, such as ASYMMETRIC LEAVES 2, gibberellin 20-oxidase 3, genes encoding gibberellin-regulated family protein, genes encoding growth-regulating factor and KNOTTED1-like homeobox were also detected in DEGs. After integrating the QTL mapping and RNA sequencing data, 33 genes, including a gene encoding auxin-responsive GH3 family protein and a gene encoding sphere organelles protein-related gene, were selected as candidates that may control leaf shape. Our findings should be valuable for studies of the genetic control of leaf morphological trait regulation in B. napus.


Introduction
Oilseed rape (Brassica napus L.) is an important oil crop worldwide that provides both edible oil and industrial materials. Many important traits are controlled by quantitative trait loci (QTLs), such as oil content [1], seed weight [2], flowering time [3], silique length [4] and yield [5]. Meanwhile, some of the candidate genes for these traits have been identified and their functions characterized. Leaf morphology, which has a significant impact on yield, is an important agronomic trait, but information on gene levels in oilseed rape is lacking. Leaves, as lateral organs and determinants of growth, develop from flanking regions of the shoot apical meristem (SAM) [6]. Their morphology involves a balance between cell proliferation and polar cell division and expansion along the proximal-distal, medial-lateral and adaxial-abaxial axes [7][8][9][10]. Moreover, plant hormones, such as strigolactone, auxin, cytokinin, and gibberellin (GA), influence leaf development and morphogenesis [11][12][13].
QTL-based approaches are frequently used to reveal loci that regulate natural phenotypic diversity. QTLs have been identified for leaf morphology in Brassica rapa [14,15], tomato [16], grape [17], oak [18], maize [19], and Arabidopsis [20], suggesting that leaf morphological variation is controlled by multiple genes. QTLs for leaf architecture have also been identified in Brassica species [21][22][23][24]. Lan and Paterson (2001) [25] identified QTLs, that explained 45% of the phenotypic variation of lamina length in Brassica oleracea, and three co-localized with QTLs affecting leaf width (LW). Ten QTLs were later identified that influenced leaf traits using three different mapping populations [21]. Lobe depth and leaf hairiness are controlled by Brassica rapa gibberellin 20-oxidase 3 (BrGA20OX3) and Brassica rapa glabra1, which were identified in a synteny analysis of the Brassica rapa and Arabidopsis thaliana genomes [26,27]. Leaf hairiness and seed coat color are controlled by the Brassica rapa transparent testa glabra1 gene, which was successfully cloned [28].
Until now, information on the genetics of leaf morphology in B. napus has been limited. Here, we used QTL mapping and RNA sequencing (RNA-Seq) to uncover the genetic architecture of leaf morphology in B. napus. In total, 31 QTLs for leaf morphology traits were identified on 12 chromosomes and differentially expressed genes were identified. These findings will be very useful for studies of the genetic control of leaf morphological trait regulation in B. napus.

Plant material, growth conditions and phenotyping
The genetic map constructed by Liu et al. (2013) [29] using an F 9 RIL mapping population derived from a cross between the parental lines GH06 (female, fewer lobes) and P174 (male, more lobes) was used for QTL mapping in our study. For phenotyping, the parental lines and 172 F 11 plants were grown in 3 replications from October 2015 to May 2016 in the breeding nursery of Southwest University in Beibei (29˚39´N, 106˚36´E), Chongqing (China).
The leaf characteristics were scored on fully developed leaves 60 days after sowing. Lamina length (LL), lamina width (LW), petiole length (PL), leaf total length (LTL), lobe number (PN) and the lamina size ratio (calculated as LL/LTL) were scored and are shown in Fig 1(A). All of the above traits were determined for five plants per line and three leaves per plant. The values of LL, LW, PL and LTL were measured using a ruler as described by Wang et al. (2015) [15]. The SAS 9.0 program (SAS Institute, Inc., Cary, NC, USA) was used for distribution analysis.

QTL analysis
The linkage map was constructed using 2,795 single-nucleotide polymorphism (SNP) markers with a total map length of 1,832.9 cM and an average distance of 0.66 cM. WinQTL cartographer version 2.5 (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm) was used to screen QTLs for leaf morphological traits by composite interval mapping (CIM). Additionally, 1000 permutations were analyzed to calculate the LOD threshold and the corresponding LOD threshold was used to identify significant QTLs for the corresponding traits. QTL positions were drawn with the software Mapchart [30]. The naming rules for QTLs described by McCouch et al. [31] were used with minor modifications. A QTL designation begins with the QTL and trait abbreviation, such as ''qLL" (q, QTL; LL, lamina length), followed by the linkage group, and finally, the serial number of the QTL in the linkage group (e.g., qLL-A2-1).

RNA isolation and transcriptome sequencing
The fresh leaves of five lines with more lobes (A) and five lines with fewer lobes (B) were selected from the RIL population, along with leaves from the dissected lines (P), as illustrated in Fig 1(B). These samples were immediately placed in liquid nitrogen and stored at -80˚C. Total RNA was extracted using the Plant RNA Mini Kit (TIANGEN, Inc., China). The Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) was used to determine the total RNA. qualities. cDNA libraries were constructed and RNA-Seq was performed on an Illumina HiSeq 2000 platform. Paired-end reads of 100-bp were generated for A, B and P samples by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China) according to the manufacturer's instructions. In addition, sequencing data were uploaded to the NIH Short Read Archive under accession number SRP079113.

RNA sequencing data analysis
After removing the adapter sequences and low quality sequences from the raw data, clean reads were mapped to the B. napus reference genome (http://www.genoscope.cns.fr/brassicanapus/) and were then assembled using TopHat 2.0.0 and Cufflinks [32]. The FPKM (fragments per kilobase of exon per million mapped fragments) values were used to estimate the gene expression levels, and differentially expressed genes (DEGs) between two samples were identified with Cuffdiff using the criteria FDR<0.01 and |log 2 (fold change)|>1. To identify possible transcription factors (TFs) and plant hormone genes, DEGs were aligned to known TFs and plant hormone genes in the A. thaliana genome that were downloaded from The Plant Transcription Factor Database (http://planttfdb.cbi.pku.edu.cn/index.php) [33] and Arabidopsis Hormone Database 2.0 (http://ahd.cbi.pku.edu.cn/) [34], respectively.

qRT-PCR analysis
To confirm the accuracy of the sequencing data, the expression patterns of 25 randomly selected genes were determined using qRT-PCR analysis. For all three samples, 1 μg of RNA was treated with Turbo DNA-free (Ambion, Austin, TX) and was reverse transcribed into cDNA according to the manufacturer's instructions (Bio-Rad). Gene-specific primers were designed with Primer3 software (http://frodo.wi.mit.edu/primer3/) and are listed in S1 Table. The 20-μL PCR reactions, each containing 10 μl of SYBR Green Supermix (Bio-Rad), 2.0 μl of cDNA, 0.4 μM each primer, and variable amounts of distilled water, were assembled, and the following cycling program was used for qRT-PCR: 95˚C for 30 s; 35 cycles of 95˚C for 5 s, followed by annealing at 56-67˚C (depending on the primers used) for 30 s. Actin was used as a reference gene to normalize template amounts. The 2 -ΔΔCT method was used to estimate the fold change according to Livak & Schmittgen [36]. The data were analyzed using Bio-Rad CFX Manager software. Three biological replicates with three technical replicates were performed for each reaction.

Analysis of leaf morphological traits
The frequency distributions of leaf morphological traits in the RIL population are summarized in Table 1 and Fig 2. We examined a wide range of leaf morphological traits in this study, which showed the transgressive and continuous distributions expected for characters displaying quantitative trait segregation. This observation suggested that there were polygenic effects on leaf morphology and that the data were suitable for QTL analysis (Table 1 and Fig 2).

Detection of QTLs that influence leaf morphology
We analyzed leaf morphology traits using a CIM approach. A total of 31 QTLs for leaf morphology traits were identified on 12 chromosomes. These QTLs explained 5.32-39.34% of the phenotypic variation (PV) ( Table 2 and Fig 3). There were 8, 6, 2, 5, 8, and 2 QTLs for PL (petiole length), PN (lobe number), LW (lamina width), LL (lamina length), LL/LTL (the lamina size ratio) and LTL (leaf total length), respectively ( Table 2). PL: Eight QTLs that were located on five chromosomes explained 86.08% of the PV, and the alleles of six of the QTLs came from P174. qPL-A01-1, a major QTL with the largest effect, was detected on chromosome A01 (Fig 3), which accounted for approximately 24.94% of the PV. An allele from GH06 increased PL by 2.49% (additive effect, Table 2). Four QTLs located on chromosome C01 explained 37.81% of the PV, and one QTL located on each of the chromosomes A03, A10 and C08 explained 6.02%-8.87% of the PV (Table 2).
LL: Five QTLs located on four chromosomes were detected for LL that explained 41.36% of the PV with 5.58%-10.49% of the PV accounted for by a single QTL. The additive effect of qLL-A02-1 came from P174 (-0.67), and the remaining four QTLs came from the GH06 parental line ( Table 2).
LTL: Two QTLs for LTL on chromosome C01 were detected; they explained 5.60% and 10.67% of the PV, respectively. Both of them decreased LTL due to a negative additive effect (Table 2). Moreover, the confidence intervals of qLTL-C01-2 and qPL-C01-3 overlapped, and the additive effects were both negative, indicating that the loci contributed to both PL and LTL.
LL/LTL: Eight QTLs on five chromosomes were detected, explaining 90.56% of the PV with 5.32%-39.34% of the PV accounted for by a single QTL. One major QTL, qLL/LTL-A01-1 on chromosome A01, explained 39.34% of the PV, and the alleles came from P174, as indicated by an additive effect of -0.04 (Table 2). Interestingly, the confidence intervals of qPN-A01-2 and qLL/LTL-A01-1 overlapped with a reversed additive effect, suggesting that these loci had pleiotropic effects (Table 2).

Illumina sequencing and global analysis of gene expression
To obtain a global view of the leaf transcriptome in the different lines, high-throughput RNA--Seq was performed on the three samples (A, B and P) (Fig 1B). In total, 147,360,468 raw reads were generated from the three samples, and more than 96% of the reads remained after low quality and adapter reads were removed, which was sufficient for a quantitative analysis of gene expression. The clean reads were mapped to the reference genome database using TopHat software. Of the total reads, 74.55%-75.16% were mapped to either multiple (3.59%-4.49%) or unique (70.06%-71.57%) genomic locations, and the remaining 24.84%-25.45% were unmapped (Table 3).
In total, 67,554 expressed genes in A, B and P libraries were detected. Additionally, 25.6%-26.9% of genes had FPKM values lower than 1.0, 18.0%-18.8% of genes had FPKM values between 1.0 and 3.0, 34.7%-35.8% of genes had FPKM values between 3.0 and 15.0, 14.0%-14.9% of genes had FPKM values between 15.0 and 60.0, and 5.5%-5.6% of genes had FPKM values more than 60.0 in the three libraries (Fig 4A). A Venn diagram was used to show the distributions of expressed genes in samples A, B and P (Fig 4B). Among these genes, 54,778 (81.1%) were expressed in all three samples; 3,778, 1,449 and 1,445 were co-expressed in A and B, B and P, and A and P, respectively; 1,868, 1,983 and 2,253 genes were unique to A, B and P, respectively.

Transcriptome changes in the three different types of leaves
To identify the candidate genes responsible for leaf development, we identified the genes that were differentially expressed between each of the samples using the criteria: |log 2 fold change|> 1.0 and FDR<0.01. Between the A and B libraries, 74 DEGs were detected, with 37 up-regulated and 37 down-regulated genes (Fig 4C). In total 1,166 DEGs were detected in the comparison of the A and P libraries, with 485 up-regulated and 681 down-regulated genes, and 1,272 DEGs were identified between the B and P libraries, with 600 up-regulated and 672 down-regulated genes ( Fig 4C). Furthermore, 24 DEGs were identified in the "A vs B", "A vs P" and "B vs P" comparisons, and 53, 38 and 897 were co-expressed in the "A vs B" and "A vs P", "A vs B" and "B vs P", "A vs P" and "B vs P" comparisons, respectively. Seven, 240 and 361 genes were specifically expressed in the "A vs B", "B vs P" and "A vs P" comparisons, respectively (Fig 4D). Moreover, the fold changes (up-or down-regulated) of most DEGs in the "A vs B", "B vs P" and "A vs P" comparisons were approximately 2-8 ( Fig 4E).

Functional classification of DEGs involved in the development of leaf shape
To determine the important roles of the DEGs described above in the development of leaf shape, GO and KEGG analyses were used to classify their functions. In total, 35, 44 and 44 functional groups were categorized for DEGs from the A vs B, A vs P and B vs P comparisons, respectively (S1 Fig). In each of the three main categories, the terms cell (GO: 0005623), binding (GO: 0005488), and cellular process (GO: 0009987) were dominant. High percentages of the genes were also associated with the terms cell part (GO: 0044464), organelle (GO: 0043226), catalytic activity (GO: 0003824) and metabolic process (GO: 0008152) (S1 Fig). In addition, 4, 5 and 5 KEGG pathways were identified in the A vs B, A vs P and B vs P comparisons, respectively (S2 Table). The top three most enriched pathways assigned to DEGs were Ribosome (ath03010), Carbon metabolism (ath01200) and Carbon fixation in photosynthetic organisms (ath00710).

Transcription factors and plant hormones involved in leaf morphology
TFs play key roles in plant organ development [37]. In our study, 31 DEGs belonging to 21 TF families were identified in the DEGs and were shared between the A vs P and B vs P comparisons ( Fig 5A).  [11,12]. In this study, 63 hormone-associated DEGs were identified in the DEGs shared between A vs P and B vs P comparisons ( Fig  5B). The top three hormones with the most associated DEGs were ABA (18), ET (15) and auxin (13), accounting for 73% of the hormone-associated DEGs. Almost half of the total hormone-associated DEGs were up-regulated, and all of the GA and BR genes were up-regulated.

Expression analysis of the orthologous genes influencing leaf morphology in Arabidopsis
It was unknown whether any of the DEGs were associated with controlling leaf morphology in B. napus. In total, 201 genes with roles in leaf development were identified in the B. napus genome based on their corresponding homologs in the A. thaliana genome. To identify pivotal genes that controlled leaf morphology, RNA-Seq data were filtered using the criteria: |log 2 fold change|!1.0 and FPKM !1 comparing P lines (for at least one sample). As illustrated in Table 4, 9 GRFs (1 GRF1, 2 GRF3, 1 GRF4, 3 GRF5 and 2 GRF8), 7 GA-related genes (2 GA20OX3 and 5 gibberellin-regulated family proteins), 2 AS2s and 2 KNATs (1 KNAT4 and 1 KNAT5) were identified. Moreover, 13 of these genes were up-regulated, including 8 GRFs, 2 AS2, 2 KNATs and 1 GA-regulated family protein.
Screening for candidate genes that control leaf morphology by the integration of QTL mapping and RNA sequencing In total, 1,205 genes were identified in QTL regions by aligning SNP marker physical locations to the oilseed reference genome. The expression profiles of these genes in the three samples were determined via RNA-Seq (S4 Table). Thirty-three candidate genes (21 up-regulated and 12 down-regulated) were identified using the criteria: |log 2 fold change|!1.0 and FPKM !1 comparing B samples (for at least one sample). These genes included two encoding auxinresponsive GH3 family proteins (BnaA06g06420D and BnaA06g06400D), the sphere organelles protein-related gene (BnaA06g07690D), the ARF-GAP domain 4 (BnaA06g06680D), and two clathrin-encoding genes (BnaA01g31000D and BnaA01g31010D) ( Table 5).

Verification of transcriptome sequencing
To confirm the results of RNA sequencing, 25 DEGs were selected randomly for qRT-PCR analysis (Fig 6). All DEGs exhibited similar expression patterns supporting the reliability of the data from RNA sequencing. Table 4. Expression analysis of the orthologous genes that influence leaf morphology in Arabidopsis. A: Extreme lines with more lobes; B: extreme lines with fewer lobes; P: dissected lines.

Discussion
In our study, QTL and RNA-Seq analyses were used to screen for candidate genes that could influence leaf morphology in B. napus. We identified 31 QTLs for leaf morphological traits that were located on 12 chromosomes. One major QTL, located on chromosome A01, controlled PL and LL/LTL (Table 2). Until now, no QTL mapping studies of leaf morphological traits in B. napus have been reported. Leaf morphology shows extensive variation among different plant species and even among varieties of the same species and is strongly influenced by genetics. Experiments in model species have identified some of the genes and mechanisms controlling leaf size and shape [7,8]. Screening of our QTL and RNA-Seq results, we identified 33 candidate genes, and the expression profiles of 25 genes were verified by qRT-PCR. Unfortunately, no evidence that these genes regulated leaf size and shape has been previously reported. Further studies, especially those focusing on BnaA06g07690D (sphere organelles protein-related), BnaA06g06420D (Auxin-responsive GH3 family protein) and BnaA06g06680D (ARF-GAP domain 4), should be conducted to confirm the roles of these genes in leaf development.
To understand the molecular mechanisms controlling the final size and shape of B. napus leaves, we used RNA-Seq to measure the expression levels of the B. napus homologs of these potential regulators in three samples with different leaf shapes. There were 2, 2, 5, 9, 1 and 1 homologs for AS2, GA20OX3, gibberellin-regulated family protein, growth-regulating factor, KNOTTED1-like homeobox gene 4 (KNAT4), and KNOTTED1-like homeobox gene 5 (KNAT5), respectively (Table 4). Coordinated cell proliferation and cell expansion are the main types of growth affecting leaf size and shape [38]. In model systems such as A. thaliana, several genes that control leaf growth have been identified. These genes determine leaf size and shape by regulating the balance between cell proliferation and cell expansion.
In our study, both AS2 and KNOX1 genes were down-regulated in the P sample compared with their expression in the A or B samples. BnaCnng20070D, which is homologous to AT5G11060 (KNAT4), and BnaA01g04870D, which is homologous to AT4G32040 (KNAT5), are KNOTTED1-like homeobox genes with important roles in promoting leaf complexity [39]. KNOX1 genes are expressed in the SAM and maintain indeterminacy of leaf stem cells [40]. KNOX1 genes also play pivotal roles in leaf initiation in most plants with compound leaves. Their influence on leaf development has been confirmed by altering their expression profiles [41][42][43][44]. AS2, a key negative regulator of KNOX genes, is expressed in leaves [45]. Defects in leaf patterning have been detected in as2, whereas in AS2-overexpressing lines, lamina outgrowths occurred on the abaxial side of the leaf. There are 36 homologs of the nine AtGRFs in B. napus. In this study, nine GRFs, including three GRF5s (homologous genes of AT3G13960), were identified. BnaC08g28950D was up-regulated while BnaA05g25480D and BnaA03g33180D were down-regulated in the P sample compared with the A or B samples. Previous studies have indicated that AtGRF5 plays important roles in leaf cell proliferation [46]. Another gene, ANGUSTIFOLIA3 (AN3), works together with AtGRF5 and is involved in leaf cell proliferation, leaf size and shape regulation [47]. The double mutants of an3 and atgrf5 exhibited narrow-leaf phenotypes because of reduced cell numbers. By contrast, cell number increased and leaves grew larger than the wild type in AN3or AtGRF5-overexpressing lines. These results indicated that functional differentiation occurred in the homologous GRF5 genes.
Moreover, plant hormones, such as GAs, interact with KNOX pathways, influencing leaf shape and final size [13]. In this study, BnaC09g48370D and BnaA10g23640D, homologs of AT5G07200 (GIBBERELLIN 20 OXIDASE 3), were identified and were up-regulated in P samples compared with A or B samples. Furthermore, five gibberellin-regulated family proteins were also identified, four of which were up-regulated in P samples compared with A or B samples. The GA-biosynthesis gene GA 20-oxidase [48] is suppressed by the Nicotiana tabacum Homeobox15 (NTH15) gene, and GA levels subsequently decreased in tobacco leaves of NTH15-overexpressing lines [49]. Previous studies have also shown that GA synthesis is regulated by KNOX genes in Arabidopsis. The altered leaf shape phenotype is partially alleviated by increasing GA in Arabidopsis plants through the ectopic expression of KNOX genes in the leaves [50]. Li et al., [26] also identified a QTL for leaf lobe and leaf hairiness containing GIBBERELLIN 20 OXIDASE 3.

Conclusions
In summary, a total of 31 QTLs for the leaf morphological traits evaluated in this study were identified on 12 chromosomes, explaining 5.32-39.34% of the phenotypic variation using a population of 172 RILs. In addition, 74, 1,166 and 1,272 DEGs were identified using RNA-Seq technologies in the 'A vs B', 'A vs P' and 'B vs P' comparisons, respectively. In addition, 1,205 genes were identified in QTL regions by aligning the SNP marker physical locations with the oilseed reference genome. Thirty-three candidate genes were identified by QTL mapping and RNA-Seq technologies. Our findings should be valuable for studies of the genetic control of the regulation of leaf morphological traits in B. napus.