bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1 Genome-wide association study for maize leaf cuticular conductance identifies 2 candidate genes involved in the regulation of cuticle development 3 Meng Lin,* Susanne Matschi,† Miguel Vasquez,† James Chamness,*,1 Nicholas Kaczmar,* 4 Matheus Baseggio,*,2 Michael Miller,*,3 Ethan L. Stewart,*,4 Pengfei Qiao,‡ Michael J. 5 Scanlon,‡ Isabel Molina,§ Laurie G. Smith,†,5 and Michael A. Gore*,5 6 *Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell 7 University, Ithaca, NY 14853, USA 8 † 9 CA 92093, USA Section of Cell and Developmental Biology, University of California San Diego, La Jolla, 10 ‡ 11 14853, USA 12 § Department of Biology, Algoma University, Sault Ste Marie, ON P6A 2G4, Canada 13 1 Present address: Department of Genetics, Cell Biology, and Development, University of 14 Minnesota, Saint Paul, MN 55108, USA 15 2 Present address: Seneca Foods Corporation, 1201 N. 4th Street, LeSueur, MN 56058, USA 16 3 Present address: Plant Biology Section, School of Integrative Plant Science, Cornell 17 University, Ithaca, NY 14853, USA 18 4 Vienna BioCenter Core Facilities, Dr. Bohr-Gasse 3, 1030 Vienna, Austria 19 5 Corresponding authors: Michael A. Gore, 358 Plant Science Building, Cornell University, 20 Ithaca, NY 14853, USA E-mail: mag87@cornell.edu 21 Laurie G. Smith, University of California San Diego, Biological Sciences #0116, La Jolla, 22 CA 92093-0116, USA E-mail: lgsmith@ucsd.edu Plant Biology Section, School of Integrative Plant Science, Cornell University, Ithaca, NY 23 1 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 24 Running title 25 GWAS of maize leaf cuticular conductance 26 Article summary 27 The cuticle serves as the major barrier to water loss when stomata are closed at night and 28 under water-limited conditions and potentially relevant to drought tolerance in crops. We 29 performed a genome-wide association study to elucidate the genetic architecture of natural 30 variation for maize leaf cuticular conductance. We identified epidermally expressed candidate 31 genes that are potentially involved in cuticle biosynthesis, trafficking and deposition, cutin 32 polymerization, and cell wall modification. Finally, we observed moderately high predictive 33 abilities for whole-genome prediction of leaf cuticular conductance. Collectively, these 34 findings may help breeders more effectively develop drought-tolerant maize. 35 Keywords 36 Cuticle, Cuticular conductance, Genome-wide association study, RNA sequencing, Whole- 37 genome prediction 2 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 38 Abstract 39 The cuticle, a hydrophobic layer of cutin and waxes synthesized by plant epidermal cells, is 40 the major barrier to water loss when stomata are closed at night and under water-limited 41 conditions. Elucidating the genetic architecture of natural variation for leaf cuticular 42 conductance (gc) is important for identifying genes relevant to improving crop productivity in 43 drought-prone environments. To this end, we conducted a genome-wide association study of 44 gc of adult leaves in a maize inbred association panel that was evaluated in four environments 45 (Maricopa, AZ, and San Diego, CA in 2016 and 2017). Five genomic regions significantly 46 associated with gc were resolved to seven plausible candidate genes (ISTL1, two SEC14 47 homologs, cyclase-associated protein, a CER7 homolog, GDSL lipase, and β-D- 48 XYLOSIDASE 4). These candidates are potentially involved in cuticle biosynthesis, 49 trafficking and deposition of cuticle lipids, cutin polymerization, and cell wall modification. 50 Laser microdissection RNA sequencing revealed that all these candidate genes, with the 51 exception of the CER7 homolog, were expressed in the zone of the expanding adult maize 52 leaf where cuticle maturation occurs. With direct application to genetic improvement, 53 moderately high average predictive abilities were observed for whole-genome prediction of 54 gc in locations (0.46 and 0.45) and across all environments (0.52). The findings of this study 55 provide novel insights into the genetic control of gc and have the potential to help breeders 56 more effectively develop drought-tolerant maize for target environments. 57 3 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 58 59 Introduction The cuticle is a hydrophobic layer covering the surface of the shoot, which limits 60 transpiration, in addition to protecting shoot tissues from UV radiation, heat, and pathogen 61 attack (Shepherd and Wynne Griffiths 2006; Xue et al. 2017). Cuticles have two major lipid 62 components: cutin and waxes. Cutin is a polymer composed of fatty acid derivatives and 63 glycerol, and is highly cross-linked to form an insoluble matrix (Pollard et al. 2008; Fich et 64 al. 2016). Soluble waxes are deposited into and on top of this matrix and consist mainly of 65 aliphatic compounds derived from very-long-chain fatty acids including alkanes, aldehydes, 66 alcohols, ketones, and wax esters (Yeats and Rose 2013). The major pathways for both wax 67 and cutin monomer biosynthesis have been elucidated via genetic and biochemical studies 68 conducted mainly on model plant systems (Yeats and Rose 2013; Lee and Suh 2015; Fich et 69 al. 2016). Many transcriptional regulators of cuticle biosynthesis have also been identified 70 (Borisjuk et al. 2014). Pathways for delivery of cuticle lipids from the intracellular 71 membranes where they are synthesized to the extracellular environment, and to their final 72 destinations external to the cell wall, have been partially elucidated. Golgi-mediated vesicle 73 trafficking (McFarlane et al. 2014), ATP BINDING CASSETTE TRANSPORTER G 74 (ABCG) family proteins (Pighin et al. 2004; Bird et al. 2007; Panikashvili et al. 2010; 75 Bessire et al. 2011; Chen et al. 2011), and extracellular lipid transfer proteins (DeBono et al. 76 2009; Kim et al. 2012) are required for accumulation of multiple classes of cuticle lipids at 77 the shoot surface, but how these transport processes work together is not well understood. 78 Cuticle composition and structure varies widely among species and tissue types (Jetter 79 et al. 2008), and its permeability to water can vary as much as three orders of magnitude 80 (Kerstiens 2006). Relationships between cuticle composition, structure, and water barrier 81 function are complex. The simple idea that cuticle impermeability to water increases with 82 increased wax load or cuticle thickness has been repeatedly refuted; instead, cuticle 4 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 83 composition and the organization of components appear to determine water barrier function 84 (Riederer and Schreiber 2001). Cuticle composition is also modulated by environmental 85 factors. Light, temperature, and osmotic stress influence both the quantity and composition of 86 cuticular waxes (Shepherd and Wynne Griffiths 2006; Kosma and Jenks 2007). While high 87 relative humidity usually suppresses wax production (Sutter and Langhans 1982; Koch et al. 88 2006), drought stress has been found to change cuticular wax composition (Panikashvili et al. 89 2007; Kosma et al. 2009) and increase wax deposition in several crop species (Shepherd and 90 Wynne Griffiths 2006; Cameron et al. 2006; Kosma and Jenks 2007). 91 In addition to modulation of its composition by drought stress, a variety of other 92 observations support a role for the cuticle in drought tolerance. Most mutations and other 93 genetic modifications affecting cuticle composition or wax load increase its permeability to 94 water, and this has often been associated with decreased drought tolerance (Zhou et al. 2013; 95 Zhu and Xiong 2013; Li et al. 2019). In a few cases, overexpression of cuticle lipid 96 biosynthetic enzymes, or their transcriptional regulators, increased cuticular lipid abundance 97 and drought tolerance (Aharoni et al. 2004; Zhang et al. 2005; Bourdenx et al. 2011; Wang et 98 al. 2012). Glaucousness (a visible trait resulting from abundant accumulation of epicuticular 99 wax crystals on shoot tissue surfaces) was selected for during the domestication of wheat and 100 involves genes regulating wax biosynthesis (Hen-Avivi et al. 2016; Bi et al. 2016); analyses 101 of heritable variation in glaucousness in wheat and barley has revealed positive correlations 102 between this trait and drought tolerance (Febrero et al. 1998; Guo et al. 2016). These findings 103 point to the potential relevance of cuticle modification for increasing drought tolerance in 104 cereal crops. 105 There is a longstanding interest in breeding for cuticle-related traits such as drought 106 tolerance, but the complexity of the relationship between cuticle composition and resistance 107 to dehydration, and the lack of methods amenable to high-throughput phenotyping for 5 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 108 cuticle-associated traits, has hampered progress in this area (Petit et al. 2017). Genome-wide 109 association studies (GWAS) (Yu et al. 2006; Zhang et al. 2010; Lipka et al. 2015), taking 110 advantage of natural variation in cuticle permeability and historical recombination events 111 within a single species, offer an attractive approach to identifying genes and alleles that can 112 decrease cuticle permeability. Additionally, findings from GWAS of cuticle-related 113 phenotypes could be used to better inform the application of genomic selection strategies 114 (Meuwissen et al. 2001; Lorenz et al. 2011) for increasing drought tolerance in crop species. 115 In this study, we utilized GWAS combined with laser microdissection RNA 116 sequencing (LM-RNAseq) of epidermal tissue samples, from zones of the expanding maize 117 leaf where the cuticle develops, to implicate candidate genes controlling the water barrier 118 function of the maize leaf cuticle. Methods directly measuring permeability of the cuticle to 119 water (Valeska Zeisler-Diehl et al. 2017) are not amenable to the high-throughput 120 phenotyping needed for GWAS involving analysis of thousands of samples. Thus, we utilized 121 a phenotyping strategy that indirectly measures cuticle water barrier function by calculating 122 the drying rates of detached leaves placed in the dark to close stomata (Ristic and Jenks 123 2002). While some water may be lost via stomata that are not completely sealed, the sealing 124 of stomata is thought to depend on cuticular flaps lining the stomatal pore that overlap when 125 stomata are closed (Zhao and Sack 1999; Kosma and Jenks 2007). Thus, we used this 126 phenotyping strategy with the expectation that it primarily measures cuticle-dependent water 127 loss, which we refer to as adult leaf cuticular conductance (gc). Maize is most sensitive to 128 drought stress at flowering (Grant et al. 1989), when juvenile leaves have already senesced 129 and only adult leaves remain. Thus, we utilized adult leaves of field-grown plants for this 130 analysis. Moreover, since cuticle composition is known to vary depending on the growth 131 environment (Shepherd and Wynne Griffiths 2006), we utilized plants grown in two 132 contrasting environments: the cooler and more humid environment of San Diego, CA and the 6 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 133 hotter and more arid environment of Maricopa, AZ. We aimed to (i) detect genomic regions 134 associated with natural variation in gc; (ii) enhance selection of candidate genes impacting 135 this trait through an LM-RNAseq analysis of the expanding maize adult leaf; and (ⅲ) 136 evaluate whole-genome prediction models to facilitate genomic selection on this trait, 137 providing tools for potential improvement of drought-tolerance in maize for target 138 environments. 139 Materials and Methods 140 Plant materials and experimental design 141 A set of 468 maize inbred lines from the Wisconsin Diversity panel (Hansey et al. 142 2011) was evaluated for adult leaf cuticular conductance (gc). The inbred lines were planted at 143 the Maricopa Agricultural Center, Maricopa, AZ, and University of California San Diego, 144 San Diego, CA, in 2016 and 2017. The layout of the experiment in each of the four 145 environments (location × year combination) was arranged as an 18 x 26 incomplete block 146 design. Each incomplete block of 18 experimental lines was augmented by the random 147 positioning of two checks (B73 and Mo17, 2016; N28Ht and Mo17, 2017). The entire 148 experiment of 468 unique inbred lines and repeated checks was grown as a single replicate in 149 each environment. Edge effects were reduced by planting a locally adapted maize inbred line 150 around the perimeter of the experiment. Experimental units were one-row plots of 3.05 m 151 (Maricopa) and 4.88 m (San Diego) in length with 1.016 m inter-row spacing. There was a 152 0.91 m alley at the end of each plot. In each plot, 12 kernels were planted, followed by 153 thinning the stand as needed to limit plant-to-plant competition. We obtained and summarized 154 meteorological data (Table S1) from automated weather stations located at a distance of ~200 155 m or less from the experimental field sites in Maricopa, AZ (Arizona Meteorological 156 Network; http://ag.arizona.edu/azmet/index.html) and San Diego, CA (Scripps Hydroclimate 157 Weather Station). 7 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 158 159 gc evaluation and phenotypic data analysis Initially, a porometer (model SG-1, Decagon Devices, Pullman, WA) was used to 160 measure diffusion conductance for adult leaves collected from the field in the morning (at the 161 same time of day that leaves were harvested for gc analysis, within a week of their anthesis 162 dates), and transferred to a dark room where they were kept with cut leaf bases immersed in 163 water throughout the analysis. Diffusion conductance was measured at a series of timepoints 164 for a selection of lines that had previously been demonstrated to have high gc, and B73 as a 165 reference standard harvested each sampling day. As shown in Figure S1, most lines reached 166 minimum conductance values, indicating stomatal closure, within 30 min, and all lines 167 reached minimum values within 90 min. Based on these findings, we concluded that 2 h in the 168 dark was sufficient to achieve stomatal closure, even for lines with high gc. 169 The findings from the porometer study were used to inform our phenotyping of the 170 Wisconsin Diversity panel. To inform sampling times and help control for maturity 171 differences, flowering time (days to anthesis, DTA) was recorded as the number of days from 172 planting to initiation of pollen shed for 50% of the plants in a plot. For each plot, gc was 173 evaluated on an adult leaf from five plants when 50% of the plants in the plot were at 174 anthesis. From each of five plants, the leaf subtending the uppermost ear, or the leaf 175 immediately above or below it, was excised at 2.5 cm below the ligule. If there were fewer 176 than five plants to measure in a plot, two, three, or four plants were evaluated to represent the 177 gc of that plot. Cut ends were immersed in water and leaves were incubated in a dark, well- 178 ventilated room for 2 h at 20-22 °C and 55-65% RH to close stomata while ensuring that the 179 leaves were fully hydrated. Subsequently, each leaf was hung on a line from a boot clip in the 180 same dark, temperature- and humidity-controlled room, after carefully wiping the leaves to 181 remove excess water. In 2017, the cut leaf end was wrapped in parafilm when hung to further 182 reduce the incidence of non-cuticular transpiration. The wet weight of each leaf was then 8 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 183 recorded every 45 min over a time period of 3 h, for a total of five measurements per leaf. 184 Leaves were subsequently dried at 60 °C for 4 d in a forced-air oven, and afterward weighed. 185 To investigate the degree to which dry weight approximates surface area, we 186 constructed an imaging table that consisted of a digital camera (Canon PowerShot S110) 187 mounted at a fixed distance and folding table top to help flatten the adult maize leaves. In 188 2017, each experimental leaf that had completed the gc phenotyping process was flattened, 189 imaged, and then oven-dried as earlier described and weighed. Each image contained 4 to 7 190 leaves, depending on their sizes. Next, images were individually analyzed with a custom 191 ImageJ (Schneider et al. 2012) macro to estimate the surface area (mm2) of each leaf. Even 192 though every attempt was made to flatten the leaves, there were frequent occurrences of 193 where wavy leaves could not be completely flattened. Very strong Pearson’s correlations 194 were observed between the plot-level averages of leaf dry weight and surface area in MA (r = 195 0.93) and SD (r = 0.92) in 2017 (Figure S2; Table S2). Given these strong correlations and 196 challenges in measuring surface area, leaf dry weight was considered to be a reasonable 197 approximation of leaf surface area, and was used in the calculation of cuticular conductance 198 as described below. 199 200 Adult leaf cuticular conductance (gc) from unit surface area was calculated as follows: gc (g∙h-1∙g-1) = - b / dry weight, 201 where b (g∙h-1) is the coefficient of the linear regression of leaf wet weight (g) on time (h), 202 and dry weight (g) is an approximation of leaf surface area. 203 To screen the phenotypic data (flowering time or gc) for significant outliers, univariate 204 mixed linear models were fitted as follows: (1) each single environment; (2) the two 205 environments in each location; and (3) all four environments. The model terms included 206 grand mean and check as fixed effects and environment, genotype, genotype-by-environment 207 (G×E) interaction (only for Models 2 and 3), incomplete block within environment, and 9 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 208 column within environment as random effects. The Studentized deleted residuals (Neter et al. 209 1996) generated from these mixed linear models were assessed and significant (α = 0.05) 210 outliers removed. For each outlier screened phenotype, an iterative mixed linear model fitting 211 procedure was conducted for each of the three (1 to 3) full models in ASReml-R version 3.0 212 (Gilmour et al. 2009). All random terms that were not significant at α = 0.05 in a likelihood 213 ratio test (Littell et al. 2006) were removed from the model, allowing a final best-fit model to 214 be obtained for each phenotype. These final models (Table S3) were used to generate a best 215 linear unbiased predictor (BLUP) for gc and DTA for each line (Table S4). 216 Variance component estimates from the fitted full models (Table S5) were used to 217 estimate heritability on a line-mean basis (Holland et al. 2010; Hung et al. 2012) for each 218 phenotype in a location (Model 2) and across all four environments (Model 3). Standard 219 errors of the heritability estimates were calculated with the delta method (Lynch et al. 1998; 220 Holland et al. 2010). Pearson’s correlation coefficients were used to evaluate the strength of 221 correlation between BLUP values for gc and DTA in all pairwise combinations. The 222 significance of correlations (α = 0.05) was tested using the function “cor.test” in R version 223 3.5.1 (R core team, 2018). 224 DNA extraction and genotyping 225 For each inbred line, total genomic DNA was isolated from a leaf tissue sample 226 consisting of bulked young leaves from a single plant. The leaf tissue samples were 227 lyophilized and ground using a GenoGrinder (Spex SamplePrep, Metuchen, NJ, USA), 228 followed by isolation of total genomic DNA using the DNeasy 96 Plant Kit (Qiagen 229 Incorporated, Valencia, CA, USA). The genotyping-by-sequencing (GBS) procedure was 230 conducted on the DNA samples following Elshire et al. (2011) with ApeKI as the restriction 231 enzyme at the Cornell Biotechnology Resource Center (Cornell University, Ithaca, NY, 10 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 232 USA). The constructed 192-plex GBS libraries were sequenced on a NextSeq 500 (Illumina 233 Incorporated, San Diego, CA, USA). 234 We called genotypes at 955,690 high-confidence single-nucleotide polymorphism 235 (SNP) loci in B73 RefGen_v2 coordinates following the method of Baseggio et al. (2019). 236 Briefly, the raw SNP genotype calls were initially filtered to exclude singleton and doubleton 237 SNPs (a minor allele observed in only a single line) and retain only biallelic SNPs having a 238 SNP call rate greater than 10% and minimum inbreeding coefficient of 0.8 per Romay et al. 239 (2013). Additionally, only lines with a call rate greater than 40% were retained. Missing SNP 240 genotypes were partially imputed using FILLIN (Swarts et al. 2014) with a set of maize 241 haplotype donors files having a 4 kb window size 242 (AllZeaGBSv2.7impV5_AnonDonors4k.tar.gz, available at panzea.org). Given that missing 243 genotype data still remained following partial imputation, the imputed genotype data set was 244 further filtered to retain SNPs with a minimum call rate of 0.6, a minimum inbreeding 245 coefficient of 0.8, and a minimum minor allele frequency of 5%. The final complete set 246 contained 235,004 high-quality SNP markers scored on 451 lines that had a BLUP value for 247 gc in one or more environments. The genome coordinates of the SNP loci were uplifted by 248 aligning 101 bp context sequences containing target SNPs (± 50 bp) to the B73 249 RefGen_AGPv4 reference genome with Vmatch (Kurtz 2003). 250 Population structure analysis 251 We estimated population structure from the 451 line × 235,004 SNP genotype matrix 252 in fastSTRUCTURE version 1.0 (Raj et al. 2014), with the number of ancestral populations 253 (K) varying from 1 to 10 using the simple prior. To complement population structure 254 inference with fastSTRUCTURE, a principal component analysis (PCA) was performed on 255 the identical genotype matrix with the prcomp function in R version 3.5.1 (R core team, 256 2018). Next, results from fastSTRUCTURE and PCA were jointly interpreted and 11 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 257 complemented with pedigree information (Hansey et al. 2011) to choose K=6 as the number 258 of subpopulations. Additionally, lines were classified following group assignments of Hansey 259 et al. (2011) (Figure S3A). In fastSTRUCTURE, the simple prior approach was used with 260 K=6 to calculate the subpopulation composition of each line (Figure S3B). Lines with an 261 assignment value of Q ≥ 0.50 were assigned to subpopulations (1, 2, 3, 4, 5, or 6), while those 262 with assignment values of Q < 0.50 for all six subpopulations were classified as admixed 263 (Figure S3C; Table S6). This information was used to inform a stratified sampling approach 264 used for whole-genome prediction (WGP). 265 Genome-wide association study 266 To identify associations between each of the 235,004 SNP markers and the BLUP 267 values of gc from the 451 inbreds, a GWAS was conducted with a univariate mixed linear 268 model that employed population parameters previously determined (P3D) approximation 269 (Zhang et al. 2010) in the R package Genome Association and Prediction Integrated Tool 270 (GAPIT) version 3.0 (Lipka et al. 2012). To control for phenotypic variation due to maturity, 271 population structure, and unequal relatedness, the fitted mixed linear models included BLUPs 272 of flowering time (DTA) of corresponding environments, principal components (PCs), and a 273 genomic relationship (kinship) matrix. The kinship matrix based on the centered IBS method 274 (Endelman and Jannink 2012) implemented in TASSEL 5.0 (Bradbury et al. 2007) was 275 calculated from a subset of 41,259 SNPs remaining after linkage disequilibrium (LD) pruning 276 (r2 ≤ 0.2) of the complete marker data set in PLINK version 1.09_beta5 (Purcell et al. 2007). 277 The optimal number of covariates (BLUPs of DTA and PCs based on the genotypic matrix) to 278 include in the mixed linear model was determined with the Bayesian information criterion 279 (BIC) (Schwarz 1978). Prior to conducting the GWAS, the remaining missing genotypes for 280 all SNP loci were conservatively imputed as a heterozygote in GAPIT. A likelihood-ratio- 281 based R2 statistic (R2LR) (Sun et al. 2010) was used to estimate the amount of phenotypic 12 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 282 variation explained by the model with or without a tested SNP. The false discovery rate 283 (FDR) was controlled at 10% with the Benjamini-Hochberg procedure (Benjamini and 284 Hochberg 1995). 285 Linkage disequilibrium 286 The squared allele-frequency correlation (r2) method (Hill and Weir 1988) was used 287 to estimate linkage disequilibrium between all pairs of SNPs on a chromosome in TASSEL 288 5.0 (Bradbury et al. 2007). The complete set of 235,004 high-quality SNPs was used for LD 289 estimation; however, the missing SNP genotypes remaining after FILLIN imputation were not 290 imputed with a heterozygote value prior to LD analysis. The estimates of LD were used to 291 approximate the physical distance at which LD decayed to genome-wide background levels 292 (r2 < 0.10) and assess the pattern of LD surrounding GWAS-detected SNPs. 293 Candidate gene identification 294 We searched for candidate genes within ±200 kb of the most significant SNP (peak 295 SNP) for each detected genomic region associated with gc. Gene functional annotations 296 (B73_RefGen_v4) of identified candidate genes were obtained from Gramene 297 (http://www.gramene.org/). Additionally, the Gramene Mart tool was used to identify 298 homologs of candidate genes in Arabidopsis thaliana (L.) Heynh. (Columbia-0 ecotype) and 299 rice (Oryza sativa L. ssp. Japonica cv. ‘Nipponbare’). Gene functional annotations for the 300 Arabidopsis and rice homologs were obtained from TAIR (https://www.arabidopsis.org/) and 301 RAP-DB (https://rapdb.dna.affrc.go.jp/), respectively. 302 Transcript abundance analysis for candidate genes 303 To investigate the patterns of transcript abundance for candidate genes potentially 304 regulating gc of maize leaves, we analyzed the LM-RNAseq dataset of Qiao et al. (2019) 305 generated from epidermal and internal tissues that were LM from seven 2-cm-long intervals 306 (six intervals from 2-14 cm, and one interval from 20-22 cm) from the expanding leaf 8 of 13 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 307 maize inbred B73 with three biological replicates and three plants per replicate. The 308 generated RNAseq reads were aligned to B73 genome RefGen_v4 with HiSAT2 (Kim et al. 309 2015) and counted with HTSeq (Anders et al. 2015). Transcript abundance levels for 310 canonical protein coding genes were characterized as counts per million (cpm) after 311 normalizing against library sizes with the R package edgeR 3.3.2 (Robinson et al. 2010). To 312 identify potential outlier samples, genes with >1 cpm in two or more samples were retained. 313 Next, a PCA was conducted on the cpm values of all the retained genes on all samples 314 sequenced (Figure S4), and the epidermal and internal samples from interval 4-6 cm in the 3rd 315 replicate were removed as outliers (circled samples in Figure S4). With the remaining 316 epidermis samples, a gene was declared expressed in the epidermis if its cpm was greater 317 than one for at least three out of the 20 samples. 318 Whole-genome prediction 319 WGP of gc was conducted using a maximum likelihood approach to ridge regression 320 in the ‘rrBLUP’ R package (Endelman 2011). The final complete set of 235,004 SNPs with 321 post-FILLIN missing genotype data imputed as a heterozygote was used for WGP. We 322 trained WGP models with gc BLUP values from either MA (two environments), SD (two 323 environments), or AllEnv (four environments), then each of the three models were separately 324 used to predict gc for MA, SD, and AllEnv. Collectively, this resulted in a total of nine 325 prediction scenarios. Five-fold cross-validation as described in Owens et al. (2014) was 326 performed to evaluate the predictive ability for gc in each of the nine scenarios by calculating 327 the Pearson’s correlation between observed BLUP values and genomic estimated breeding 328 values. We used a stratified sampling approach to account for the presence of population 329 structure, enabling each fold to have the proportion of subpopulations (1 to 6 and admixed) 330 observed for the whole panel (Table S6). The procedure was conducted 50 times for each 331 scenario, and the predictive ability was calculated as a mean of the correlations. 14 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 332 333 Data availability The raw GBS sequence data were deposited in the National Center of Biotechnology 334 Information (NCBI) Sequence Read Archive (SRA) under accession number SRP160407 and 335 in BioProject under accession PRJNA489924. The BLUP values of the phenotypes are 336 provided in Table S4. The FILLIN partially imputed SNP genotype data for the 235,004 loci 337 scored on the 451 inbred lines, weather data from Maricopa and San Diego in 2016 and 2017, 338 and leaf image data from Maricopa and San Diego in 2017 are available at CyVerse 339 (http://datacommons.cyverse.org/browse/iplant/home/shared/GoreLab/dataFromPubs/Lin_Le 340 afCuticle_2019). The ImageJ macro for image analysis is available at Github 341 (https://github.com/GoreLab/Maize_leaf_cuticle/blob/master/GWAS_CE/cliveMacro_2.4.txt) 342 . RNAseq reads of epidermal cells along the developmental gradient of the expanding leaf 8 343 of maize inbred B73 were deposited in the NCBI SRA under accession number SRP116320 344 and in BioProject under accession PRJNA400334, as described in Qiao et al. (2019). 345 Supplemental material available at FigShare: URL to be added 346 Results 347 Phenotypic variability 348 To establish the feasibility of GWAS for gc, the extent of phenotypic variation for 349 adult leaf gc was assessed in the Wisconsin Diversity panel of more than 450 maize inbred 350 lines that was grown in two field locations (SD, San Diego, CA; and MA, Maricopa, AZ) in 351 2016 and 2017. The climatically contrasting environments of the SD and MA locations 352 (Table S1) had the lowest and highest average gc estimates, respectively (Table 1). There was 353 a moderately strong correlation for gc between years in a location (r = 0.52 and 0.54), 354 whereas correlations were slightly weaker (r = 0.36 to 0.48) between locations within and 355 across years (Figure S5). With respect to flowering time, gc tended to have a moderately 356 weak negative correlation with days to anthesis (DTA) in each environment (r = -0.14 to - 15 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 357 0.50), with the exception of a negative correlation of moderate strength (-0.50) in the MA17 358 environment (Figure S5). Comparatively, the first two principal components (PCs) but not the 359 third PC calculated from the SNP genotype matrix (PC1, 3.68%; PC2, 1.83%; and PC3, 360 1.62%), which captures patterns of population structure, showed comparable correlations 361 with gc across all four environments (|r| = 0.14 to 0.31). Heritability on a line-mean basis 362 across (AllEnv) locations was 0.71, with a range of 0.63 to 0.67 in locations (Table 1). These 363 findings suggest that most phenotypic variation for gc is under genetic control. 364 Genome-wide association study 365 The genetic basis of natural variation for gc was examined in the Wisconsin Diversity 366 panel that had been genotyped with 235,004 high-quality SNP markers. We conducted 367 GWAS in (MA or SD) and across locations (AllEnv) using a mixed linear model that controls 368 for population structure and unequal relatedness. In addition, flowering time, which was also 369 found to be moderately heritable (0.88 to 0.93) for MA, SD, and AllEnv (Table S7), was 370 retained as a covariate for GWAS of gc from MA and AllEnv based on the BIC values, thus 371 attenuating the confounding effect of maturity when estimating allelic effects. A total of nine 372 unique SNPs localized to five different genomic regions were significantly associated with gc 373 from MA and/or AllEnv at a genome-wide FDR of 10% (Figures 1 and S6). Of these nine 374 SNPs, two were associated with gc in both MA and AllEnv, for a total of 11 significant 375 marker-trait associations (Table 2). No significant associations were detected for SNP 376 markers with gc from SD at the 10% FDR level (Figure S6). 377 In the Wisconsin Diversity panel, genome-wide LD decayed to background levels (r2 378 < 0.1) by 213 kb at the 90% percentile of the r2 distribution (Figure S7). This percentile 379 cutoff of the distribution was used to sample the large variance in LD structure presumably 380 caused by rare variants (Wallace et al. 2014) and allow for the inclusion of potential distant 381 regulatory elements as has been observed for cloned quantitative trait loci (QTL) in maize 16 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 382 (Salvi et al. 2007; Studer et al. 2011). Therefore, the genomic search space to identify 383 candidate genes was restricted to ± 200 kb window of the five peak SNPs that each tagged a 384 unique genomic interval, resulting in the inclusion of 76 annotated genes (Table S8). Of these 385 76 genes, 57 are canonical protein coding genes, while the other 19 encode non-coding 386 RNAs. 387 The strongest association signal was identified on chromosome 4 (Figure 1), with 388 SNPs S4_30201436 and S4_30201599 comprising the peak associations with gc in MA (P- 389 value 2.68 × 10-7) and AllEnv (P-value 1.04 × 10-6), respectively. The peak SNP 390 S4_30201599 resides within a microRNA locus (zma-MIR2275a) and is located 1,146 bp 391 from a gene (Zm00001d049479) encoding an INCREASED SALT TOLERANCE1-LIKE1 392 (ISTL1) protein (Figure 2; Table S8). ISTL1 is the most probable candidate gene underlying 393 this association signal because it is predicted to function together with the endosomal sorting 394 complex required for transport machinery (ESCRT) and ATPase VACUOLAR PROTEIN 395 SORTING4 (VPS4) in the formation of multivesicular bodies (MVBs) (Hill and Babst 2012; 396 Buono et al. 2016). 397 The genomic region associated with gc on chromosome 1 was detected in MA, with 398 the signal peak defined by SNP S1_27055151534 (P-value 4.37 × 10-7) (Table 2). The peak 399 SNP is located 8,190 bp from the nearest gene (Zm00001d033835), but this gene encodes an 400 uncharacterized protein. Two genes adjacent to this one encode closely related SEC14 401 homologs, separated by about 50 kb (Zm00001d033836, ~50 kb from peak SNP; and 402 Zm00001d033837, ~100 kb from peak SNP) (Figure S8; Table S8). SEC14 proteins in plants 403 and other eukaryotes function in the transfer of phosphoinositides between different cellular 404 membranes (Huang et al. 2016) and are noteworthy candidates for regulators of gc because of 405 their function in stimulating vesicle formation from the trans Golgi network in yeast (de 406 Campos and Schaaf 2017), given the role of Golgi-dependent vesicle trafficking in cuticle 17 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 407 biosynthesis in plants (McFarlane et al. 2014). Additional plausible candidates on 408 chromosome 1 encode a homolog of ECERIFERUM7 (CER7) (Zm00001d033842) and a 409 CYCLASE-ASSOCIATED PROTEIN (CAP) (Zm00001d033830) located (+) 179,198 bp 410 and (-) 152,746 bp from the peak SNP, respectively (Figure S8). CER7, a 3′-to-5′ 411 exoribonuclease family protein, is a core subunit of the RNA degrading exosome that 412 positively regulates the expression of CER3/WAX2/YRE, a key wax biosynthesis gene in 413 Arabidopsis (Kurata et al. 2003; Chen et al. 2003). CAP, a regulator of actin cytoskeleton 414 dynamics, is highly conserved in plants, yeast, flies, and mammals (Balcer et al. 2003; Ono 415 2013). 416 An additional association with gc was detected in MA on chromosome 8, defined by 417 S8_152967952 as the peak SNP (P-value 1.92 x 10-6) along with two neighboring SNPs 418 (S8_152967929 and S8_152967966). The closest gene (Zm00001d011655) is located 4,067 419 bp from the peak SNP, and it encodes a MITOGEN-ACTIVATED PROTEIN KINASE 420 KINASE KINASE18. A gene encoding a GLY-ASP-SER-LEU ESTERASE/LIPASE (GDSL 421 lipase) (Zm00001d011661) was identified 170,643 bp from the peak SNP (Figure S9; Table 422 S8), and it was identified as the most plausible candidate in this region based on the known 423 functions of GDSL lipases in lipid biosynthesis (Girard et al. 2012; Yeats et al. 2012b). 424 Finally, in AllEnv, gc was also significantly associated with genomic regions on 425 chromosomes 7 and 10 with peak SNPs S2_231152916 (P-value 7.88 x 10-7) and 426 S10_144686998 (P-value 1.08 x 10-6), respectively. The peak SNP on chromosome 7 is 427 located within a gene encoding CYTOKININ-N-GLUCOSYLTRANSFERASE1 428 (Zm00001d019250) (Figure S10; Table S8). The peak SNP on chromosome 10 is located 429 within a gene that encodes β-D-XYLOSIDASE 4 (Zm00001d026415), a cell wall 430 biosynthetic enzyme (Figure S11; Table S8). The latter is a plausible candidate given that the 431 outer epidermal cell wall and cuticle are closely connected and polysaccharides can be found 18 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 432 embedded in some parts of the cuticle (Yeats and Rose 2013). Therefore, sugar- or cell wall- 433 modifying enzymes might directly or indirectly influence cuticle composition and 434 permeability (Müller et al. 2013). 435 Transcript abundance analysis of identified candidate genes 436 We analyzed the transcript abundance of the 57 candidate genes that encode a protein 437 to help further prioritize which of them may be involved in the genetic control of gc. Given 438 that the leaf cuticle is produced by epidermal cells, we analyzed LM-RNAseq data of 439 epidermal cells from seven 2-cm intervals excised from the expanding leaf 8 of maize inbred 440 B73 (Qiao et al. 2019), representing sequential stages in cuticle maturation (Bourgault et al. 441 2019). Of the 57 candidate genes, 24 were found to be expressed in at least one of the seven 442 sampled intervals (Figure S12). 443 Of the seven promising candidate genes contributing to the genetic regulation of 444 natural variation for gc through potential influence on cuticle development, those encoding 445 CAP, ISTL1 protein, GDSL lipase, β-D-XYLOSIDASE 4, and the two SEC14-like proteins 446 were found to be epidermally expressed (Figure 3), while the homolog of CER7 was not 447 declared to be expressed in the leaf epidermis based on its transcript abundance (Table S9). In 448 general, transcript abundance increased for the ISTL1 protein and GDSL lipase along the leaf 449 developmental gradient from the base (youngest, 2-4 cm) towards the tip (oldest, 20-22 cm) 450 of the leaf, with the exception of a reduced transcript abundance for GDSL lipase in the 20- 451 22 cm interval, at which time cuticle maturation is already completed (Bourgault et al. 2019). 452 In contrast, CAP showed the opposite trend, with decreases in transcript abundance across the 453 gradient from 2-4 cm to 20-22 cm intervals. The two genes encoding the SEC14-like 454 proteins, Zm00001d033836 and Zm00001d033837, had relatively constant transcript 455 abundances from the base to the tip, with a slightly increased abundance at the 20-22 cm 456 interval for Zm00001d033837. The gene encoding β-D-XYLOSIDASE 4 had a peak 19 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 457 abundance at the 6-8 cm interval, followed by a progressive decrease in abundance towards 458 the leaf tip. These observed transcript abundance patterns are compatible with roles in cuticle 459 maturation. 460 Whole-genome prediction of gc 461 Whole-genome prediction (WGP) was performed with the complete marker data set 462 of 235,004 SNPs to assess the feasibility of implementing genomic selection for the genetic 463 improvement of gc in maize breeding populations. Predictive ability was evaluated on gc from 464 MA, SD, and AllEnv in a scheme that used each of the three as a training set. The overall 465 predictive ability was 0.48 across all of the nine training-prediction combinations, with 466 individual predictive abilities ranging from 0.40 to 0.54 (Figure 4; Table S10). Not 467 unexpectedly, the highest average predictive ability (0.52) was achieved when AllEnv was 468 the training set, followed by MA (0.46) and SD (0.45). Indicative of contrasting 469 environmental conditions between field locations, MA and SD showed lower predictive 470 abilities for each other compared to when MA or SD was used to predict itself. However, 471 AllEnv as the training set resulted in more similar prediction abilities for gc in MA (0.53) and 472 SD (0.49). 473 Discussion 474 The cuticle, a crucial barrier to plant water loss (Shepherd and Wynne Griffiths 2006; 475 Xue et al. 2017), has been widely studied using single gene loss- and gain-of-function 476 mutants, and interspecies crosses (Kosma and Jenks 2007; Yeats et al. 2012a; Beisson et al. 477 2012; Yeats and Rose 2013). However, no prior study has explored natural variation in a 478 single species to elucidate the genetic control of cuticle function as a water barrier on a 479 genome-wide scale. In this study, we focused on leaf cuticular conductance (gc), a trait of 480 potential agronomic value related to drought tolerance. To that end, we evaluated the 481 phenotypic variation of gc using detached adult leaves in a maize diversity panel in two 20 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 482 climatically contrasting field locations (Table S1) and conducted a GWAS and WGP to 483 explore the genetic architecture and the potential genetic gains that could be expected under 484 selection in a breeding program for gc. 485 Analysis was performed to evaluate the repeatability of gc in and across locations, as 486 well as how gc related to flowering time (DTA). Moderately strong correlations (r) of gc 487 measurements were observed between years in a location, while relatively weaker 488 correlations existed between locations (Figure 5). Given the lower relative humidity and 489 higher solar radiation and air temperature of MA compared to SD during maize growing 490 seasons (Table S1), the weaker correlative patterns between locations imply that gc is 491 moderately responsive to environmental conditions. Noticeably, weak to moderately strong 492 negative correlations (MA16, -0.33; MA17, -0.50; SD16, -0.14; SD17, -0.20) existed 493 between gc and flowering time collected from the same year in each location. These findings 494 suggest low to moderate levels of confounding between gc and maturity depending on the 495 environment, with the prolonged exposure of the developing cuticle of plants to more 496 extreme weather variables possibly explaining the stronger correlations in the MA field 497 location. Even though these results implicated the need to assess for the confounding effect of 498 flowering time when conducting GWAS, the ~2-fold range in gc variation, calculated as the 499 ratio between the maximum and minimum BLUP values, in (MA and SD) and across 500 (AllEnv) locations, and the moderately high heritabilities of gc (Table 1) indicate that this 501 phenotype should be very amenable to genetic dissection and prediction in maize. 502 We conducted a GWAS that resulted in the identification of five genomic regions 503 associated with gc in the Wisconsin Diversity panel. Of these five genomic regions, two were 504 detected only in MA (chromosomes 1 and 8), two in AllEnv (chromosomes 7 and 10), and 505 one (chromosome 4) was detected in both MA and AllEnv (Figure 1; Table 2). Notably, no 506 significant associations at the 10% FDR level were detected only in SD, implying that G×E 21 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 507 exists for gc as shown by the between-location correlation analysis (Figure S5). This is also 508 further supported by comparing the ratio of genetic variance to G×E variance, which was 509 3.95, 2.57, and 1.81 for SD, MA, and AllEnv, respectively (Table S5). All together, these 510 observations are consistent with those of previous studies showing that both the content and 511 composition of cuticular waxes can be influenced by environmental factors such as ultraviolet 512 (UV) radiation, temperature, humidity, and drought stress (Shepherd and Wynne Griffiths 513 2006; Xue et al. 2017). It has been reported that drought stress and low humidity increased 514 wax loads in many plant species (Sutter and Langhans 1982; Shepherd and Wynne Griffiths 515 2006; Koch et al. 2006; Kosma and Jenks 2007), but a different combination of UV intensity 516 and temperature resulted in more complex patterns of wax products (Baker 1974; Riederer 517 and Schneider 1990; Shepherd et al. 1997). Although it is not straightforward to explain how 518 each environmental factor affected gc through wax biosynthesis and deposition, the 519 differences in GWAS results in MA and SD suggests that other environment-specific factors 520 in MA besides the high temperature, elevated UV radiation, and low humidity (Supplemental 521 Table S1) might be inducing compositional features that affect water loss relative to plants in 522 SD. 523 As the major barrier of water vapor diffusion from the intercellular air space to the 524 atmosphere when stomata are closed, the composition and structure of the cuticle can affect 525 water evaporation resistance, thus causing variation in gc. The candidate genes we propose 526 for gc via GWAS encode proteins with functions that can be related to cuticle formation. 527 Three of the candidate genes implicated by our GWAS in control of gc encode 528 proteins with likely functions in membrane trafficking, a process needed for deposition of 529 cuticle lipids at the cell surface as well as for delivery of lipid transport proteins to the plasma 530 membrane (McFarlane et al. 2014). Two adjacent genes, likely resulting from a recent 531 tandem duplication, encode closely related SEC14 homologs. These belong to a subset of 22 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 532 SEC14 homologs in maize lacking the GOLD (Golgi dynamics) and nodulin domains present 533 in SEC14s that have been functionally characterized in plants (Huang et al. 2016; de Campos 534 and Schaaf 2017). Nevertheless, SEC14 domains are well known to function in movement of 535 phosphoinositides (signaling lipids) between cellular membranes; in yeast this function 536 stimulates vesicle formation from the trans-Golgi network (de Campos and Schaaf 2017). 537 Thus, the maize SEC14 proteins we identified may function in vesicle trafficking as well. 538 Another membrane-trafficking-related candidate encodes an IST1 family protein 539 closely related to Arabidopsis ISTL1 (Buono et al. 2016). Yeast IST1 functions with the 540 ATPase VPS4 to regulate the assembly and function of ESCRTIII complexes, which 541 facilitate the formation of vesicles inside MVBs (Hill and Babst 2012). MVBs are late 542 endosomal compartments that function in degradation of proteins retrieved via endocytosis 543 from the plasma membrane. A function for ISTL1 proteins in degradation of plasma 544 membrane proteins and formation of MVBs in Arabidopsis has been demonstrated (Buono et 545 al. 2016). Thus, maize ISTL1 could function in an endosomal pathway critical for 546 maintenance of plasma membrane-associated proteins that deliver cuticle lipids to the 547 extracellular environment such as ABCG transporters and lipid transfer proteins. Moreover, 548 MVB-like organelles have been implicated as the source of extracellular vesicles that deliver 549 certain proteins and small RNAs to the extracellular environment in plants (Rutter and Innes 550 2018). This suggests the additional possibility that maize ISTL1 could impact gc because of a 551 role in the formation of MVB vesicles containing cuticle lipids or GDSL lipases required for 552 cutin polymerization, which may be released into the extracellular environment via fusion of 553 MVBs with the plasma membrane. Confirmation of a role for MVBs in cuticle formation 554 would reveal new aspects of membrane trafficking supporting cuticle formation. 555 CYCLASE-ASSOCIATED PROTEIN1, encoded by the candidate gene on 556 chromosome 1 near peak SNP S1_270551534, has not been associated with cuticle 23 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 557 development in previous studies. CAPs are actin monomer-binding proteins that regulate 558 actin dynamics (Qualmann et al. 2000; Hubberstey and Mottillo 2002). Functions for CAPs 559 in actin regulation and actin-dependent processes in Arabidopsis have been demonstrated 560 through genetic and biochemical studies (Barrero et al. 2002; Chaudhry et al. 2007; Deeks et 561 al. 2007). Actin dynamics have been implicated in multiple aspects of membrane trafficking 562 in eukaryotic cells including vesicle formation and vesicle movement (Lanzetti 2007). Thus, 563 CAP could impact gc in maize via regulation of membrane trafficking events supporting 564 cuticle formation such as Golgi-mediated secretion (McFarlane et al. 2014; Luo et al. 2019) 565 and/or endocytosis. 566 Two of the candidate genes we identified have predicted functions in cuticle 567 formation. One of these is CER7, which functions in Arabidopsis to regulate the biosynthesis 568 of alkanes, one of the major classes of cuticle waxes in this species (Jenks et al. 2002) and 569 also in maize adult leaves(Bourgault et al. 2019). CER7 encodes an exosomal 3′-to-5′ 570 exoribonuclease that regulates the expression level of CER3 (Hooker et al. 2007). The CER3 571 protein acts together with CER1 to catalyze the formation of alkanes (Bernard et al. 2012). 572 The other candidate with a predicted function in cuticle formation encodes a GDSL lipase 573 (Zm00001d011661). GDSL lipases catalyze the hydrolysis of mono-, di- and triglycerols and 574 release free fatty acids and alcohols (Angkawidjaja and Kanaya 2006). Although GDSL 575 lipases are involved in a large number of biological processes in plants, their function in the 576 biosynthesis of cutin is well-recognized(Takahashi et al. 2010; Girard et al. 2012; Yeats et al. 577 2012b). Moreover, a homolog of Zm00001d011661 in Arabidopsis, AT3G16370, was down- 578 regulated along with other cuticle associated genes in a transgenic DESPERADO/ABCG11 579 silenced line with changes in levels of cutin monomers (Panikashvili et al. 2010), whereas its 580 homolog in Citrus clementina was reported as the top differentially expressed gene in the 581 fruit epidermis(Matas et al. 2010). 24 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 582 The final candidate gene identified encodes the cell wall biosynthesis enzyme β-D- 583 XYLOSIDASE 4 (Zm00001d026415). Sometimes the cuticle is described as a specialized 584 lipidic modification of the cell wall, and polysaccharides are known to be deposited in some 585 regions of the cuticle (Yeats and Rose 2013). Transmission electron microscopy and non- 586 destructive polarization modulation-infrared reflection-absorption spectroscopy have 587 challenged the long-standing notion that polysaccharide compounds, mostly pectin, are 588 confined to the cuticular layer of the cuticle, adjacent to the cell wall. Other polysaccharides, 589 hemicelluloses such as xylan and xyloglucan, were found close to the surface of the cuticle, 590 within the cuticle proper of several species (Guzmán et al. 2014; Hama et al. 2019). 591 Therefore, cell wall-modifying enzymes could have an influence on the polysaccharide 592 content of the cuticle, or might impact cuticle composition or organization by affecting the 593 transit of cuticle lipids across the wall to the cuticle (Müller et al. 2013). 594 Genes encoding CAP, ISTL1, GDSL lipase, BETA-D-XYLOSIDASE 4, and both 595 SEC14 homologs were declared expressed in epidermal cells of developing maize leaves 596 (Figure 3), where the cuticle matures, indicating that these genes are involved in metabolic 597 and other cellular processes in maize leaf epidermal cells and possibly associated with cuticle 598 development. The abundance of the CER7 homolog transcript, a core subunit of the exosome 599 that regulates wax biosynthesis (Hooker et al. 2007), was too low to be declared expressed, 600 with cpm values ranging from 0.033 to 0.137 in 6 out of the 20 samples. Similarly, very low 601 abundances of this transcript were observed by RNAseq in a number of different tissues from 602 maize inbred B73 (Stelpflug et al. 2016). The increasing/decreasing trends in transcript 603 abundance along the developmental gradient (Figure 3) are consistent with regulatory roles in 604 cuticle maturation; however, these trends are equally consistent with roles in other 605 developmental changes taking place at the same time. Additional work will be needed to 606 definitively link changes in expression levels of these genes with variation in gc. 25 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 607 In our GWAS results, even after controlling for the confounding impact of flowering 608 time, each genomic region significantly associated with gc explained 4-6% (R2LR-SNP - R2LR) of 609 the variation in the panel (Table 2), thus only accounting for a minor fraction of the estimated 610 heritability. If the genetic architecture of gc is predominated by a large number of loci with 611 small allelic effect sizes, then it is possible that our association population was not of a 612 sample size sufficient to offer the statistical power needed to detect functional variants with 613 small effects (Long and Langley 1999). With only a total of 235,004 SNPs with common 614 minor allele frequencies (5% or greater) scored on the diversity panel, the “missing” 615 heritability could be also partially attributed to not having a density of SNP markers required 616 to achieve strong LD (r2 > 0.80) with unscored causal variants that are potentially rare in 617 frequency (Myles et al. 2009; Buckler et al. 2009). Therefore, increasing SNP marker density 618 and the mapping population size are likely important for explaining the missing heritability of 619 gc, which is possibly controlled by large numbers of small effect alleles in maize considering 620 the complex cuticle biosynthesis network and transport mechanisms (Yeats and Rose 2013), 621 as well as factors that affect leaf water trafficking. 622 In general, the predictive abilities of WGP were moderately high for the tested nine 623 training-prediction combinations (Figure 4; Table S10), indicating that genetic gain for gc 624 could be accelerated with WGP in maize breeding programs. Presumably due to influence of 625 contrasting environmental conditions on the phenotypic data collected at the two field 626 locations (Figure S5; Tables S1 and S4), the prediction abilities were lowest when the MA 627 and SD training sets predicted gc each other (Figure 4; Table S10). However, AllEnv showed 628 similar predictive abilities for MA and SD as for themselves, because AllEnv contained 629 phenotypic information from both locations. Given the range in phenotypic variation 630 observed for gc across locations, it is highly recommended that phenotyping and selection 631 occur in the specific target breeding environments, especially if the focus is on hot, arid 26 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 632 environments such as MA. Furthermore, if gc does ultimately have a polygenic genetic basis, 633 genomic selection would be the optimal breeding strategy for gc and inevitably drought 634 tolerance in maize breeding populations rather than a marker-assisted selection approach 635 designed for only a few loci with large effects (Meuwissen et al. 2001; Lorenz et al. 2011). 636 637 Acknowledgements 638 We especially thank Aaron Waybright, Christine Caron, Angel Mendoza, Albert Nguyen, 639 Indira Oueralta Castillo, Anastasia Zagordo, and the BICD 101 students at UCSD in the 640 summer of 2016 and 2017, for collecting phenotypic data. We thank Mark Millard and 641 Candice Gardner of the USDA-ARS North Central Regional Plant Introduction Station 642 (NCRPIS) in Ames, Iowa, and Candice Hirsch at the University of Minnesota for providing 643 seed of the Wisconsin Diversity panel. We also thank Bill Luckett, Andrew French, Kelly 644 Thorp, Alison Thompson, John Dyer and others at the U.S. Arid-Land Agricultural Research 645 Center in Maricopa, AZ, for their assistance with planting and providing a facility for 646 phenotyping. Additionally, we thank Clint Jones, Greg Main, Russell Noon, Rick Ward and 647 others at the University of Arizona, Maricopa Agricultural Center in Maricopa, AZ, for the 648 management of the Arizona field trials. This research was supported by the National Science 649 Foundation IOS1444507. 650 27 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 651 Literature Cited 652 Aharoni, A., S. Dixit, R. Jetter, E. Thoenes, G. van Arkel et al., 2004 The SHINE clade 653 of AP2 domain transcription factors activates wax biosynthesis, alters cuticle 654 properties, and confers drought tolerance when overexpressed in Arabidopsis. 655 Plant Cell 16: 2463–2480. 656 657 658 659 660 661 662 Anders, S., P. T. Pyl, and W. Huber, 2015 HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics 31: 166–169. Angkawidjaja, C., and S. Kanaya, 2006 Family I.3 lipase: bacterial lipases secreted by the type I secretion system. Cell. Mol. Life Sci. 63: 2804–2817. Baker, E. A., 1974 The influence of environment on leaf wax development in Brassica oleracea var. gemmifera. New Phytol. 73: 955–966. Balcer, H. I., A. L. Goodman, A. A. Rodal, E. Smith, J. Kugler et al., 2003 663 Coordinated regulation of actin filament turnover by a high-molecular-weight 664 Srv2/CAP complex, cofilin, profilin, and Aip1. Curr. Biol. 13: 2159–2169. 665 Barrero, R. A., M. Umeda, S. Yamamura, and H. Uchimiya, 2002 Arabidopsis CAP 666 regulates the actin cytoskeleton necessary for plant cell elongation and division. 667 Plant Cell 14: 149–163. 668 Baseggio, M., M. Murray, M. Magallanes-Lundback, N. Kaczmar, J. Chamness et al., 669 2019 Genome-wide association and genomic prediction models of 670 tocochromanols in fresh sweet corn kernels. Plant Genome 12: 1-17. 671 Beisson, F., Y. Li-Beisson, and M. Pollard, 2012 Solving the puzzles of cutin and 672 suberin polymer biosynthesis. Curr. Opin. Plant Biol. 15: 329–337. 673 Benjamini, Y., and Y. Hochberg, 1995 Controlling the false discovery rate: A practical 674 and powerful approach to multiple testing. J. R. Stat. Soc. Series B Stat. 675 Methodol. 57: 289–300. 28 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 676 Bernard, A., F. Domergue, S. Pascal, R. Jetter, C. Renne et al., 2012 Reconstitution of 677 plant alkane biosynthesis in yeast demonstrates that Arabidopsis 678 ECERIFERUM1 and ECERIFERUM3 are core components of a very-long- 679 chain alkane synthesis complex. Plant Cell 24: 3106–3118. 680 Bessire, M., S. Borel, G. Fabre, L. Carraça, N. Efremova et al., 2011 A member of the 681 PLEIOTROPIC DRUG RESISTANCE family of ATP binding cassette 682 transporters is required for the formation of a functional cuticle in Arabidopsis. 683 Plant Cell 23: 1958–1970. 684 Bi, H., S. Luang, Y. Li, N. Bazanova, S. Morran et al., 2016 Identification and 685 characterization of wheat drought-responsive MYB transcription factors 686 involved in the regulation of cuticle biosynthesis. J. Exp. Bot. 67: 5363–5380. 687 Bird, D., F. Beisson, A. Brigham, J. Shin, S. Greer et al., 2007 Characterization of 688 Arabidopsis ABCG11/WBC11, an ATP binding cassette (ABC) transporter that 689 is required for cuticular lipid secretion. Plant J. 52: 485–498. 690 691 692 Borisjuk, N., M. Hrmova, and S. Lopato, 2014 Transcriptional regulation of cuticle biosynthesis. Biotechnol. Adv. 32: 526–540. Bourdenx, B., A. Bernard, F. Domergue, S. Pascal, A. Léger et al., 2011 693 Overexpression of Arabidopsis ECERIFERUM1 promotes wax very-long-chain 694 alkane biosynthesis and influences plant response to biotic and abiotic stresses. 695 Plant Physiol. 156: 29–45. 696 Bourgault, R., S. Matschi, M. Vasquez, P. Qiao, A. Sonntag et al., 2019 Constructing 697 functional cuticles: Analysis of relationships between cuticle lipid composition, 698 ultrastructure and water barrier function in developing adult maize leaves. 699 Annals of Botany. 29 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 700 Bradbury, P. J., Z. Zhang, D. E. Kroon, T. M. Casstevens, Y. Ramdoss et al., 2007 701 TASSEL: software for association mapping of complex traits in diverse 702 samples. Bioinformatics 23: 2633–2635. 703 Buckler, E. S., J. B. Holland, P. J. Bradbury, C. B. Acharya, P. J. Brown et al., 2009 704 The genetic architecture of maize flowering time. Science 325: 714–718. 705 Buono, R. A., J. Paez-Valencia, N. D. Miller, K. Goodman, C. Spitzer et al., 2016 Role 706 of SKD1 regulators LIP5 and IST1-LIKE1 in endosomal sorting and plant 707 development. Plant Physiol. 171: 251–264. 708 Cameron, K. D., M. A. Teece, and L. B. Smart, 2006 Increased accumulation of 709 cuticular wax and expression of lipid transfer protein in response to periodic 710 drying events in leaves of tree tobacco. Plant Physiol. 140: 176–183. 711 Chaudhry, F., C. Guérin, M. von Witsch, L. Blanchoin, and C. J. Staiger, 2007 712 Identification of Arabidopsis cyclase-associated protein 1 as the first nucleotide 713 exchange factor for plant actin. Mol. Biol. Cell 18: 3002–3014. 714 Chen, X., S. M. Goodwin, V. L. Boroff, X. Liu, and M. A. Jenks, 2003 Cloning and 715 characterization of the WAX2 gene of Arabidopsis involved in cuticle 716 membrane and wax production. Plant Cell 15: 1170–1185. 717 Chen, G., T. Komatsuda, J. F. Ma, C. Nawrath, M. Pourkheirandish et al., 2011 An 718 ATP-binding cassette subfamily G full transporter is essential for the retention 719 of leaf water in both wild barley and rice. Proc. Natl. Acad. Sci. USA 108: 720 12354–12359. 721 722 de Campos, M. K., and G. Schaaf, 2017 The regulation of cell polarity by lipid transfer proteins of the SEC14 family. Curr. Opin. Plant Biol. 40: 158–168. 30 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 723 DeBono, A., T. H. Yeats, J. K. C. Rose, D. Bird, R. Jetter et al., 2009 Arabidopsis 724 LTPG is a glycosylphosphatidylinositol-anchored lipid transfer protein required 725 for export of lipids to the plant surface. Plant Cell 21: 1230–1238. 726 Deeks, M. J., C. Rodrigues, S. Dimmock, T. Ketelaar, S. K. Maciver et al., 2007 727 Arabidopsis CAP1 - a key regulator of actin organisation and development. J. 728 Cell Sci. 120: 2609–2618. 729 Elshire, R. J., J. C. Glaubitz, Q. Sun, J. A. Poland, K. Kawamoto et al., 2011 A robust, 730 simple genotyping-by-sequencing (GBS) approach for high diversity species. 731 PLoS ONE 6: e19379. 732 733 734 735 736 Endelman, J. B., 2011 Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4: 250-255. Endelman, J. B., and J.-L. Jannink, 2012 Shrinkage estimation of the realized relationship matrix. G3 2: 1405–1413. Febrero, A., S. Fernandez, J. L. Molina-Cano, and J. L. Araus, 1998 Yield, carbon 737 isotope discrimination, canopy reflectance and cuticular conductance of barley 738 isolines of differing glaucousness. J. Exp. Bot. 49: 1575–1581. 739 Fich, E. A., N. A. Segerson, and J. K. C. Rose, 2016 The plant polyester cutin: 740 biosynthesis, structure, and biological roles. Annu. Rev. Plant Biol. 67: 207– 741 233. 742 Gilmour, A. R., B. J. Gogel, B. R. Cullis, R. Thompson, D. Butler et al., 2009 ASReml 743 user guide release 3.0. VSN International Ltd, Hemel Hempstead, UK. 744 Girard, A.-L., F. Mounet, M. Lemaire-Chamley, C. Gaillard, K. Elmorjani et al., 2012 745 Tomato GDSL1 is required for cutin deposition in the fruit cuticle. Plant Cell 746 24: 3119–3134. 31 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 747 748 749 Grant, R. F., B. S. Jackson, J. R. Kiniry, and G. F. Arkin, 1989 Water deficit timing effects on yield components in maize. Agron. J. 81: 61–65. Guo, J., W. Xu, X. Yu, H. Shen, H. Li et al., 2016 Cuticular wax accumulation is 750 associated with drought tolerance in wheat near-isogenic lines. Front. Plant Sci. 751 7: 1809. 752 Guzmán, P., V. Fernández, M. L. García, M. Khayet, A. Fernández et al., 2014 753 Localization of polysaccharides in isolated and intact cuticles of eucalypt, 754 poplar and pear leaves by enzyme-gold labelling. Plant Physiol. Biochem. 76: 755 1–6. 756 Hama, T., K. Seki, A. Ishibashi, A. Miyazaki, A. Kouchi et al., 2019 Probing the 757 molecular structure and orientation of the leaf surface of Brassica oleracea L. 758 by polarization modulation-infrared reflection-absorption spectroscopy. Plant 759 Cell Physiol. 60: 1567–1580. 760 Hansey, C. N., J. M. Johnson, R. S. Sekhon, S. M. Kaeppler, and N. de Leon, 2011 761 Genetic diversity of a maize association population with restricted phenology. 762 Crop Sci. 51: 704. 763 Hen-Avivi, S., O. Savin, R. C. Racovita, W.-S. Lee, N. M. Adamski et al., 2016 A 764 metabolic gene cluster in the wheat W1 and the barley Cer-cqu loci determines 765 β-diketone biosynthesis and glaucousness. Plant Cell 28: 1440–1460. 766 Hill, C. P., and M. Babst, 2012 Structure and function of the membrane deformation 767 AAA ATPase Vps4. Biochim. Biophys. Acta 1823: 172–181. 768 Hill, W. G., and B. S. Weir, 1988 Variances and covariances of squared linkage 769 disequilibria in finite populations. Theor. Popul. Biol. 33: 54–78. 770 771 Holland, J. B., W. E. Nyquist, and C. T. Cervantes-Martínez, 2010 Estimating and interpreting heritability for plant breeding: an update. Plant Breed. Rev. 9–112. 32 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 772 Hooker, T. S., P. Lam, H. Zheng, and L. Kunst, 2007 A core subunit of the RNA- 773 processing/degrading exosome specifically influences cuticular wax 774 biosynthesis in Arabidopsis. Plant Cell 19: 904–913. 775 Huang, J., R. Ghosh, and V. A. Bankaitis, 2016 Sec14-like phosphatidylinositol 776 transfer proteins and the biological landscape of phosphoinositide signaling in 777 plants. Biochim. Biophys. Acta 1861: 1352–1364. 778 779 780 Hubberstey, A. V., and E. P. Mottillo, 2002 Cyclase-associated proteins: CAPacity for linking signal transduction and actin polymerization. FASEB J. 16: 487–499. Hung, H.-Y., L. M. Shannon, F. Tian, P. J. Bradbury, C. Chen et al., 2012 ZmCCT and 781 the genetic basis of day-length adaptation underlying the postdomestication 782 spread of maize. Proc. Natl. Acad. Sci. USA 109: E1913–21. 783 784 785 786 787 788 789 790 Jenks, M. A., S. D. Eigenbrode, and B. Lemieux, 2002 Cuticular waxes of Arabidopsis. Arabidopsis Book 1: e0016. Jetter, R., L. Kunst, and A. L. Samuels, 2008 Composition of plant cuticular waxes. Bio. Plant Cuticle 23: 145–181. Kerstiens, G., 2006 Water transport in plant cuticles: an update. J. Exp. Bot. 57: 2493– 2499. Kim, D., B. Langmead, and S. L. Salzberg, 2015 HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12: 357–360. 791 Kim, H., S. B. Lee, H. J. Kim, M. K. Min, I. Hwang et al., 2012 Characterization of 792 glycosylphosphatidylinositol-anchored lipid transfer protein 2 (LTPG2) and 793 overlapping function between LTPG/LTPG1 and LTPG2 in cuticular wax 794 export or accumulation in Arabidopsis thaliana. Plant Cell Physiol. 53: 1391– 795 1403. 33 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 796 Koch, K., K. D. Hartmann, L. Schreiber, W. Barthlott, and C. Neinhuis, 2006 797 Influences of air humidity during the cultivation of plants on wax chemical 798 composition, morphology and leaf surface wettability. Environ. Exp. Bot. 56: 799 1–9. 800 Kosma, D. K., B. Bourdenx, A. Bernard, E. P. Parsons, S. Lü et al., 2009 The impact 801 of water deficiency on leaf cuticle lipids of Arabidopsis. Plant Physiol. 151: 802 1918–1929. 803 Kosma, D. K., and M. A. Jenks, 2007 Eco-physiological and molecular-genetic 804 determinants of plant cuticle function in drought and salt stress tolerance, pp. 805 91–120 in Advances in Molecular Breeding Toward Drought and Salt Tolerant 806 Crops, edited by M. A. Jenks, P. M. Hasegawa and S. M. Jain. Springer, 807 Dordrecht. 808 Kurata, T., C. Kawabata-Awai, E. Sakuradani, S. Shimizu, K. Okada et al., 2003 The 809 YORE-YORE gene regulates multiple aspects of epidermal cell differentiation 810 in Arabidopsis. Plant J. 36: 55–66. 811 812 Kurtz, S., 2003 The Vmatch large scale sequence analysis software. Ref Type: Computer Program 412: 297. 813 Lanzetti, L., 2007 Actin in membrane trafficking. Curr. Opin. Cell Biol. 19: 453–458. 814 Lee, S. B., and M. C. Suh, 2015 Advances in the understanding of cuticular waxes in 815 816 817 818 Arabidopsis thaliana and crop species. Plant Cell Rep. 34: 557–572. Li, L. Li, Y. Du, C. He, C. R. Dietrich et al., 2019 Maize glossy6 is involved in cuticular wax deposition and drought tolerance. J. Exp. Bot. 70: 3089–3099. Lipka, A. E., C. B. Kandianis, M. E. Hudson, J. Yu, J. Drnevich et al., 2015 From 819 association to prediction: statistical methods for the dissection and selection of 820 complex traits in plants. Curr. Opin. Plant Biol. 24: 110–118. 34 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 821 822 823 Lipka, A. E., F. Tian, Q. Wang, J. Peiffer, M. Li et al., 2012 GAPIT: genome association and prediction integrated tool. Bioinformatics 28: 2397–2399. Littell, R.C., G. A. Milliken, W. W. Stroup, R.D. Wolfinger, and O. Schabenberger, 824 2006 Appendix 1: Linear mixed model theory, pp. 733–756 in SAS for Mixed 825 Models. SAS Institute Inc., Cary, NC. 826 Long, A. D., and C. H. Langley, 1999 The power of association studies to detect the 827 contribution of candidate genetic loci to variation in complex traits. Genome 828 Res. 9: 720–731. 829 Lorenz, A. J., S. Chao, F. G. Asoro, E. L. Heffner, T. Hayashi et al., 2011 Genomic 830 selection in plant breeding: knowledge and prospects, pp. 77–123 in Advances 831 in agronomy, edited by D. L. Sparks. Elsevier Inc., Amsterdam. 832 Luo, Z., P. Tomasi, N. Fahlgren, and H. Abdel-Haleem, 2019 Genome-wide 833 association study (GWAS) of leaf cuticular wax components in Camelina 834 sativa identifies genetic loci related to intracellular wax transport. BMC Plant 835 Biol. 19: 187. 836 837 838 Lynch, M., B. Walsh, and Others, 1998 Genetics and analysis of quantitative traits. Sinauer Sunderland, MA. Matas, A. J., J. Agustí, F. R. Tadeo, M. Talón, and J. K. C. Rose, 2010 Tissue-specific 839 transcriptome profiling of the citrus fruit epidermis and subepidermis using 840 laser capture microdissection. J. Exp. Bot. 61: 3321–3330. 841 McFarlane, H. E., Y. Watanabe, W. Yang, Y. Huang, J. Ohlrogge et al., 2014 Golgi- 842 and trans-golgi network-mediated vesicle trafficking is required for wax 843 secretion from epidermal cells. Plant Physiol. 164: 1250–1260. 844 845 Meuwissen, T. H., B. J. Hayes, and M. E. Goddard, 2001 Prediction of total genetic value using genome-wide dense marker maps. Genetics 157: 1819–1829. 35 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 846 Müller, K., G. Levesque-Tremblay, A. Fernandes, A. Wormit, S. Bartels et al., 2013 847 Overexpression of a pectin methylesterase inhibitor in Arabidopsis thaliana 848 leads to altered growth morphology of the stem and defective organ separation. 849 Plant Signal. Behav. 8: e26464. 850 Myles, S., J. Peiffer, P. J. Brown, E. S. Ersoz, Z. Zhang et al., 2009 Association 851 mapping: critical considerations shift from genotyping to experimental design. 852 Plant Cell 21: 2194–2202. 853 854 855 Neter, J., M. H. Kutner, C. J. Nachtsheim, and W. Wasserman, 1996 Applied linear statistical models. McGraw-Hill, Boston, MA. Ono, S., 2013 The role of cyclase-associated protein in regulating actin filament 856 dynamics--more than a monomer-sequestration factor. J. Cell Sci. 126: 3249– 857 3258. 858 Owens, B. F., A. E. Lipka, M. Magallanes-Lundback, T. Tiede, C. H. Diepenbrock et 859 al., 2014 A foundation for provitamin A biofortification of maize: genome- 860 wide association and genomic prediction models of carotenoid levels. Genetics 861 198: 1699–1716. 862 Panikashvili, D., S. Savaldi-Goldstein, T. Mandel, T. Yifhar, R. B. Franke et al., 2007 863 The Arabidopsis DESPERADO/AtWBC11 transporter is required for cutin and 864 wax secretion. Plant Physiol. 145: 1345–1360. 865 Panikashvili, D., J. X. Shi, S. Bocobza, R. B. Franke, L. Schreiber et al., 2010 The 866 Arabidopsis DSO/ABCG11 transporter affects cutin metabolism in reproductive 867 organs and suberin in roots. Mol. Plant 3: 563–575. 868 869 Petit, J., C. Bres, J. P. Mauxion, and B. Bakan, 2017 Breeding for cuticle-associated traits in crop species: traits, targets, and strategies. J. Exp. Bot. 68: 5369–5387. 36 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 870 871 Pighin, J. A., H. Zheng, L. J. Balakshin, I. P. Goodman, T. L. Western et al., 2004 Plant cuticular lipid export requires an ABC transporter. Science 306: 702–704. 872 Pollard, M., F. Beisson, Y. Li, and J. B. Ohlrogge, 2008 Building lipid barriers: 873 biosynthesis of cutin and suberin. Trends Plant Sci. 13: 236–246. 874 Purcell, S., B. Neale, K. Todd-Brown, L. Thomas, M. A. R. Ferreira et al., 2007 875 PLINK: a tool set for whole-genome association and population-based linkage 876 analyses. Am. J. Hum. Genet. 81: 559–575. 877 Qiao, P., R. Bourgault, M. Mohammadi, L. G. Smith, M. A. Gore et al., 2019 Network 878 analyses implicate a role for PHYTOCHROME-mediated light signaling in the 879 regulation of cuticle development in plant leaves. bioRxiv. doi: 10.1101/812107 880 (Preprint posted October 21, 2019). 881 Qualmann, B., M. M. Kessels, and R. B. Kelly, 2000 Molecular links between 882 endocytosis and the actin cytoskeleton. J. Cell Biol. 150: F111–6. 883 R Core Team (2018). R: A language and environment for statistical computing. R 884 Foundation for 885 Statistical Computing, Vienna, Austria. URL https://www.R-project.org/ 886 Raj, A., M. Stephens, and J. K. Pritchard, 2014 fastSTRUCTURE: variational 887 inference of population structure in large SNP data sets. Genetics 197: 573– 888 589. 889 Riederer, M., and G. Schneider, 1990 The effect of the environment on the 890 permeability and composition of Citrus leaf cuticles. Planta 180: 154-165. 891 Riederer, M., and L. Schreiber, 2001 Protecting against water loss: analysis of the 892 893 894 barrier properties of plant cuticles. J. Exp. Bot. 52: 2023–2032. Ristic, Z., and M. A. Jenks, 2002 Leaf cuticle and water loss in maize lines differing in dehydration avoidance. J. Plant Physiol. 159: 645–651. 37 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 895 Robinson, M. D., D. J. McCarthy, and G. K. Smyth, 2010 edgeR: a Bioconductor 896 package for differential expression analysis of digital gene expression data. 897 Bioinformatics 26: 139–140. 898 Romay, M. C., M. J. Millard, J. C. Glaubitz, J. A. Peiffer, K. L. Swarts et al., 2013 899 Comprehensive genotyping of the USA national maize inbred seed bank. 900 Genome Biol. 14: R55. 901 902 903 Rutter, B. D., and R. W. Innes, 2018 Extracellular vesicles as key mediators of plant– microbe interactions. Curr. Opin. Plant Biol. 44: 16–22. Salvi, S., G. Sponza, M. Morgante, D. Tomes, X. Niu et al., 2007 Conserved 904 noncoding genomic sequences associated with a flowering-time quantitative 905 trait locus in maize. Proc. Natl. Acad. Sci. USA 104: 11376–11381. 906 907 Schneider, C. A., W. S. Rasband, and K. W. Eliceiri, 2012 NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9: 671-675. 908 Schwarz, G., 1978 Estimating the dimension of a model. Ann. Stat. 6: 461–464. 909 Shepherd, T., G. W. Robertson, D. W. Griffiths, and A. N. E. Bircht, 1997 Effects of 910 environment on the composition of epicuticular wax esters from kale and 911 swede. Phytochemistry 46: 83–96. 912 913 914 Shepherd, T., and D. Wynne Griffiths, 2006 The effects of stress on plant cuticular waxes. New Phytol. 171: 469–499. Stelpflug, S. C., R. S. Sekhon, B. Vaillancourt, C. N. Hirsch, C. R. Buell et al., 2016 915 An expanded maize gene expression atlas based on RNA sequencing and its use 916 to explore root development. Plant Genome 9: 1-16. 917 Studer, A., Q. Zhao, J. Ross-Ibarra, and J. Doebley, 2011 Identification of a functional 918 transposon insertion in the maize domestication gene tb1. Nat. Genet. 43: 919 1160–1163. 38 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 920 921 922 Sun, G., C. Zhu, M. H. Kramer, S.-S. Yang, W. Song et al., 2010 Variation explained in mixed-model association mapping. Heredity 105: 333–340. Sutter, E., and R. W. Langhans, 1982 Formation of epicuticular wax and its effect on 923 water loss in cabbage plants regenerated from shoot-tip culture. Can. J. Bot. 60: 924 2896–2902. 925 Swarts, K., H. Li, J. A. Romero Navarro, D. An, M. C. Romay et al., 2014 Novel 926 methods to optimize genotypic imputation for low-coverage, next-generation 927 sequence data in crop plants. Plant Genome 7: 1-12. 928 Takahashi, K., T. Shimada, M. Kondo, A. Tamai, M. Mori et al., 2010 Ectopic 929 expression of an esterase, which is a candidate for the unidentified plant 930 cutinase, causes cuticular defects in Arabidopsis thaliana. Plant Cell Physiol. 931 51: 123–131. 932 Valeska Zeisler-Diehl, V., B. Migdal, and L. Schreiber, 2017 Quantitative 933 characterization of cuticular barrier properties: methods, requirements, and 934 problems. J. Exp. Bot. 68: 5281–5291. 935 Wallace, J. G., P. J. Bradbury, N. Zhang, Y. Gibon, M. Stitt et al., 2014 Association 936 mapping across numerous traits reveals patterns of functional variation in 937 maize. PLoS Genet. 10: e1004845. 938 Wang, Y., L. Wan, L. Zhang, Z. Zhang, H. Zhang et al., 2012 An ethylene response 939 factor OsWR1 responsive to drought stress transcriptionally activates wax 940 synthesis related genes and increases wax production in rice. Plant Mol. Biol. 941 78: 275–288. 942 Xue, D., X. Zhang, X. Lu, G. Chen, and Z.-H. Chen, 2017 Molecular and evolutionary 943 mechanisms of cuticular wax for plant drought tolerance. Front. Plant Sci. 8: 944 621. 39 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 945 Yeats, T. H., G. J. Buda, Z. Wang, N. Chehanovsky, L. C. Moyle et al., 2012a The 946 fruit cuticles of wild tomato species exhibit architectural and chemical 947 diversity, providing a new model for studying the evolution of cuticle function. 948 Plant J. 69: 655–666. 949 Yeats, T. H., L. B. B. Martin, H. M.-F. Viart, T. Isaacson, Y. He et al., 2012b The 950 identification of cutin synthase: formation of the plant polyester cutin. Nat. 951 Chem. Biol. 8: 609–611. 952 953 954 Yeats, T. H., and J. K. C. Rose, 2013 The formation and function of plant cuticles. Plant Physiol. 163: 5–20. Yu, J., G. Pressoir, W. H. Briggs, I. Vroh Bi, M. Yamasaki et al., 2006 A unified 955 mixed-model method for association mapping that accounts for multiple levels 956 of relatedness. Nat. Genet. 38: 203–208. 957 Zhang, J.-Y., C. D. Broeckling, E. B. Blancaflor, M. K. Sledge, L. W. Sumner et al., 958 2005 Overexpression of WXP1, a putative Medicago truncatula AP2 domain- 959 containing transcription factor gene, increases cuticular wax accumulation and 960 enhances drought tolerance in transgenic alfalfa (Medicago sativa). Plant J. 42: 961 689–707. 962 Zhang, Z., E. Ersoz, C.-Q. Lai, R. J. Todhunter, H. K. Tiwari et al., 2010 Mixed linear 963 model approach adapted for genome-wide association studies. Nat. Genet. 42: 964 355–360. 965 966 Zhao, L., and F. D. Sack, 1999 Ultrastructure of stomatal development in Arabidopsis (Brassicaceae) leaves. Am. J. Bot. 86: 929–939. 967 Zhou, L., E. Ni, J. Yang, H. Zhou, H. Liang et al., 2013 Rice OsGL1-6 is involved in 968 leaf cuticular wax accumulation and drought resistance. PLoS One 8: e65139. 40 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 969 Zhu, X., and L. Xiong, 2013 Putative megaenzyme DWA1 plays essential roles in 970 drought resistance by regulating stress-induced wax deposition in rice. Proc. 971 Natl. Acad. Sci. USA 110: 17790–17795. 972 41 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 973 Figure 1. 974 Manhattan plot of results from a genome-wide association study (GWAS) of adult maize leaf 975 cuticular conductance (gc) conducted in Maricopa (MA) and across all four environments 976 (AllEnv). The -log10 P-value of each SNP tested in a mixed linear model analysis of gc is 977 plotted as a point against its physical position (B73 RefGen_v4) for the 10 chromosomes of 978 maize. The least significant single-nucleotide polymorphism (SNP) at a genome-wide false 979 discovery rate of 10% in MA and AllEnv is indicated by a dashed horizontal orange line and 980 a dotted horizontal orange line, respectively. SNPs significantly associated with gc in MA and 981 AllEnv are represented by blue circles and red triangles, respectively. The most plausible 982 candidate genes within ± 200 kb of the significantly associated SNPs are listed above their 983 corresponding GWAS signal. 984 Figure 2. 985 Association of SNP markers with adult maize leaf cuticular conductance (gc) across a 986 genomic region on chromosome 4. Scatter plot of association results from a mixed linear 987 model analysis of gc conducted across all four environments (AllEnv) and linkage 988 disequilibrium (LD) estimates (r2) for a genomic region that contains a gene encoding an 989 INCREASED SALT TOLERANCE1-LIKE1 (ISTL1) protein (Zm00001d049479). The -log10 990 P-values of tested single-nucleotide polymorphisms (SNPs) are represented by vertical lines. 991 Blue vertical lines are SNPs that are statistically significant at a false discovery rate (FDR) of 992 10%. The r2 values of each SNP relative to the peak SNP (indicated by a solid orange 993 triangle) at 31,733,973 bp (B73 RefGen_v4) are indicated by triangles. Open orange triangles 994 represent SNPs with r2 > 0.1 relative to the peak SNP. The least significant SNP at a genome- 995 wide FDR of 10% is indicated by a dashed horizontal orange line. The black dashed vertical 996 line indicates the genomic position of the ISTL1 protein. 997 Figure 3. 42 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 998 Transcript abundance of six candidate genes in maize leaf epidermal cells. Transcript 999 abundance (counts per million, cpm) of each gene is plotted against the developmental 1000 gradient (six intervals from 2–14 cm, and one interval from 20–22 cm) of expanding adult 1001 leaf 8 from maize inbred B73. Error bars represent the standard deviation. 1002 Figure 4. 1003 Predictive abilities of whole-genome prediction for adult maize leaf cuticular conductance 1004 (gc) in locations (Maricopa, MA; and San Diego, SD) and across all four environments 1005 (AllEnv) in a scheme that used each of the three as a training set. Grey, orange and blue 1006 represent the predictive abilities for gc in MA, SD, and AllEnv, respectively. 1007 43 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1008 Table 1. 1009 Averages and ranges for best linear unbiased predictions (BLUPs) of adult maize leaf 1010 cuticular conductance (gc) in a location (MA and SD) and across all four environments 1011 (AllEnv), and estimated heritabilities on a line-mean basis. 1012 Table 2. 1013 Significant SNPs detected at false discovery rate (FDR) of 10% through a genome-wide 1014 association study of adult maize leaf cuticular conductance (gc) in a location (Maricopa, MA; 1015 San Diego, SD) and across all four environments (AllEnv). 1016 44 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1017 Supplemental information 1018 Supplemental Figure S1. 1019 Diffusion conductance for adult maize leaves transferred from daylight into the dark. Adult 1020 maize leaves were excised from field-grown plants in the same manner routinely used for leaf 1021 cuticular conductance (gc) phenotyping, and transferred to a dark room where their cut bases 1022 were submerged in water. Diffusion conductance was measured immediately (time zero), and 1023 at successive 30 min intervals, to determine the amount of time needed to reach minimum 1024 conductance values in the dark. Panels A-I display results for each genotype tested (indicated 1025 at the top of the panel) with each line representing a different plant (1-3 plants tested for most 1026 genotypes and 7 plants for B73, a reference standard harvested each sampling day). All 1027 analyzed inbred lines except for B73, which was the reference standard, were in the top 15% 1028 for gc across all four environments (AllEnv). 1029 Supplemental Figure S2. 1030 The degree of relationship (Pearson’s correlation, r) between the plot-level averages of leaf 1031 surface area (mm2) and leaf dry weight (g) collected from the same adult maize leaves used 1032 for phenotyping cuticular conductance (gc) in Maricopa and San Diego in 2017. 1033 Supplemental Figure S3. 1034 Population structure analysis of the Wisconsin Diversity panel. (A) Principal component 1035 analysis (PCA) plot colored by group classification of lines [Flint, Iodent, Non-Stiff Stalk 1036 (NSS), Stiff Stalk Synthetic (SSS), Popcorn, Sweet corn, Tropical, Unknown (i.e., genotyped 1037 but ambiguous group), and Not available (i.e., not genotyped)] following that of Hansey et al. 1038 (2011). (B) Distruct plot of fastSTRUCTURE subpopulation membership with K = 6. (C) 1039 PCA plot colored according to the fastSTRUCTURE subpopulation membership with K=6. 1040 Individuals with assignment values (Q) less than 50% were classified as admixed. 1041 Populations 1, 2, and 4 predominantly contained lines belonging to the NSS group; 45 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1042 population 3 mostly had lines from the popcorn and sweet corn groups; population 5 mainly 1043 included lines from the tropical group; and population 6 primarily had lines from the SSS 1044 group. 1045 Supplemental Figure S4. 1046 Principal component analysis (PCA) of transcript abundance for the 21 samples collected 1047 from leaf 8 of maize inbred B73 in the LM-RNAseq analysis. The first two principal 1048 components (PCs) correspond to developmental stage (PC1) and tissue type (PC2). Each 1049 point corresponds to an RNAseq sample. From left to right, each color corresponds to a 1050 specific developmental stage (youngest to oldest part of the leaf). From bottom to top, circles 1051 represent epidermal tissues, and triangles represent internal tissues. The orange dot and 1052 triangle points that are encompassed by a black circle were considered an outlier for stage 2 1053 (interval 4 - 6 cm), thus these samples were removed. 1054 Supplemental Figure S5. 1055 Correlation matrix for best linear unbiased predictors (BLUPs) of maize adult leaf cuticular 1056 conductance (gc; g∙h-1∙g-1) and flowering time (FT; days to anthesis) from each of four single 1057 environments [Maricopa (MA) and San Diego (SD) in 2016 and 2017; MA16, SD16, MA17, 1058 SD17], and the first three principal components calculated from the SNP genotype matrix. 1059 Pearson’s correlation coefficients (r) are presented in the upper right triangle, while the 1060 corresponding P-values for the significance of associations (α = 0.05) are displayed below the 1061 diagonal. 1062 Supplemental Figure S6. 1063 Genome-wide association study (GWAS) of adult maize leaf cuticular conductance (gc) 1064 conducted in locations (Maricopa, MA; and San Diego, SD) and across all four environments 1065 (AllEnv). The -log10 P-value of each SNP tested in a mixed linear model analysis of gc is 1066 plotted as a point against its physical position (B73 RefGen_v4) for the 10 chromosomes of 46 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1067 maize. The least significant single-nucleotide polymorphism (SNP) at a genome-wide false 1068 discovery rate of 10% is indicated by a dashed horizontal black line. 1069 Supplemental Figure S7. 1070 Linkage disequilibrium (LD) estimates in the Wisconsin Diversity panel. Mean decay of LD 1071 measured as pairwise r2 from the 235,004 SNP markers over physical distance. Black lines 1072 indicate the distribution of SNP markers at different percentile cutoffs. The grey horizontal 1073 line indicates r2 of 0.1. 1074 Supplemental Figure S8. 1075 Association of SNP markers with adult maize leaf cuticular conductance (gc) across a 1076 genomic region on chromosome 1. Scatter plot of association results from a mixed linear 1077 model analysis of gc conducted in Maricopa (MA) and linkage disequilibrium (LD) estimates 1078 (r2) for a genomic region that contains genes encoding a CYCLASE-ASSOCIATED 1079 PROTEIN1 (CAP, Zm00001d033830), two SEC14-like proteins [SEC14 (3836), 1080 Zm00001d033836; SEC14 (3837), Zm00001d033837], and a homolog of ECERIFERUM7 1081 (CER7, Zm00001d033842). The -log10 P-values of tested single-nucleotide polymorphisms 1082 (SNPs) are represented by vertical lines. Blue vertical lines are SNPs that are statistically 1083 significant at a false discovery rate (FDR) of 10%. The r2 values of each SNP relative to the 1084 peak SNP (indicated by a solid orange triangle) at 275,318,146 bp (B73 RefGen_v4) are 1085 indicated by triangles. Open orange triangles represent SNPs with r2 > 0.1 relative to the peak 1086 SNP. The least significant SNP at a genome-wide FDR of 10% is indicated by a dashed 1087 horizontal orange line. The black dashed vertical line indicates the genomic positions of the 1088 CAP, two SEC14-like proteins, and a CER7 homolog. 1089 Supplemental Figure S9. 1090 Association of SNP markers with adult maize leaf cuticular conductance (gc) across a 1091 genomic region on chromosome 8. Scatter plot of association results from a mixed linear 47 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1092 model analysis of gc conducted in Maricopa (MA) and linkage disequilibrium (LD) estimates 1093 (r2) for a genomic region that contains a gene encoding a protein of the GLY-ASP-SER-LEU 1094 ESTERASE/LIPASE (GDSL lipase, Zm00001d011661) on chromosome 8. The -log10 P- 1095 values of tested single-nucleotide polymorphisms (SNPs) are represented by vertical lines. 1096 Blue vertical lines are SNPs that are statistically significant at a false discovery rate (FDR) of 1097 10%. The r2 values of each SNP relative to the peak SNP (indicated by a solid orange 1098 triangle) at 157,645,473 bp (B73 RefGen_v4) are indicated by triangles. Open orange 1099 triangles represent SNPs with r2 > 0.1 relative to the peak SNP. The least significant SNP at a 1100 genome-wide FDR of 10% is indicated by a dashed horizontal orange line. The black dashed 1101 vertical line indicates the genomic position of GDSL lipase. 1102 Supplemental Figure S10. 1103 Association of SNP markers with adult maize leaf cuticular conductance (gc) across a 1104 genomic region on chromosome 7. Scatter plot of association results from a mixed linear 1105 model analysis of gc conducted across all four environments (AllEnv) and linkage 1106 disequilibrium (LD) estimates (r2) for a genomic region that contains the peak SNP. The - 1107 log10 P-values of tested single-nucleotide polymorphisms (SNPs) are represented by vertical 1108 lines. Blue vertical lines are SNPs that are statistically significant at a false discovery rate 1109 (FDR) of 10%. The r2 values of each SNP relative to the peak SNP (indicated by a solid 1110 orange triangle) at 23,984,279 bp (B73 RefGen_v4) are indicated by triangles. Open orange 1111 triangles represent SNPs with r2 > 0.1 relative to the peak SNP. The least significant SNP at a 1112 genome-wide FDR of 10% is indicated by a dashed horizontal orange line. 1113 Supplemental Figure S11. 1114 Association of SNP markers with adult maize leaf cuticular conductance (gc) across a 1115 genomic region on chromosome 10. Scatter plot of association results from a mixed linear 1116 model analysis of gc conducted across all four environments (AllEnv) and linkage 48 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1117 disequilibrium (LD) estimates (r2) for a genomic region that contains a gene encoding a β-D- 1118 XYLOSIDASE 4 (Zm00001d026415). The -log10 P-values of tested single-nucleotide 1119 polymorphisms (SNPs) are represented by vertical lines. Blue vertical lines are SNPs that are 1120 statistically significant at a false discovery rate (FDR) of 10%. The r2 values of each SNP 1121 relative to the peak SNP (indicated by a solid orange triangle) at 145,399,603 bp (B73 1122 RefGen_v4) are indicated by triangles. Open orange triangles represent SNPs with r2 > 0.1 1123 relative to the peak SNP. The least significant SNP at a genome-wide FDR of 10% is 1124 indicated by a dashed horizontal orange line. The black dashed vertical line indicates the 1125 genomic position of β-D-XYLOSIDASE 4. 1126 Supplemental Figure S12. 1127 Heatmap of average log2 transformed transcript abundance (counts per million; cpm) for the 1128 24 candidate genes declared to be expressed in epidermal cells of expanding adult leaf 8 of 1129 maize inbred B73. A logarithmic scale was used due to the large range of cpm values across 1130 genes, and a constant with a value of 1 was added to all cpm values (cpm + 1) to prevent the 1131 generation of negative values after log2 transformation. Candidate genes are grouped by the 1132 five genomic regions associated with gc, which is indicated by the peak SNP of each genomic 1133 region. 49 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1134 Supplemental Table S1. 1135 Summary of weather conditions for each field experiment conducted in Maricopa (MA) and 1136 San Diego (SD) in 2016 and 2017 (MA16, SD16, MA17 and SD17) from the day of planting 1137 to the last day of phenotyping leaf cuticular conductance (gc). 1138 Supplemental Table S2. 1139 Plot-level averages of leaf dry weight and leaf surface area for 451 maize inbred lines from 1140 Maricopa and San Diego in 2017 (MA17 and SD17). For each inbred line, an average of each 1141 trait was calculated based on all measured leaves from each plot. 1142 Supplemental Table S3. 1143 Best fitted models used to calculate best linear unbiased predictors (BLUPs) for adult maize 1144 leaf cuticular conductance (gc) and flowering time (days to anthesis) for each of four single 1145 environments [Maricopa (MA) and San Diego (SD) in 2016 and 2017; MA16, SD16, MA17, 1146 SD17], in a location (MA and SD), and across all four environments (AllEnv) according to a 1147 likelihood ratio test (α = 0.05). The star (*) indicates that a random effect term was retained 1148 in the mixed linear model, whereas the x indicates that a fitted random effect term was not 1149 significant and removed from the mixed linear model. The hyphen (-) indicates that a random 1150 effect term was not fitted in a mixed linear model. 1151 Supplemental Table S4. 1152 Best linear unbiased predictors (BLUPs) of adult maize leaf cuticular conductance (gc; g∙h- 1153 1 1154 single environments [Maricopa (MA) and San Diego (SD) in 2016 and 2017; MA16, SD16, 1155 MA17, SD17], in a location (MA and SD), and across all four environments (AllEnv). 1156 Supplemental Table S5. 1157 Variance component estimates from the fitted full mixed linear model used to calculate 1158 heritability on a line-mean basis for leaf cuticular conductance (gc) and flowering time (days ∙g-1) and flowering time (FT; days to anthesis) for 451 maize inbred lines from each of four 50 bioRxiv preprint doi: https://doi.org/10.1101/835892. this version posted November 9, 2019. The copyright holder for this preprint (which was not certified by peer review) is the author/funder. It is made available under a CC-BY-NC-ND 4.0 International license. 1159 to anthesis) in a location (Maricopa, MA; and San Diego, SD) and across all four 1160 environments (AllEnv). 1161 Supplemental Table S6. 1162 Population structure analysis of the Wisconsin Diversity panel. Included in the table are the 1163 first three principal components (PCs) calculated from the SNP genotype matrix (PC1, PC2, 1164 and PC3), assignment value (Q) for each subpopulation (subpopulations 1 to 6), group 1165 classification of lines following assignments of Hansey et al. (2011), subpopulation with the 1166 highest assignment value (K), greatest Q assignment value, and the subpopulation (subpop. 1, 1167 2, 3, 4, 5, 6, or admixed) to which each line was assigned (Category) in our study. 1168 Supplemental Table S7. 1169 Averages and ranges for best linear unbiased predictors (BLUPs) of flowering time (days to 1170 anthesis) in a location (Maricopa, MA; and San Diego, SD) and across all four environments 1171 (AllEnv), and estimated heritabilities on a line-mean basis. 1172 Supplemental Table S8. 1173 Genomic information (RefGen_v4) for the 76 candidate genes within ± 200 kb of the peak 1174 SNPs identified in the genome-wide association study of adult maize leaf cuticular 1175 conductance (gc) in a location (Maricopa, MA; San Diego, SD) and across all four 1176 environm