The genetic architecture of leaf stable carbon isotope composition in Zea mays and the effect of transpiration efficiency on leaf elemental accumulation

Abstract With increased demand on freshwater resources for agriculture, it is imperative that more water-use efficient crops are developed. Leaf stable carbon isotope composition, δ13C, is a proxy for transpiration efficiency and a possible tool for breeders, but the underlying mechanisms effecting δ13C in C4 plants are not known. It has been suggested that differences in specific leaf area (SLA), which potentially reflects variation in internal CO2 diffusion, can impact leaf δ13C. Furthermore, although it is known that water movement is important for elemental uptake, it is not clear how manipulation of transpiration for increased water-use efficiency may impact nutrient accumulation. Here, we characterize the genetic architecture of leaf δ13C and test its relationship to SLA and the ionome in five populations of maize. Five significant QTL for leaf δ13C were identified, including novel QTL as well as some that were identified previously in maize kernels. One of the QTL regions contains an Erecta-like gene, the ortholog of which has been shown to regulate transpiration efficiency and leaf δ13C in Arabidopsis. QTL for δ13C were located in the same general chromosome region, but slightly shifted, when comparing data from two different years. Our data does not support a relationship between δ13C and SLA, and of the 19 elements analyzed, only a weak correlation between molybdenum and δ13C was detected. Together these data add to the genetic understanding of leaf δ13C in maize and suggest that improvements to plant water use may be possible without significantly influencing elemental homeostasis.


Introduction
The impacts of global population growth and climate change on natural resources indicate that the future of food security will depend on increasing both the productivity and sustainability of agriculture systems (National Academies of Sciences, Engineering, and Medicine 2018). Improving crop water-use efficiency (WUE) would ameliorate the effects of the increasing frequency and severity of droughts (Sheffield and Wood 2008;Chapman et al. 2012;Leakey 2019) . Agronomic WUE can be defined as the amount of yield, whether grain or biomass, produced per the total amount of water utilized by the crop (Condon et al. 2004). Many factors can affect WUE including transpirational water loss through the stomatal pores on the leaf's surface. In C 3 plants the amount of carbon available for assimilation is limited by stomatal and mesophyll conductances to CO 2 (Flexas et al. 2016) and therefore correlated to the rate of transpiration. For example, yield was shown to be positively associated with cumulative transpiration in soybean (Purcell et al. 2007), and higher net carbon assimilation was accompanied by higher transpiration in rice (Adachi et al. 2017). However, higher rates of biomass yield do not always correspond to higher transpiration rates in C 4 plants due to the evolution of the carbon concentrating mechanism. The uncoupling of CO 2 assimilation and transpiration has been demonstrated in the field and greenhouse-grown maize (Walker 1986;Kolbe et al. 2018a). Thus, there is the potential to increase transpiration efficiency, or carbon gain per amount of water transpired, without reducing productivity in C 4 species (Leakey 2019). A large amount of variation is present in the transpiration rates of C 4 crop species, including sorghum (Hammer et al. 1997), and maize (Bunce 2010), suggesting that existing occurring alleles could be exploited for optimizing WUE.
Although increasing transpiration efficiency provides a strategy to avoid the negative effects of water limitation on plant growth and development (Passioura 1996;Chaves et al. 2002;Jaleel et al. 2009), there is the possibility of pleiotropic side effects given the fundamental requirement for water movement in plants. A potential impact of reducing transpiration could be a corresponding reduction in the uptake and mobilization of watersoluble nutrients. As water is absorbed by roots, nutrients in solution come in contact with the root surface in a process known as mass flow (Barber et al. 1963). Most nutrients are acquired by mass flow, although phosphorus is a notable exception that contacts the root through diffusion (Barber 1962). Reducing transpiration may also affect nutrient uptake facilitated by symbiosis with mycorrhizal fungi (Marschner and Dell 1994). Therefore, the manipulation of basic plant processes such as transpiration for improved WUE must also consider potential impacts on the availability of essential plant nutrients. Essential elements contribute to cellular function in numerous ways including in biochemical reactions as catalytic cofactors, for charge balance in cellular and subcellular compartments, as well as in DNA and protein, the building blocks of life (Baxter 2009). Previous research has shown in the C 4 plant sorghum that total leaf mineral content is positively correlated with transpiration efficiency (Masle et al. 1992). While meta-analyses of high CO 2 grown plants with reduced transpiration have shown a drastic reduction in nutrient accumulation in several tissues including the leaves and grain of C 3 crops (McGrath and Lobell 2013), sorghum showed no difference and maize had similar levels of zinc, protein, and phytate, but a decrease in iron accumulation in the grain (Myers et al. 2014). Although part of the reduction in nutrient content can be explained by dilution, due to increased growth at high CO 2 , this does not completely account for the observed reduction. An ionomics (high-throughput elemental profiling) approach has been used in maize to assess kernel nutrient content (Baxter et al. 2014). A similar ionomics approach in leaf tissue could be used to assess the effect of transpiration on nutrient accumulation in the leaves.
The difficulty and labor-intensive nature of accurately quantifying the amount of water that an individual plant transpires have been a major limitation to breeding for transpiration efficiency. This has resulted in the selection for drought tolerance rather than applying a direct selection for water use (Cooper et al. 2014;Gaffney et al. 2015). One alternative method is the use of leaf stable carbon isotopes as a proxy for transpiration efficiency. The stable carbon isotope composition, d 13 C, reflects the amount of 13 C present in plant tissue relative to a standard (Keeling 1979). Enzymes in the process of carbon fixation discriminate differently against the heavier 13 C atoms in a process known as fractionation (Farquhar et al. 1982;O'Leary 1988). It has been widely shown that stable carbon isotopes can be used as a proxy trait for quantifying a plant's transpiration efficiency in C 3 plants (Farquhar et al. 1989a(Farquhar et al. ,1989bCondon et al. 1990;Virgona et al. 1990;Condon et al. 1993;Barbour et al. 2010) and in C 4 plants (Henderson et al. 1998;von Caemmerer et al. 2014;Ellsworth et al. 2017Ellsworth et al. , 2020Twohey III et al. 2019). Due to the differences in the photosynthetic pathways, d 13 C is positively correlated with WUE in C 3 species, but negatively correlated with d 13 C in C 4 species. Studies have also shown that d 13 C can be influenced by environmental factors such as light intensity and drought, as well as anatomical traits and bundle-sheath leakiness to CO 2 (reviewed in Cernusak et al. 2013). However, genes underlying variation for d 13 C remains unknown in C 4 species. Kolbe et al. (2018b) showed that leaf d 13 C did not correlate with the leaf activity of any of the photosynthetic enzymes previously posited to influence d 13 C variation. In addition, a transcriptome analysis was unable to identify a clear candidate gene (Kolbe et al. 2018b). Quantitative genetic approaches have the potential to reveal the genetic control of d 13 C in C 4 species because genomic locations are tested for associations with the trait of interest, without a priori knowledge of the mechanism underlying the variation. Maize is ideal for use in mapping studies due to its high level of recombination and low linkage disequilibrium (Yu and Buckler 2006). Mapping methods have been successfully used for decades to identify genes controlling complex traits in maize, with evolving approaches to tackle more difficult traits (Wallace et al. 2014). In addition, maize is both a model organism with available populations and genomic data, and one of the three most important global crops contributing to 30% of the total calories consumed by humans (Shiferaw et al. 2011).
There have been several previous studies that have used quantitative genetics to investigate d 13 C in C 3 species (Teulat et al. 2002;Masle et al. 2005;Rebetzke et al. 2008;Xu et al. 2009). In Arabidopsis, the gene ERECTA was identified in a QTL study for isotopic discrimination and was found to alter transpiration efficiency by altering stomatal density (Masle et al. 2005). Genetic mapping of leaf d 13 C has also been performed in the C 4 species Setaria viridis (Feldman et al. 2018;Ellsworth et al. 2020) and kernel d 13 C has been mapped in the C 4 maize (Gresset et al. 2014;Avramova et al. 2019). Although the QTL found for C 4 species still require fine-mapping to identify the causative gene, no correlation was observed between kernel d 13 C and leaf d 13 C (Foley 2012). The lack of correlation may be the result of post-photosynthetic fractionation (Badeck et al. 2005;Cernusak et al. 2013), and therefore mapping QTL for d 13 C in leaves may reveal additional loci not found using kernels. In this study, we focus on leaf d 13 C in maize and its association with leaf elemental composition. We leverage results from previous studies to select biparental mapping populations. Specifically, the NAM founder lines CML103, CML333, and Tx303 consistently contrast B73 with respect to leaf d 13 C (Kolbe et al. 2018b;Twohey III et al. 2019). The founder line NC358 had a moderate leaf d 13 C value (Kolbe et al. 2018b;Twohey III et al. 2019) and was also included in this study. In addition to leaf d 13 C, we also investigated variation in specific leaf area (SLA) and its potential relationship to leaf d 13 C by CO 2 diffusion. Characterization of the genetic architecture of leaf d 13 C will provide a better context for understanding what drives d 13 C, which will allow breeders to utilize this trait in crop improvement.

Plant material
All experiments were planted at the University of Illinois Crop Sciences Research Farm, Urbana, IL, USA and were subject to natural conditions without supplemental irrigation. Plants did not show signs of stress prior to sample collection. Weather data collected less than 1.5 km from the field site can be found in Supplementary Table S5 (Water and Atmospheric Resources Monitoring Program, Illinois Climate Network 2020). NAM RIL families CML103, CML333, NC358, and Tx303 and NAM founder parents (McMullen et al. 2009) are publicly available through the Maize Genetic Cooperative Stock Center. The NAM RIL families were planted in the summer of 2015 and a subset was planted in 2019 using an augmented incomplete block design. Fifteen kernels were planted in each 3.7 m row with 0.8 m spacing between rows and 0.9 m alleys. The families were randomized together, with each block consisting of 20 lines and 2 checks (B73 and one of the other founder lines). Ten percent of the plots were dedicated to checks with the common parent B73 appearing in each block with one of the four founder lines. All lines used for the GWAS experiment are publicly available through the USDA Germplasm Resources Information Network (GRIN). This experiment was planted in the summer of 2016. Twenty kernels were planted in each 3.7 m row with 0.8 m spacing between rows and 0.9 m alleys, and then thinned to 15 plants per row. The reference line B73 was replicated 22 times within the experiment. A complete list of lines used can be found in Supplementary Tables S1-S3.

Tissue sampling
Samples for d 13 C analysis from the 2015 NAM RIL populations were collected 6 weeks after planting (approximately V9) as follows. A rectangular piece of tissue approximately 7.5 cm Â 5 cm was taken from the middle of the leaf blade (excluding the midrib) of the uppermost fully expanded leaf from four plants in each plot. Samples were placed in a coin envelope and dried at 65 C for at least 7 days. After drying, four-hole punches (each 0.058532 cm 2 ) were taken and placed in a 6 mm Â 4 mm tin capsules (OEA Laboratories # C11350.500P) for analysis using a Delta PlusXP (Washington State University) isotope ratio mass spectrometer. Leaf samples for SLA measurements were collected from four plants in each of the plots (preferentially but not necessarily the same plants as were collected for d 13 C) using a 1.6 cm diameter cork borer. Leaf discs were dried at 65 C for at least 7 days prior to weighing on an analytical balance (Model MS204S). SLA was calculated as the area of a leaf disc divided by its dry weight. These same leaf discs were then used for ionomics analyses as described in Pauli et al. (2018). Leaf samples for d 13 C analysis from the GWAS panel and the 2019 NAM RIL population were collected at approximately V10 (6-7 weeks after planting) using the hole punch method and processed as described previously (Twohey III et al. 2019). Due to the high level of diversity in this panel, some lines were flowering when samples were collected, which resulted in tissue being collected from the flag leaf. These samples were analyzed using a Costech instruments elemental combustion system and a Delta V Advantage isotope ratio mass spectrometer.

Statistical analysis
All analyses were completed using custom scripts (available upon request) and statistical packages in R (R Core Team 2017). A mixed model approach was utilized to obtain the best linear unbiased estimation (BLUE) of the fixed genotypic effect (Henderson 1975). For both years 2015 and 2019, the experimental design of the RIL panel was blocked containing partial repeated entries of the parental genotypes, B73 and CML103. A mixed-effect model was created using the R package "lme4" to account for the variation among repeated entries and to estimate d 13 C for each genotype (Bates et al. 2015). In the model, genotypes were treated as fixed effects with year and block terms treated as random effects. In addition, the block term was nested within year as the genotypes were randomly assigned to blocks each year. The R package "emmeans" utilized Least-Square Means to provide BLUEs for each genotype (Searl et al. 1980).

Correlation analysis
Pearson correlations using phenotype mean values were calculated with corr.test() in R package "psych" (Revelle 2017) using complete observations and Holm's method (Holm 1979) to adjust P-values for multiple testing. The correlation matrix was visualized using pairs.panel() in the R package "psych" (Revelle 2017).

Single family QTL mapping
The analysis was completed using NAM_phasedImputed_1cM_ AllZeaGBSv2.3 dataset. The file contains fully imputed and phased genotypes for most of the RILs in the NAM population (Zhao et al. 2006;Lipka et al. 2015). This HapMap format file was converted to numeric format where 0 is the B73 homozygote reference, 1 is a heterozygote, and 2 is the homozygote alternative parent. Phenotypic means were regressed onto genotype. Lowest P-values from the ANOVA values of the linear model were recorded (i.e., pvalues[i] ¼ anova(lm(mypheno$geno[i , ])). The previously identified marker was added to the model and re-run in a stepwise regression procedure. The final model included all identified QTL. Significance thresholds were determined by 200 permutations and alpha was set at 0.05. All analysis was completed using custom scripts in R (R Core Team 2017). Results were then compared to composite interval QTL mapping completed in R package "r/QTL" (Broman et al. 2003). QTL locations are indicated using B73 version 2 positions.

Joint linkage mapping
The analysis was completed using HapMapv2 (Chia et al. 2012). The genotypic dataset consisted of 836 markers were scored on 624 RILs from four biparental families with B73 as a common parent. The marker subset is composed of markers that could be placed unambiguously on the physical map, previously described in Brown et al. (2011). Unambiguous markers are defined by those anchored in CDS positions of genes that have held consistent over genome versions verified by MaizeGDB cross-reference tables. The markers are approximately evenly spaced across the genome with an average spacing of 1.6 cM. Missing data were imputed as previously described in Tian et al. (2011). Joint linkage models were constructed using custom script in R (R Core Team 2017) by a stepwise regression procedure. In general, we used linkage to test every marker across all four families to find the most significant QTL. The model has a family term and a marker: family term. The family term accounts for differences in mean phenotype between families. Inclusion of the marker: family term means that for each QTL we are assigning a separate effect to each family. The family term was included in the model and each of the 836 possible marker-by-family terms were assessed following the method of Ogut et al. (2015). The lowest P-values from the ANOVA values of the linear regression model were recorded (i.e., JL_pvalues[i] ¼ anova(lm(my_pheno$familyþgeno[, i]: family)). All 836 marker-by-family terms were tested. SNP effects were nested in families to reflect the potential for unique QTL allele effects within each family. Significance thresholds were determined using 1000 permutations for each family independently and alpha was set at 0.05. The lowest resulting P-value was recorded for each permutation.  Supplementary Table S2. Hirsch et al. (2014) collected RNA from whole seedling tissue which was sequenced via IlluminaHiSeq and filtered to create a working set of 485,179 SNPs that is available at https://datadryad.org//resource/ doi : 10.5061/dryad.r73c5. The 413 lines were grown in 2016 and tissue was sampled when B73 was at the developmental stage V10. Isotopic analysis is described above. A genome-wide association analysis was run using R package "GAPIT" (Lipka et al. 2012) on leaf d 13 C. Removal of SNPs with a minor allele frequency of less than 0.05 resulted in a subset of 438,222 SNPs being used in this analysis. A MLM model was used with model selection set to true to find the optimum number of principal components to account for population structure (Lipka et al. 2012). Significance thresholds were calculated using the Bonferroni correction of familywise error rate. An alternative significance test was calculated using the Benjamini-Hochberg procedure for controlling the false discovery rate (Benjamini and Hochberg 1995).

Data availability
Genotypic datasets were downloaded from Panzea CyVerse iPlant Data Storage Commons (http://datacommons.cyverse.org/ browse/iplant/home/shared/panzea). All phenotypic datasets were quality controlled for complete technical replicates and availability of genotypic data. A list of all genotypes used in each analysis is provided in Supplementary Tables S1-S4 and have been uploaded to figshare. Briefly, the d 13 C analysis was completed with 640 RILs; including 156 CML103 RILs, 160 CML333 RILs, 159 NC358 RILs, and 165 Tx303 RILs (Supplementary Table  S1). The element analysis was completed using a total of 704 RILs; including 175 CML103 RILs, 181 CML333 RILs, 175 NC358 RILs, and 173 Tx303 RILs (Supplementary Table S1 Figure S1 shows the distribution of leaf d 13 C for each of the NAM RIL families. Supplementary Figure S2 presents the correlation matrix for the elemental analysis, and Supplementary Figure S3 shows the chromosomes where significant QTL were identified for each element. Supplementary Figure  S4 is the LOD plot from the GWAS mapping of leaf d 13 C. Supplementary material is available at figshare: https://doi.org/ 10.25387/g3.14114234.

Single-family QTL mapping
The RIL families generated from four NAM founder lines (CML103, CML333, NC358, and Tx303) were grown for linkage analysis. Consistent with previous studies, both the CML103 and CML333 parent lines had a significantly less negative leaf d 13 C than B73 (P < 0.05), when grown as replicated controls among the RILs. However, the Tx303 and NC358 parental lines were not found to be significantly different from B73. Transgressive segregation was observed in all four RIL families (Supplementary Figure S1).
Stepwise regression analyses found significant QTL for leaf d 13 C in the NAM RIL families CML103, CML333, and Tx303 but not in NC358 ( Figure 1A). Interestingly, none of these QTL were shared between RIL families in this analysis. The CML103 RIL family had two significant QTL, one on chromosome 5 at 211.7 Mb and another on chromosome 7 at 142.4 Mb. Combined these two QTL accounted for 21.36% of the total phenotypic variance (R 2 ¼ 0.2136; Table 1). The RIL family CML333 had one significant QTL on chromosome 3 at 183.9 Mb, which accounted for 8.37% of the total phenotypic variation explained (Table 1). Finally, the Tx303 RIL family had a significant QTL on chromosome 2 at 13.5 Mb explaining 9.48% of phenotypic variation (Table 1). No significant QTL for leaf d 13 C were identified in the NC358 RIL family.
SLA was used as a proxy trait to test for a relationship between leaf thickness and leaf d 13 C. No significant correlation was observed between SLA and leaf d 13 C (P ¼ 0.304, r ¼ À0.0414). In addition to testing for a correlation with leaf d 13 C, QTL mapping was performed for SLA to identify any possible overlaps with genomic regions identified for leaf d 13 C. Mapping of SLA in the four RIL families identified two significant QTL. In the CML103 RIL family, a QTL was identified on chromosome 5 at 86.1 Mb and in the Tx303 RIL family a QTL on chromosome 9 at 107.8 Mb. Neither of the SLA QTL identified overlapped with QTL for leaf d 13 C ( Figure  1C).
To test a potential link between transpiration and nutrient uptake, an elemental analysis was performed on leaf samples from each of the four RIL families. Samples were analyzed for 19 different elements using an ICP-MS. A full correlation matrix shows that some elements are highly correlated with each other (Supplementary Figure S2), but no strong correlations (r > 6 0.7) were identified with d 13 C. However, there was a weak but significant correlation (P ¼ 6.745E-05, r ¼ 0.18) between leaf d 13 C and Mo (Figure 2). Subsequent QTL mapping of the 19 element concentrations identified 28 QTL across 12 different elements (Figure 3, Supplementary Figure S3). Significant QTL were found for B, Mg, P, S, K, Fe, Mn, Co, Cu, Rb, Sr, and Mo (Supplementary Table S4). None of the elemental QTL overlapped with those found for leaf d 13 C or SLA. However, in some cases multiple elements had common QTL, such as Mg and Mn on chromosome 10 in the CML103 RIL family and Co and Cu on chromosome 3 in the NC358 RIL family. In addition, common QTL for an element were found across families, as in the case of Mg in the NC358 and Tx303 RIL families.

Joint linkage QTL mapping
A joint linkage analysis was performed for leaf d 13 C to test whether any additional QTL would be identified by combining the four RIL families into a single analysis. The joint linkage analysis identified the same significant QTL for leaf d 13 C on chromosomes 2, 3, and 5 (Table 2) as in the single-family stepwise regression analysis. Although the QTL on chromosome 7 was not found using the joint linkage approach, an additional significant QTL was identified on chromosome 1. Given that no significant QTL for leaf d 13 C were identified in the NC358 RIL family, we tested whether removing this family from the joint linkage analysis would change the outcome. When the joint linkage analysis was rerun excluding family NC358, the same four QTL were reidentified with decreased P-values, and the total phenotypic variation explained (R 2 value) increased in later steps of the model. However, no new QTL was identified with this approach.

QTL are consistent across years
One of the biparental mapping populations was grown in a second field trial to investigate the repeatability of the QTL between years. The CML103 RIL family had a significant QTL on chromosome 5 in both years ( Figure 1B). The peak of the QTL shifted slightly (<1 Mb) between years, such that the 1 LOD intervals did not overlap (Table 1). However, there was overlap when a 2 LOD interval was calculated. A Best Linear Unbiased Estimate (BLUE) was calculated using the blocking information to account for spatial effects in each grow out. Again a QTL was identified in the same region, although shifted from what was observed when analyzing each year independently. The second QTL observed on chromosome 7 in the 2015 CML103 RIL family was not significant in 2019, nor was it significant when using BLUE for mapping.

Genome-wide association study
Once significant QTL intervals were identified for leaf d 13 C using a biparental mapping strategy (Figure 1, A and C), we performed a genome-wide association study to try and narrow down the intervals to specific genic regions. The Wisconsin Diversity Panel was chosen because it represents a large portion of variation found within maize and has a robust publicly available 485,179 SNP set. A subset of 413 of the possible lines was chosen due to seed availability, and were grown in a single randomized block. No significant SNP associations with leaf d 13 C were identified (Supplementary Figure S4).

Discussion
Leaf d 13 C has a moderately high heritability in maize (Twohey III et al. 2019), which facilitates the use of quantitative genetics approaches to pinpoint the genomic locations controlling this trait. Here we characterized the genetic control of d 13 C in maize using leaf tissue collected at vegetative stage V9-V10 to reflect the photosynthetic pool during active growth. We were able to identify several significant QTL for leaf d 13 C across three NAM RIL families using stepwise regression. Using these populations we were also successful in identifying QTL for SLA and 12 different elements. Contrary to our hypothesis, no significant correlation was observed between leaf d 13 C and SLA or elemental composition.
We strategically picked NAM RIL families based on the founder parents that had the largest differences in their d 13 C for single family and joint linkage analyses. However, we also included a parent which was not extremely different from B73. Interestingly, we were unable to identify significant QTL in the NC358 RIL family despite transgressive segregation. We interpret  this result as an indication the NC358 contains only small effect QTL relative to B73 that were not detected in this study. Alternatively, NC358 leaf d 13 C may be more sensitive to the growing environment with a smaller genetic component. Twohey III et al. (2019) noted that while several maize lines were stable when tested in greenhouse and field environments, there were other lines that had highly variable isotopic signatures. A large amount of environmental influence over this trait in some backgrounds would obscure the genetic contribution and our ability to detect significant QTL. When we compared the regions identified here with regions previously mapped in S. viridis no obvious overlap was observed (Ellsworth et al., 2020). However, of the QTL identified for kernel d 13 C in maize (Gresset et al. 2014; Avramova et al 2019), our QTL for leaf d 13 C overlapped those on chromosomes 1, 3, and 7. This result demonstrates that some QTL for d 13 C may be shared between tissues, and that these QTL are identified across several populations and environments. Indeed, when the CML103 RIL family was grown over two seasons the major effect QTL on chromosome 5 was observed in both growing environments.
Although the QTL analyses presented here do not provide gene-level resolution, we were able to look for candidate genes within the intervals. The 2015 chromosome 5 QTL includes an Erecta-like gene (er1, GRMZM2G463904,211.8 Mb). However, er1 was located outside of the 1 LOD interval of the 2019 chromosome 5 QTL. Furthermore, mapping of BLUE further shifted the chromosome 5 QTL peak away from er1 (Table 1). Unfortunately, stomatal density data were not collected on these populations, which would further support the role of er1 in variation of leaf d 13 C. This would be an interesting area of future research given that this gene was found to effect d 13 C in Arabidopsis by changing stomatal density (Masle et al. 2005). We also looked for genes that have been previously shown to directly influence transpiration efficiency in maize (reviewed in Leakey 2019). However, none of these were found to be located in our QTL intervals.
The linkage analyses using biparental mapping populations identified several significant QTL, but none of the single-family QTL were independently identified in more than one family ( Figure 1). This result indicates that leaf d 13 C can be controlled by different factors depending on the genetic background. Furthermore, if in fact leaf d 13 C is controlled by many small-effect QTL, this may explain why the GWAS did not identify any significant SNP associations with leaf d 13 C. Identifying rare alleles with small to moderate effect size is a known weakness of the GWAS method (Bazakos et al. 2017). A better understanding of the mechanisms influencing leaf d 13 C would allow future analyses to move beyond single marker tests and instead look at SNPs in genes representing a particular pathway or process that could be collectively significant. This approach was successfully used to study maize lipid biosynthesis (Li et al. 2019a).
The diffusion of CO 2 into mesophyll cells is a potential source of variation in leaf d 13 C, which could be linked to stomatal density or leaf thickness. Previous work in maize has shown that stomatal density is not correlated with leaf d 13 C in a small diversity panel of maize (Foley 2012). SLA has not been linked to d 13 C in maize, but in rice d 13 C and SLA have shared QTL . In this study, we were able to test SLA and d 13 C in four RIL families, and no correlation was observed. Likewise, a comparison of the QTL analyses showed no overlapping regions for SLA and those mapped for leaf d 13 C. This result suggests two possibilities. First, it is possible that differences in SLA observed in these populations are not due to leaf thickness, but rather composition. Identifying the causative genes underlying the QTL would give insight into the mechanism. A second possibility is that leaf anatomical traits other than leaf thickness influence leaf d 13 C. A variety of anatomical traits could affect d 13 C and would not be directly captured by measurement of SLA, such as stomatal ratios, interveinal distance, mesophyll surface area, chloroplast placement, and cell wall thickness. Variation in anatomical traits has the potential to influence CO 2 diffusion as well as bundle sheath leakiness to CO 2 , both of which could influence leaf d 13 C.
With our data, we were able to indirectly test the relationship between nutrient uptake and transpiration. If reducing transpiration limits nutrient uptake, transpiration efficiency as a trait for increasing WUE would have limited application. The 19 elements tested here were previously reported to have narrow-sense heritabilities ranging from 0.11 to 0.66 (Baxter 2014). The only element found to be significantly correlated to leaf d 13 C was Molybdenum. Molybdenum is required for several vital biological processes related to nitrogen and water (Baxter 2008). Because molybdenum is a required cofactor for ABA synthesis, maize plants overexpressing molybdenum cofactor sulfurase gene have increased drought tolerance (Lu et al. 2013). However, in our study, we observed a positive correlation between leaf d 13 C and molybdenum, which is contrary to expectation given that an increase in d 13 C signifies a decrease in WUE. Overall, it is encouraging that the majority of elements sampled were not associated with d 13 C. This suggests that breeding for leaf d 13 C as a means to reduce transpiration would be unlikely to result in plants with nutrient uptake deficiencies.
Although the main focus of this study was to investigate leaf d 13 C and its relationship to SLA and nutrient accumulation, the QTL mapping of the analyzed elements was an interesting biproduct. Mapping the leaf ionome of the four RIL families resulted in many significant QTL, including some overlapping intervals for different elements. Multi element QTLs are common, and are thought to be due to loci affecting processes such as the acidification of the rhizosphere or altering the permeability of the casparian strip (Baxter 2015). Interestingly there was not much overlap between the ionomic QTL identified here and a previous study on kernels (Table 4; Baxter et al. 2014). The only overlapping QTL was for rubidium on chromosome 3. There are several possible causes for the limited overlap between these methods. There could be differences between the leaf and grain ionome due to differential mobilization of nutrients from vegetative tissues into kernels during grain fill. In addition, the ionome is strongly influenced by genotype by environment interactions, with many of the QTL identified in previous studies being location specific (Asaro et al. 2016). Environmental interactions may also explain the unexpected high correlation between Fe and Al. One possible mechanism for this correlation is the variance in soil pH, either due to field effects or to the acidification of the rhizosphere by the plants roots. Both Fe and Al are more available at lower pHs. Interestingly, the Fe QTL did overlap a Zinc and Iron transporter (ZmZIP5; Zm00001d036965) on chromosome 6 (Li et al. 2019b).

Figure 3 Element Single Family
Stepwise Regression QTL Mapping. QTL mapping identified 28 QTL across 12 different elements. Significant QTL (alpha ¼ 0.05) for each element are plotted. QTL location is shown across the 10 maize chromosomes (cM) on the x-axis. Dashes indicate a significant QTL, with the NAM RIL family in which the QTL was found designated by color; CML103 (black), CML333 (orange), Tx303 (gray), NC358 (blue). All dashes are the same length for visibility.