Allelic diversity of NAC18.1 is a major determinant of fruit firmness and harvest date in apple (Malus domestica)

Softening is a hallmark of ripening in fleshy fruits, and has both desirable and undesirable implications for texture and postharvest stability. Accordingly, the timing and extent of ripening and associated textural changes are key targets for improving fruit quality through breeding. Previously, we identified a large effect locus associated with harvest date and firmness in apple (Malus domestica) using genome-wide association studies (GWAS). Here, we present additional evidence that polymorphisms in or around a transcription factor gene, NAC18.1, cause variation in these traits. First, we confirmed our previous findings with new phenotype and genotype data from ~800 apple accessions. In this population, we compared NAC18.1 to three other ripening-related markers currently used by breeders (ACS1, ACO1, and PG1), and found that the effect of the NAC18.1 genotype on both traits greatly exceeded that observed for the other markers. By sequencing NAC18.1 across 18 accessions, we revealed two predominant haplotypes containing the SNP previously identified using GWAS, as well as dozens of additional SNPs and indels in both the coding and promoter sequences. NAC18.1 encodes a protein with high similarity to the NON-RIPENING (NOR) transcription factor, an early regulator of ripening in tomato (Solanum lycopersicum). To test whether these genes are functionally orthologous, we introduced NAC18.1 transgenes into the tomato nor mutant and showed that both haplotypes complement the nor ripening deficiency. Taken together, these results indicate that polymorphisms in NAC18.1 underlie substantial variation in apple firmness and harvest time through modulation of a conserved ripening program. Highlight NAC18.1 is a member of a family of conserved transcriptional regulators of ripening that underlies variation in fruit firmness and harvest date in diverse apple accessions.


Highlight: 23
NAC18.1 is a member of a family of conserved transcriptional regulators of ripening that 24 underlies variation in fruit firmness and harvest date in diverse apple accessions. 25 26

Abstract: 27
Softening is a hallmark of ripening in fleshy fruits, and has both desirable and 28 undesirable implications for texture and postharvest stability. Accordingly, the timing and extent 29 of ripening and associated textural changes are key targets for improving fruit quality through 30 breeding. Previously, we identified a large effect locus associated with harvest date and firmness 31 in apple (Malus domestica) using genome-wide association studies (GWAS). Here, we present 32 additional evidence that polymorphisms in or around a transcription factor gene, NAC18.1, cause 33 variation in these traits. First, we confirmed our previous findings with new phenotype and 34 genotype data from ~800 apple accessions. In this population, we compared NAC18.1 to three 35 other ripening-related markers currently used by breeders (ACS1, ACO1, and PG1), and found 36 that the effect of the NAC18.1 genotype on both traits greatly exceeded that observed for the 37 other markers. By sequencing NAC18.1 across 18 accessions, we revealed two predominant 38 haplotypes containing the SNP previously identified using GWAS, as well as dozens of 39

Introduction 66
Despite their diverse structure, ontogeny, and biochemical composition, fleshy fruits 67 from a taxonomically broad range of species undergo coordinated ripening programs that have 68 many features in common. Ripening involves numerous physiological and biochemical changes 69 that render the fruit attractive and nutritious for consumption by seed-dispersing animals, or 70 human consumers in the case of cultivated crops. These include the ripening-associated 71 accumulation of sugars, pigments, and flavor or aroma compounds, as well as a loss of flesh 72 firmness due in large part to the controlled modification and depolymerization of cell wall 73

Apple Firmness and Harvest Date 133
The Apple Biodiversity Collection (ABC) is an apple orchard located in Kentville, Nova 134 Scotia, which contains 1,113 accessions. It was first established in spring 2012, when trees were 135 budded onto M.9 rootstocks. In fall 2012, the trees were uprooted and kept in cold storage, 136 before planting in spring 2013. The trees were planted in an incomplete block design, which 137 includes 1 of 3 standards per grid, allowing for correction of positional effects using a REstricted 138 Maximum Likelihood (REML) model, described in Migicovsky et al. (2018). 139 In 2017, we evaluated harvest date for 1,348 trees and fruit firmness for 1,328 trees 140 within the ABC orchard. Due to the diversity of apples within the collection, a variety of 141 methods were used to determine the appropriate time to harvest. First, we observed if the tree 142 had dropped fruit or, for red apples, if the fruit were a deep red color. Next, a sample apple was 143 taken from each tree and touched to assess firmness, tasted to assess starch and sweetness, cut in 144 half to check browning of seeds, and then sprayed with iodine solution to assess starch content 145 (Blanpied and Silsby, 1992). 146 When fruit were determined to be mature, they were harvested and evaluated for 147 firmness. We recorded the firmness (kg cm -2 ) of 5 fruit per tree using a penetrometer with a 1 cm 148 diameter (Fruit Texture Analyzer, GS-14, Güss Manufacturing). A small section of skin was 149 removed using a vegetable peeler, and each fruit was placed on the penetrometer platform so that 150 the piston entered the middle of the apple where the skin had been removed. Data were 151 automatically recorded into a spreadsheet. 152 Due to the number of trees, harvesting the ABC orchard often lasted more than one day. 153 Therefore, differences in harvest date within a week were simply due to the time required to 154 8 harvest, and the harvest date for each tree was recorded as the Monday of that week. We used the 155 "lmer" function in the R package lme4 (Bates et al., 2015) to fit a REML model for harvest and 156 firmness data. Next, we calculated the least squares mean using the "lsmeans" function in the 157 lsmeans R package (Lenth, 2016), resulting in one value per unique accession. We fit this model 158 for 866 unique accessions with harvest dates and 863 accessions with firmness measurements. 159 160

High Resolution Melting (HRM) Genotyping 161
DNA was extracted using silica columns from fresh leaf tissue collected from the ABC 162 orchard, quantified using PicoGreen (Thermo) and normalized to a concentration of 20 ng µL -1 . 163 Genotyping was conducted using PCR and high resolution melting (HRM) on a LightScanner 164 HR384 (BioFire). Primers are listed in Supplemental Table S1.

Marker-phenotype associations 182
We evaluated the ability of markers at NAC18.1, ACO1, ACS1, and PG1 to predict both 183 harvest date and firmness of accessions in the ABC orchard. The number of accessions with 184 phenotype and genotype information ranged from 754 to 852 depending on the trait/marker 185 combination. The association test was conducted using Spearman's rank correlation test. We 186 visualized results using the "geom_boxplot" function in ggplot2 in R (Wickham, 2016). Lastly, 187 we tested which mode of inheritance best fit the data using SNPstats (Solé et al., 2006), and 188  Table S1) and 194 Phusion® High-Fidelity PCR Master Mix with HF Buffer (NEB). PCR product size and purity 195 was confirmed by agarose gel electrophoresis, and the remaining product was purified using 196  Individual colonies were selected for complete sequencing of the cloned amplicon using 199 the primers NAC18F2, NAC18F3, NAC18F4, NAC18R1, and NAC18R2 (Supplemental Table  200 S1). For accessions homozygous for the D5Y SNP, the NAC18.1 amplicon from a single clone 201 was sequenced. For heterozygous accessions, two clones representing each D5Y allele were 202 selected based on partial sequencing of the D5Y region, followed by complete sequencing of the 203 2.3 kb amplicon, as described above. The nucleotide sequences were aligned by MUSCLE 204

Evaluation of markers for firmness and harvest date in apple 242
We sought to re-evaluate the ability of D5Y and other published markers to predict fruit 243 firmness and harvest date using new phenotype and genotype data: specifically, we aimed to 244 address some experimental limitations of our previous work. First, the GWAS that identified 245 D5Y, and failed to find associations for previously published markers, made use of only 8,000 12 SNPs (Migicovsky et al., 2016). Given the rapid LD decay in apple, our low SNP density was 247 insufficient to conclude that we had exhaustively searched the apple genome for loci involved in 248 firmness and harvest date (Migicovsky et al., 2016). Secondly, the historic phenotype data we 249 used were imprecise. For example, firmness was recorded as either "firm" or "soft" rather than 250 as a biomechanical measurement, and harvest date was recorded as simply "early", "mid", or 251 "late". These factors would have limited the detection of additional loci that underlie firmness 252 and/or harvest date. 253 In order to directly compare the NAC18.1 D5Y marker that we identified with previously 254 published markers, we developed HRM genotyping assays for the D5Y SNP, as well as markers 255 at ACS1, ACO1, and PG1. We used these assays to genotype accessions from the ABC orchard 256 (Supplementary Table S2). While this orchard largely consists of clones of individuals present in 257 the population used in our previous GWAS (i.e. the USDA germplasm collection), it also 258 contains a number of additional accessions. More importantly, for the purposes of phenotyping, 259 the ABC orchard is planted in an incomplete block design that allows for modeling of location 260 effects using REML methods (Migicovsky et al., 2018). Using this approach, we generated 261 firmness and harvest date measures for over 800 unique accessions for which HRM genotype 262 data were available. 263 We evaluated the mode of inheritance for each genetic marker and phenotype 264 combination and found that, while the effects of the markers at ACS1, ACO1, and PG1 were all 265 dominant, the D5Y marker at NAC18.1 had a codominant effect. As a result, accessions that 266 were heterozygous for D5Y had firmer fruit (an increase of 1.2 kg cm -2 ) and a later harvest date 267 (20.13 days) than those that were homozygous for the 'soft' A allele. The combined effect of two 268 C alleles was notably stronger than the effect of a single C allele, increasing firmness by 2.24 kg 269 13 cm -2 and harvest date by a month (30.52 days), when compared to AA genotypes. The difference 270 between firm/late and soft/early genotypes for D5Y was at least 4 times higher for firmness and 271 3 times higher for harvest date than for the markers at ACS1, ACO1, and PG1 (Figure 1). To evaluate the distribution of firmness/harvest date markers across commercial 280 cultivars, we present the genotypes for all markers genotyped here in nine of the top ten apple 281 cultivars sold in the USA (Table 1)  Since the D5Y SNP is predicted to result in an amino acid change in the NAC18.1 306 protein and LD decay is generally high in apple (r 2 decayed to < 0.2 within 100 bp), it is possible 307 that D5Y is causative for the early/firm phenotype (Migicovsky et al., 2016). However, GBS 308 sequence data are too sparse to rule out additional polymorphisms in LD with D5Y contributing 309 to the phenotype. To better understand allelic diversity of NAC18.1, we selected 6 cultivars each 310 of the D5Y A/A, A/C, and C/C genotypes, based on our HRM genotyping data. For homozygous 16 cultivars, single alleles were sequenced, while for the heterozygous samples we resequenced 312 both alleles, resulting in NAC18.1 sequences from 24 haplotypes. These were compared to two 313 reference sequences from the GDDH13 v1.1 apple genome (Daccord et al., 2017): 314 Md03g1222600 and Md11g1239900, representing NAC18.1 and its closest paralog, respectively. 315 In addition to confirming the expected D5Y genotype in all individuals, a number of 316 additional SNPs and indels were revealed within both coding and non-coding regions of the 317 NAC18.1 allele. Multiple sequence alignment and subsequent phylogenetic analysis indicated 318 two major clades corresponding to the D5Y A and C genotypes (Figure 2A). SNPs and small 319 indels were observed throughout the sequenced region. Notably, in addition to the D5Y amino 320 acid change, several additional amino acid changes were observed between sequences from D5Y 321 genotypes and the reference sequence of NAC18.1. For example, near the site of the D5Y 322 polymorphism, all "A" haplotypes also had a 12 nucleotide insertion that introduced the amino 323 acid sequence QPQP ( Figure 2B).  19 Expression of all three marker genes was enhanced in the NAC18.1 transgenic lines 371 relative to the nor mutant control, although not to the same extent as observed in WT ripe fruit 372 ( Figure 3B-D). In contrast to the consistent level of NAC18 expression observed in each line, the 373 marker genes were more variable in their expression level between lines. A similar pattern was 374 observed for all marker genes, with the NAC18.1 C #9 line showing the smallest induction of 375 marker gene expression relative to nor. In the case of PG2 and ACS2, the difference in 376 expression in NAC18.1 C #9 was not statistically significant relative to nor (p = 0.09 and 0.18, 377 respectively). Consistent with these results, fruit from this line also exhibited the lowest amount 378 of color development ( Figure 3B). Taken together, these results indicate that a canonical ripening 379 program can be induced in the tomato nor mutant through the heterologous expression of either 380 apple NAC18.1 allele. 381