Population design
The CUBIC design is a modified Multi-parent Advanced Generation Inter-Cross (MAGIC) population consisting of 1404 progeny. It was derived from 24 elite Chinese maize inbred lines from 4 divergent heterotic groups (Additional file 1), which included 4 founders from LvDaHongGu germplasm (旅大红骨种质): LV28 (旅28), E28, DAN340 (丹340), and F349; 4 founders from ZI330 germplasm (自330亚群种质): ZI330 (自330), ZONG3 (综3), ZONG31 (综31), and HUANGC (黄C); 15 founders from SiPingTou germplasm (四平头种质): HZS (黄早四), HYS (黄野四), TY4 (天涯4), YUANGFH (原辐黄), CHANG7-2 (昌7-2), K12, XI502 (西502), LX9801, H21, SHUANG741 (双741), Q1261, JI853 (吉853), JI53 (冀53), 5237, and 81515; and 1 founder from Yugoslavia-improved germplasm (南斯拉夫群体选系): NX110 (农系110). These 24 founders were crossed under a Complete Diallel Cross type IV (CDC IV) mating design, omitting parents and reciprocal crosses, in the summer of 2004. Thirty F1s with acceptable agronomic performance (early flowering, short ear height, and big ear size) were selected to cross further in the CDC IV design, and an additional 110 F1s were randomly selected to open pollinate in isolation in the winter of 2004. A total of 200 and 400 ears resulting from the best and the random crosses, respectively, were harvested. Seeds were mixed together in a 2:1 ratio (biased to improve population performance), and individuals were planted under open pollination from the summer of 2006 to the winter of 2008 in isolated regions for 6 generations. About 2000 ears of the most phenotypically diverse individuals were retained, and seeds were mixed in equal proportions in each generation. The population was then self-pollinated by single seed descent beginning in the summer of 2009 for another 6 generations. During the open- and self-pollinated processes, the greatest possible variation for several traits was constantly maintained to maximize the phenotypic diversity. A total of 1664 inbred lines were obtained, of which 1404 were successfully sampled and sequenced, and thus used in further analysis.
Phenotyping
All 1404 inbreds were planted in the year of 2014 at 5 locations (Additional file 3: Figure S2; N 43° 42′, E 125° 18′, Yushu City, Jilin Province; N 42° 03′, E 123° 33′, Shenyang City, Liaoning Province; N 40° 10′, E 116° 21′, Changping District, Beijing City; N 38° 39′, E 115° 51′, Baoding City, Hebei Province; N 35° 27′, E 114° 01′, Xinxiang City, Henan Province) in Northern China, where the 24 elite founders that served as the parents of the population are the most adapted. The lines were planted with a random order at each environment to eliminate the potential confounding effect. About 17 individual plants were planted for each line, and the line Chang7-2 was planted after every 50th entry, whose phenotypes are used to correct for trend and spatial heterogeneity. Twenty-three phenotypes were measured at each location, including 6 flowering time traits: days to tasseling (DTT, measured as the interval from sowing to the day the tassel appeared in half of the individuals per line), days to anthesis (DTA, measured as the interval from sowing to the day of pollen shed for half of the individuals), days to silking (DTS, measured as the interval from sowing to the day silks emerged for half the individuals), and the intervals between them: interval between anthesis and tasseling (ATI), interval between silking and tasseling (STI), interval between silking and anthesis (SAI); 8 developmental agronomic traits: plant height (PH, vertical height from the ground to the top of the tassel, with an accuracy of 0.5 cm), ear height (EH, vertical height from the ground to the node where the top ear arises), ear leaf length (ELL, straight length of the first ear leaf), ear leaf width (ELW, width of the ear leaf measured at the widest point), leaf number above ear (LNAE, leaf or node number from ear leaf to the top; counted including the ear leaf), leaf number below ear (LNBE, leaf or node number below the ear leaf; not including the ear leaf), tassel branch number (TBN, only the primary tassels were considered), tassel length (TL, straight length of the main branch); and 9 ear traits: ear weight (EW), ear diameter (ED), ear length (EL), ear row number (ERN), kernel number per row (KNPR), kernel number per ear (KNPE), kernel weight per ear (KWPE), cob weight (CW), and length of the barren tip (LBT). For flowering times, all individuals in a line were considered, while only 5 consecutive plants in the middle of each line were measured and the average value was used for developmental agronomic traits. For measuring ear traits, 5 approximately equidimensional ears were randomly selected to represent each line.
For each trait, the variance of genotype and environment on phenotype were estimated, and the line mean-based broad-sense heritability for each trait was calculated as: \( {H}^2={\delta}_g^2/\left({\delta}_g^2+{\delta}_e^2/n\right) \), where \( {\delta}_g^2 \) is the genetic variance, \( {\delta}_e^2 \) is the residual variance, and n is the number of environments. The estimates of \( {\delta}_g^2 \) and \( {\delta}_e^2 \) were obtained by the mixed linear model, treating genotype and environment as random effects. The best linear unbiased predictor (BLUP) value for each inbred line was calculated across all environments using the mixed linear model in the R package “lme4” [38]. The BLUP values were used in subsequent analyses, including basic phenotypic statistics, correlation analysis, and GWAS. The phenotypic variance or heritability explained by either sQTL or hQTL for a given trait was estimated in an analogous way by using linear regression, and the total variances explained by integrated methods were estimated jointly.
Whole-genome sequencing (WGS)
Sample collection and DNA extraction
Young leaves in the vegetative growth stage were collected about 5 weeks after planting and flash-frozen in liquid nitrogen. Genomic DNA from each sample was extracted with the cetyltrimethylammonium bromide (CTAB) method [39] and dissolved in double-distilled water. The DNA was checked for quality and quantity on agarose gels and a Qubit fluorometer (Invitrogen). Samples having at least 1 μg of total DNA, A260/A230 above 2.0, and A260/A280 above 1.8 were selected to construct the sequencing library.
Library construction
High-quality genomic DNA was column-purified to remove protein, carbohydrates, salt, and other impurities. Purified genomic DNA (1 μg) was fragmented using a Biorupter UCD-200 (Diagenode). The fragmented DNA was inspected for quality via agarose gel electrophoresis, and fragments in the range of 300–500 bp were concentrated. DNA was repaired by adding adenine nucleotides (A) and a sequencing adapter to remove gaps at the 3′ end of the fragmented DNA; repaired DNA was then further purified by 2% agarose gel electrophoresis. DNA fragments of around 500 bp were recovered by PCR amplification and purification to obtain the final library, which was inspected again for quality by agarose gel electrophoresis and qPCR with the StepOne Plus Real-Time PCR system (ABI).
Sequencing
DNA libraries were sequenced with the Illumina HiSeq 2500 platform using V4 reagents with 125-bp paired-end reads, generating a total of 5 Tb of sequenced base pairs. The library preparation and sequencing work described above was accomplished by the BerryGenomics Company (Beijing, China).
Variant discovery
Quality control of raw reads
Adaptor sequences were removed as the first step of the analysis, and low-quality reads were trimmed by Trimmomatic (Version 0.33) [40], using the following parameters: TRAILING = 3, MINLEN = 50, and HEADCROP = 10.
Read alignment
Clean reads were mapped to the maize B73 reference genome (v3.25, downloaded from http://plants.ensembl.org) using the BWA aligner (Version 0.7.12) [41], which displays a good balance between running time, memory usage, and accuracy. The FM-index for the reference genome was produced first, and reads for each sample were aligned against the reference using the “aln/sampe” option with default parameters. The unique mapped reads extracted from SAM files were sorted and indexed using Picard (Version 1.119) [42].
Insertion/deletion realignment and base quality recalibration
The RealignerTargetCreator and IndelRealigner modules from the Genome Analysis Toolkit (GATK, Version 3.5) [43] were used to perform local realignment around InDels to correct mapping artifacts. SAMtools (version 0.1.19) [44] and UnifiedGenotyper from GATK were first used to generate a high-quality SNP set as known sites to build the covariance model and estimate empirical base qualities for each individual. Next, BaseRecalibrator and PrintReads tools from GATK were used to recalibrate base quality scores in order to correct sequencing errors and other potential experimental artifacts.
Variant calling (SNPs and short insertions and deletions, InDels)
The recalibrated BAM files were processed with SAMtools and UnifiedGenotyper from GATK. For GATK, the parameter “-glm” was set as “BOTH” to obtain SNPs and InDels simultaneously. Variant calls from the SAMtools mpileup package were identified using default parameters. Variants identified by both programs were only kept if they satisfied mapping quality (MQ ≥ 20.0) and sequencing coverage (DP ≥ 2 and DP ≤ 100). The filtered variants of each sample were further combined by GATK CombineVariants, and another re-calling process of every individual using GATK UnifiedGenotyper was performed to obtain initially integrated genotypes at the population level. Next, variants with a missing rate greater than 98% were removed, and heterozygous sites with a “Lowqual” label were excluded. The resulting variant calls were stored in a variant call format (VCF).
Variant calling pipeline evaluation
To estimate the error rate of the variant calling pipeline for the present study (i.e,. for the maize genome), a number (10×) of pair-end sequencing reads were simulated from the reference genome with 0 miss-match and the same pipeline was utilized to call variants. Theoretically, no differences should be identified; such variants could only be derived from incorrect mapping and calling, which may occur due to intrinsic duplication of genome fragments. Only 12 heterozygous SNPs and 0 InDels were called using the simulated sequences. This indicated that our maize variant calling process applied above is reliable, and systematic errors (caused by the complex genome structure) are limited and most likely to happen for heterozygous loci. With the expectation that this kind of error would occur most frequently at specific sites, we further simulated the same procedure 10 times with the aim of covering the majority of this type of error. The variant sites identified from this analysis were removed, if present, from the variant set called from real data.
Genotyping by 200K Array
A subset of 194 lines, including 24 founders and 170 randomly selected progenies, were further genotyped using the Affymetrix Axiom Maize Genotyping Array 200K (from Peking University, China). The Affymetrix Power Tools (APT) and SNPolisher package [45] were used to analyze the genotyping data in the “cel” format. Briefly, the DQC value and call rate for each sample were set larger than or equal to 0.82 and 0.97, respectively, and acceptable poly high resolution (PHR) values were retained leading to a dataset of 53,831 high confidence SNPs. The average concordance rate between the 200K Array and the whole-genome sequencing was 97%. All SNPs obtained from re-sequencing and the 200K array were integrated to perform genotype imputation for all lines, and those inconsistent genotypes were considered as unknown.
Genotype imputation
Missing genotypes were imputed using Beagle (version 4.0) [46]. Exploration of imputation accuracy with various missing data rates was done via the comparison between imputed calls and randomly masked genotypes (~ 5000,000 for site × individuals) for chromosome 10 only. The best parameter combination for the present study was determined to be as follows: window = 50,000, overlap = 5000, original missing rate < 75% (sites from the 200K array were always kept), and using the imputation-then-removing strategy of [47]. With this method, the average imputation accuracy was 98%. Ultimately, a set of 14,126,424 high-confidence SNPs and 439,670 InDels were retained for further analysis, covering about 99.13% of predicted maize genes (38,805 genes).
Mosaic map tracing IBD origins from 24 founders for each progeny line
The hidden Markov model (HMM) raised in [13] was used to make a multipoint probabilistic reconstruction of the genome of each progeny line as a mosaic of the founder identity-by-state (IBD). Each progeny genome is made up of IBD segments of the founder genomes, with a transition between founders occurring whenever a recombination has occurred (Additional file 3: Figure S5). The biallelic SNPs cannot distinguish between all founders, so the HMM used the information of neighboring SNPs to compute the posterior probability, \( {P}_{s,f}^i \), that at a given SNP locus s, the CUBIC progeny line i is descended from founder f. The IBD state is built as the founder with the maximum posterior probability \( {\max}_f\left\{{P}_{s,f}^i\right\} \) at SNP locus s for progeny i, only if the maximum posterior exceeds twofold of expected value by chance (1/24); otherwise, it was considered an unknown state. Segments that failed to be traced to any specific founder covered approximately 3.1% of the maize genome. In these regions, the identity-by-state between multiple founders was highly comparable compared to the flanking regions, presumably reflecting the co-ancestral origins.
Upon the CUBIC design, the HMM made the following approximations and assumptions: (i) the genome of each progeny line is homozygous due to advanced selfing generations, so we filtered out the SNPs of high heterozygous genotypes (> 10%) or missing calling before imputation (> 25%), leaving ~ 240K high-quality SNPs for IBD reconstruction; (ii) the effective number of generations (G) approximately be 9, since the eight rounds of inter-cross generation (dialle cross or open pollination), plus one round of approximate cross-generation from six generations of selfing, because on average, only half recombination events are accumulated per selfing [48]; (iii) the identity of the founder in a given IBD segment in the mosaic is uncorrelated with other segments for that individual, and the length of segment follows an exponential distribution with a mean length ρ/G, where ρ is the genetic length of the segment, corresponding to a consensus map (Additional file 14) built with ten bi-parental linkage maps [15].
The power of the HMM method in correctly tracing the IBD origins of CUBIC progeny lines was evaluated by simulating pseudoprogeny lines by reshuffling 24 founder genomes. First, we set the number of recombinant segments of the progeny line as 180 (the median of recombination per line in real data) and to be proportional to the chromosome length. Second, we set the location of the break point based on the empirical recombination ρ from the consensus genetic map (Additional file 14). The 24 founders were then randomly assigned to the simulated segments for each line. The procedure repeated 100 times to simulate 100 pseudoprogeny lines. The simulated recombination count of the 100 pseudoprogeny lines appeared to apparently decreased recombination in centromere relative to the arms of chromosomes (Additional file 3: Figure S19a). The SNP genotypes of 100 pseudoprogeny lines were projected from the founder SNP genotypes in every simulated IBD segment. Then, for each pseudoline, the HMM with the same criteria in real data was run to trace back to the 24 founders. The power was defined as the proportion of SNPs across the genome that have the identical IBD origin to what it was simulated to be. It was found that the current procedure of the HMM approach had the power of 93–97% (along different chromosomes) to correctly identify the founder origins (Additional file 3: Figure S19b).
Single-variant-based GWAS
In total, 11.8 million high-quality SNPs (MAF ≥ 0.02) were used to perform sGWAS. The relatedness matrix (K, random effect) was calculated using the “pylmmKinship.py” script from pylmm [49], and the top ten PCs (which explained 8.76% of the variance) were generated by GCTA [50] and applied as a fixed effect. The mixed linear model (MLM) corrected by a relatedness matrix and the top ten PCs was implemented by the software Tassel 3.0 [51]. The significance threshold for the association was set to 1.23E−8, which was equal to 0.05/Ne, where Ne is the effective number of independent tests [52]. To interpret GWAS results, significantly associated SNPs for each trait were first grouped into one locus in which two consecutive SNPs were less than 20 Kb, and in which at least two significant SNPs existed for each locus. The adjacent loci were further merged into a single locus if any significant SNPs between adjacent loci were in LD (r2 ≥ 0.2). To reduce the probability of missing causal genes in some narrow QTL intervals due to minor QTL, another 25 Kb was extended to both sides for those QTL with intervals less than 50 Kb and significance of the peak SNPs smaller than 100 times cutoff value (i.e., 1.23E−10). Ultimately, the significant loci were treated as sQTL, the peak SNP defined the significance of the sQTL, and the extended region of significant SNPs were defined as the sQTL interval.
IBD-based GWAS (hGWAS)
The IBD mosaic map of all CUBIC lines was collapsed into 27,005 bins based on all identified recombination break points. In each bin, there was only 1 IBD state for a given line but 16~24 IBD states available across all lines. Treating each bin as 1 variable, the hGWAS was performed via the “JLM” script described previously [53]. In hGWAS, the mixed linear model was built and the restricted maximum likelihood (REML) was used to test the significance of each bin, by treating bin and polygenic effects as random effects and the top 10 PCs as fixed effects. The covariance structure of the polygenic effects was inferred with a bin-based kinship matrix, which was calculated using the method raised by [54] through reformatting bins as the dummy variables. The general analytic differences for the present population to previous study [53] included 2 aspects: (1) the fixed effect, the previous method (JLM) intended to jointly analyze the data with 10 independent populations; thus, a categorical variable or design matrix was included as fixed effect to indicate sample-population relationships; instead, the populations structure of CUBIC was inferred using PCA, which was fitted as fixed effects in hGWAS model (similarly what we do in sGWAS); (2) random effect, it was modeled by a kinship matrix under linear mixed model framework. Previously the kinship matrix was simply calculated by genome-wide SNP data, instead for CUBIC population, we inferred it using the whole-genome bin marker data. We used a permutation test of 500 iterations to determine the threshold for the likelihood ratio test (LRT) scores for each trait [55]. The resulting LRT threshold ranged from 6.8 to 7.4 (α = 0.05); for simplicity, we chose the average LRT of 7.1 as the overall cutoff for all traits. Following the previous procedures, significant bins were merged into a locus (or hQTL) with nearby positions (≤ 1 Mb) or located within ≤ 5 bins; the interval of each hQTL was defined as the physical position range delimited by the significant bin [53].
GWAS simulation
To test and interpret the difference between the 2 GWAS methods, we simulated 3 possible QTL architecture scenarios using the CUBIC population, using biallelic quantitative trait nucleotide (QTN) expression as an SNP marker. These were as follows: (1) biallelic QTL, expressed as a single QTN at a given QTL; (2) 4-allelic QTL, attributed to 2 incompletely linked QTN at a given QTL; and (3) 9-allelic QTL, attributed to 3 incompletely linked QTN at a given QTL. We randomly selected 30 bins on chromosome 1 and set 10 bins as harboring single QTN for each bin; another 10 bins were assigned 2 QTNs per bin; and 3 QTNs assigned to each bin for the last 10 bins. We added a series of modest additive effects to QTNs, following the geometric distribution of an (a = 0.95), to simulate the polygenic nature of the complex traits [56]. The narrow-sense heritability (h2) and phenotypic variance (Vp) were defined as 0.8 and 1 for the target trait, respectively. The simulated phenotype of each line was obtained by adding the sum of additive values to a residual error, following a normal distribution with 0 mean and variance of 0.2 (i.e., h2 = 0.8). The 2 GWAS methods described above were used to detect QTL in this simulation, which was replicated 1000 times to evaluate the statistical power and FDR for each QTL scenario.
IBD-based QTL refinement by a linkage-like analysis
The major QTLs (those simultaneously detected by both GWAS methods with > 10% variance explained for the trait under study) were validated and refined. The major QTL region was defined as the union of the sQTL and hQTL intervals. The progeny lines were categorized into parental IBD groups according to the peak bin. Inferred from the simulation analysis that compared the statistical power of sGWAS and hGWAS, it may be a fact that multiple parents may share the same QTL allele. Thus, identifying or approximating the true QTL allele can reduce the model complexity (decrease the degree of freedom), which should increase the statistical power. But note that the IBD clustering for a QTL has to be phenotype dependent, which can reduce the impact of random variants on the clustering. Thus, for a QTL, the parental IBD groups were collapsed into phenotypically distinguishable clusters based on multiple comparisons of phenotype between the parent IBD groups (P < 0.05) for each trait. These IBD clusters, called functional alleles, reflected potentially true alleles dependent on specific traits, thus expecting to more closely predict trait variation than the parent IBD, independent of traits. At the QTL region, the inferred allelic types were projected onto all progeny lines, and the allelic combinations at the QTL region were re-clustered for the whole population into several major haplotypes. In this way, the pairwise comparison of phenotype between inferred haplotypes (as determined using the Student’s t test, P < 0.01) enabled the verification of the QTL credibility and further narrowed the QTL interval, in a manner similar to the traditional bi-parental fine-mapping procedure [57].
For fine mapping major GWAS loci, SNP and InDel polymorphisms located in the QTL interval were used to enable identification of the gene and the causal variant responsible for the GWAS identified locus. To help determine the likely functional genes, we ranked all genes via functional annotation (predicted by ensemble VEP program [58]) of polymorphisms located within the QTL interval, based on the generic feature format file version 3 (gff3) in maize following the listed procedure [59].
Identification of epistasis in trait variance
Variants were filtered for frequency (MAF ≥ 2%) and linkage disequilibrium (r2 < 0.5 within each window of 50 Kb, with a step of five SNPs) before testing for epistasis. For each trait, the plink program [60] was employed in fast epistasis analysis using the supported boost test [61] with --fast-epistasis. Pairs of loci with P value > 6E−16 (roughly 0.01/Ne2, where Ne is 4,057,944, the independent number of SNPs estimated by Genetic type 1 Error Calculator (version 0.2) [52]) were first removed, and only the most significant ones per pair were retained for those loci that interacted with multiple other loci within 100 Kb. The P value and corresponding variance explained for each putative epistatic pair were further adjusted by linear regression with controlling of population structure and additive effect, and those with adjusted P values < 1E−12 were kept as instances of significant epistasis. To estimate the trait variance jointly explained by epistasis, the population structure and additive effects for all remaining variants were first regressed out, and the residuals for each trait were further regressed against all interacting items. An extended interval of 50 Kb for interacting variants was used to determine the co-occurrence of identified epiQTL with previously identified networks.
RNA sequencing, eQTL mapping, and epQTL analysis
A subset of 391 progenies were randomly selected from the CUBIC population for RNA sequencing. These lines and the founder parents were grown at the Hainan field station in the winter of 2016. At the V9 stage (the stage with the fastest leaf tissue growth), total RNA was extracted from the tissue of the 11th leaf, collected from the pool of 3 plants of similar growth status per line, following the standard protocol of Quick RNA isolation Kit, #0416-50 (Huayueyang Biotechnology CO., LTD. Beijing, China). A library with insert sizes ranging from 200 to 500 bp was prepared using the commercial Illumina library preparation kits (TruSeq Stranded mRNA LT-SetA. RS-122-1201). A 150-bp paired-end Illumina sequencing was performed using the HiSeq X-Ten protocols. Each sample had on average ~ 20 million raw reads. The reads with low sequencing quality and sequencing adapter were removed using the software Trimmomatic-0.36 [40]. The paired-end reads were mapped onto the B73 AGPv3.25 reference using the software STAR [62] with a maximum intron size of 50,000 bp; only those uniquely mapped reads were used to quantify gene expression levels using the HTSeq [63]. The expression data for each gene was normalized using the software DESeq2 [64] before the subsequent analyses.
Only genes expressed in more than 60% of the lines were retained in eQTL mapping. Top ten PEER [65] factors, together with the top ten genotypic PCs, were utilized to account for covariates. The software EMMAX [66] was used in eQTL mapping, taking 1.89E−8 as a significant standard. The co-localization of pQTL and eQTL was measured with the evidence that the peak variant of pQTL was included in significant SNPs from eQTL. Those loci affecting trait variance by regulating those genes outside the pQTL region and not due to LD (> 10 M) were considered candidates with trans-effect. The OSCA program [32] was used to perform expression-phenotype associations (epQTL) under a linear mixed model, with 1E−3 as a significant threshold.
Cross-population analysis
The mapping on a maize association mapping panel (AMP) containing 508 unrelated inbred lines [67] was used to cross-validate and narrow down the candidate list. Since the mapping results of AMP were based on the B73 v4 genome, we first extracted the candidate genes located in 776 QTL in the present CUBIC population and transformed the B73 v3 gene ID into B73 v4 gene ID. All CUBIC QTLs were applied to candidate association mapping in AMP along corresponding traits and considered as significant when P value < 1E−4. Meanwhile, the cross-validation in the maize association population makes the number of candidate genes greatly reduced.
Metabolic analysis
A metabolic profile was conducted to help interpret how genes identified in this study may be involved in regulating the phenotype. From the 391 lines used in RNA-seq, a subset was selected based on the phenotype and inferred allelic types in the specific QTL under study. For the major QTL identified for ELW on chromosome 4, a total of 71 lines were chosen based on the genotype of peak SNP, with 35 vs. 36 for different alleles of the 1 bp deletion in ZmGalOx1. These lines and the founder parents were grown at the Hainan field station in the winter of 2016. At the V9 stage, we collected tissue from the 10th leaf, collected from 3 plants per genotype, for metabolic analysis. The 7th leaf of the CRISPR-Cas9 edited T0 plants was also collected at the V9 stage from plants grown in the greenhouse of the China National Seed Group Co., LTD. (Wuhan, China), in the autumn of 2017. The primary metabolites were extracted from maize leaves following the procedures described in detail in previous studies [68, 69]. In brief, 50 mg of lyophilized leaf powder was extracted with 1 ml of methanol: methyl-t-butyl-ether (MTBE) solution (1: 3 v/v). A 300-μl aliquot from the lower polar phase was taken and dried in a vacuum. The dried residue was derivatized with N-methyl-N-(trimethylsilyl) trifluoroacetamide (MSTFA) and further analyzed by GC-MS (7890A-5975C, Agilent, Santa Clara, USA) following the protocol described in previous studies [70, 71].
Hormones were extracted from maize leaves according to the protocols published previously [72], and the hormone compounds were separated by reversed-phase ultra-fast LC (Shimadzu, Kyoto, Japan). These were detected by electrospray ion source of a tandem triple quadrupole MS analyzer (API4000, AB SCIEX, Singapore) and quantified in multiple reaction monitoring (MRM) mode using optimized MS/MS conditions. The MS conditions were as follows: source, Turbo IonSpray; ion polarity, negative; IonSpray voltage, − 4500 V; source temperature, 550 °C; gas, nitrogen; curtain gas, 30 psi; nebulizing gas (GS1), 55 psi; collision gas (GS2), 55 psi; scan type, MRM; Q1 resolution: unit; and Q3 resolution: unit. The Analyst 1.5.2 software (AB SCIEX, Foster City, CA, USA) was used to control the instrument and to acquire and process all MS data. Among all samples, ten samples were tested twice for technical replications. The Pearson correlations between two replications indicated the high quality of metabolic data (r = 0.9927 ± 0.0083).
CRISPR/Cas9 editing experiment
In order to confirm that ZmGalOx1 is the causal gene for the ear leaf width QTL identified via GWAS on chromosome 4, we edited the gene sequence using the CRISPR-Cas9 system. Two guide RNAs (gRNAs) targeting the 14th exon of ZmGalOx1 were designed using the CRISPR-P [73]. The 2 gRNAs were integrated into a highly efficient vector [74] and transformed into immature embryos of ZZC01 by Agrobacterium infection by the China National Seed Group Co., LTD. A total of 55 T0 CRSPR-Cas9 edited seedlings were obtained (including positive and negative control plants), and DNA was extracted from each individual at the seedling stage. Primers were designed to amplify about 800 bp in the vicinity of (and including) the 2 gRNAs within the gene. The PCR products were used for Sanger sequencing to examine the variations along the sequence. By comparing the sequences with those of the unedited lines (ZZC01) using the BioEdit software [75], all sequence changes that could lead to amino acid changes were identified. For all 55 T0 plants, ear leaf widths in the R1 stage (after pollination) were measured for statistical analysis.
Leaf cytological experiment
Lines with the most extreme ear leaf width, and having either the allelic type with or without the 1-bp insertion in ZmGalOx1, were planted in Wuhan in the summer of 2017. Ear leaves were collected in the R1 stage to prepare for the cytological experiment. In the autumn of 2017, the ear leaves of CRISPR/Cas9 edited plants (including positive and negative control plants) were also collected at the R1 stage from greenhouse-grown plants. The abaxial leaf epidermis tissue was removed 20 cm from the leaf tip and affixed to a microscope slide. Images were taken under × 10 field in an Olympus DP72 compound microscope with × 10 eyepiece. Counting in a direction perpendicular to the leaf vein, the number of cells per field was counted by the Image-Pro Plus software [76] to determine the average width of a single cell. The total number of cells along the leaf from left to right was calculated as leaf width divided by the average single cell width. For each plant group, 12 plants were investigated for cell size and cell number with 3 leaves per plant. The 3 values were averaged for each plant. Student’s t test was used to test the significance of difference between the positively and negatively edited plants and between different ZmGalOx1 genotypic plants (P < 0.01).