Dynamic QTL analysis for developmental behavior of cell wall components and forage digestibility in maize (Zea mays L.)

Cell wall architecture plays a key role in stalk strength and forage digestibility. Lignin, cellulose and hemicellulose are the three main components of the plant cell wall and 20 can impact stalk quality by affecting cell wall structure and strength. To explore cell 21 wall development during secondary cell wall lignification in maize stalks, conventional 22 and conditional genetic mappings was used to identify the dynamic quantitative trait 23 locus (QTL) for cell wall components and digestibility traits in five growth stages after 24 silking. found not associated with stalk strength. A simultaneous improvement in forage 40 digestibility and lodging resistance can be achieved by pyramiding multiple effective 41 QTL identified in the present study.


Background 45
Maize (Zea mays L.) can be bio-refined to provide sustainable bioproducts and 46 bioenergy [1]. More importantly, the stover of this staple crop is usually treated as a 47

Phenotyping methods 117
The silking date was recorded when 50% of the silks emerged from all plants in each 118 line. The 2nd-5th internodes above the ground from six plants of each line were 119 collected using garden shears at 10, 20, 30, 40, and 50 DAS, respectively. All samples 120 were immediately enzyme-deactivated at 105°C for 30 min and air-dried for 10-14 days. 121 Dried stalk samples were ground using an automatic hammer mill herb grinder and 122 filtered through a screen with a mesh size of 0.1 mm. Before measurements, the stalk 123 samples were dried at 45°C for 48 h to exclude the influence of moisture. CEL, ADL, 124 ADF, NDF, and IVDMD were estimated using near-infrared reflectance spectroscopy. 125 The samples were scanned using a near-infrared reflectance spectrophotometer (Vector 126 22/N, Bruker Optik, Ettlingen, Germany). A modified partial least squares approach 127 implemented in OPUS 6.0 Bruker software was used to fit the calibration equations. 128 The coefficients of determination for cross-validation ( R 2 ) ranged from 90.2% 129 (IVDMD) to 94.0% (CEL), and the coefficients of external validation (R 2 ) were 92.7% 130 for ADL, 96.7% for CEL, 94.6% for ADF, 96.5% for NDF, and 91.2% for IVDMD. 131

Phenotypic data analyses 132
A mixed linear model implemented with the "lmer" function in the "lme4" package in 133 R version 3.6.0 (R development Core Team, 2019) was fitted to calculate the best linear 134 unbiased prediction (BLUP) value for each line: = + + + , where 135 represents the phenotype of the "i"th line, is the grand mean value of the target trait 136 in all environments, represents the genetic effect, is the environmental effect 137 (replications in each environment were also treated as environmental effects in the 138 BLUP mixed model), and is the random error. The grand mean was fitted as a fixed 139 effect, and genotype and environment were considered as random effects. The A standard analysis of variance was conducted using the base "aov" function in R. 143 The model used for the analysis of variance was = + + ( ) + + 144 , where is the environmental effect of the "l"th environment, ( ) is 145 the effect of replications within environments, represents the genetic effect of the 146 "i"th line, ( ) is the interaction effect between genetic and environment effects, and 147 is the residual error. All the effects were considered as random effects. Broad-sense 148 heritability was calculated as 2 = 2 /( 2 + 2 / + 2 / ), where 2 represents 149 the genetic variance, 2 is the variance of interaction between the genotype and 150 environment, 2 is the residual error variance item, and and r are the number of 151 environments and replications in each environment, respectively. The 95% confidence 152 intervals of 2 were calculated following the method of Knapp et al. [32]. 153

Genotyping and genetic map construction 154
Leaf tissues were collected from all 215 RILs and their parental lines and freeze-dried 155 at -60°C. Genomic DNA was extracted using the modified CTAB method [33] and used 156 for genotyping with the maize 55K single-nucleotide polymorphism (SNP) array [34]. 157 The quality of each SNP was manually controlled as described by Yan

Phenotypic variation and correlation between traits 184
Each line of the RIL population was sampled at five stages, including 10-50 DAS 185 (denoted as stages I, II, III, IV, and V), and evaluated for CEL, ADL, ADF, NDF, and  186 IVDMD. For convenience, the phenotypic value of each trait at each stage was denoted 187 as "Trait_Stage." For example, "ADF_I" represents the ADF measured at stage I. The 188 cell wall components and digestibility traits were evaluated in five stages and two 189 environments in the present study. The phenotypic values of ADF, NDF, CEL, and ADL 190 in Zheng58 were slightly lower than those in HD568 at stage I, nearly equal to those in 191 HD568 at stage II, and surprisingly surpassed those in HD568 at stages III and IV. At 192 the final stage, the values of these traits in Zheng58 showed a significant decrease, and 193 they were lower than those in HD568. For the RIL population, the average IVDMD  (Table 1). 206

Unconditional QTL mapping 207
After 1,000 permutation tests, the empirical threshold LOD values for the genome-wide 208 significance (p < 0.05) were determined, ranging from 3.8 to 4.1 for various stages and 209 traits. A total of 72 unconditionally detected QTL for cell wall components and 210 digestibility traits at five sampling stages were identified on nine chromosomes. The 211 percentage of trait variation explained by each QTL ranged from 3.48% to 24.04%. 212 Approximately 22.2% (16/72) QTL accounted for more than 10% of the phenotypic 213

variation. 214
For ADF, 14 QTL were detected on chromosomes 1, 2, 4, 7, 8, and 9 at the five 215 stages, individually explaining 5.20-13.48% of the phenotypic variance (Table 2). 216 Among these QTL, two major effect QTL were detected on chromosomes 2 and 9. On 217 chromosome 1, a QTL located at 201.5-207.8 cM was detected at multiple stages 218 (stages I, II, and IV) and was responsible for 6.67%, 7.01%, and 6.25% of the 5.56% of the phenotypic variance at stages I, II, and III, respectively. Moreover, two 222 overlapping QTL were detected at 62-66 cM and 42-45 cM on chromosomes 7 and 9, 223 respectively. 224 For NDF, 16 QTL were detected on chromosomes 1, 2, 7, 8, and 9. Each QTL 225 could explain 3.48-13.97% of the total phenotypic variance at different stages. Thirteen Twelve unconditional QTL for ADL were identified during the five development 230 stages. Among these, 11 QTL were distributed on chromosomes 1, 2, and 6. Similar to 231 that for ADF, the QTL located at 201.3-207.8 cM on chromosome 1 were identified for 232 ADL at stages I and II. During stages I, II, and V, two consensus QTL related to ADL 233 were found at 140.7-141.2 cM on chromosome 2 and 20.7-21.3 cM on chromosome 6. 234 Moreover, the overlapped QTL on chromosome 2 contributed to 24.04%, 15.18%, and 235 16.11% of the phenotypic variance at stages I, II, and V, respectively, which seems to 236 be a stable and major effect QTL for ADL. 237 During the five stages, 11 unconditional QTL for CEL were detected. These QTL 238 were distributed on chromosomes 1, 2, 6, and 7 (Table 2, Figure 3). Each QTL could 239 explain 5.12-18.97% of the total variance. The QTL located at 91.1-92.8 cM on 240 chromosome 2 was repeatedly detected at stages I, II, III, and IV and contributed to 241 5.42-18.97% of the total phenotypic variance. 242 For IVDMD, 19 QTL were detected on eight chromosomes. Among these, 12 243 overlapping QTL were integrated into four consensus QTL. Among these overlapping 244 QTL, four were identified repeatedly at 202.7-207.8 cM on chromosome 1. These QTL 245 contributed to 6.66%, 12.99%, 6.76%, and 5.81% of the total phenotypic variance at 246 stages II, III, IV, and V, respectively. In addition, the QTL located at approximately 247 92.0-92.8 cM on chromosome 2 was correlated with IVDMD at stages I, II, and III. On IVDMD at stages I, II, and V. Moreover, another overlapping QTL for IVDMD was 250 detected on chromosome 10, which was located at 49.5-50.2 cM and contributed to 251 5.07% and 8.16% of the total phenotypic variance at stages I and V, respectively. 252

Conditional QTL mapping 253
In total, 26 conditional QTL were identified at five stages for cell wall components and 254 digestibility traits (Table 3) Five ADF-related conditional QTL were identified as distributed on chromosomes 259 1, 2, 3, and 6. When the ADF at stage V was conditioned on the ADF at stage IV, 260 conqADF1b and conqADF2 were detected on chromosomes 1 and 2; these QTL were 261 responsible for 10.15% and 11.99% of the total variance, respectively. At the IV|III 262 stage, conqADF1a and conqADF3 were found to contribute to 1.02% and 3.89% of the 263 total variance, respectively. On chromosome 6, conqADF6 was detected with a 264 contribution of 5.40% when the ADF at stage III was conditioned on the ADF at stage 265

II. 266
Four conditional QTL for ADL were detected on chromosomes 1, 2, and 10. QTL 267 conqADL1a and conqADL1b overlapped the QTL located on chromosome 1 and 268 contributed to 5.34% and 7.98% of the total variance for ADL at the IV|III and V|IV 269 stages, respectively. In addition to that on chromosome 1, a conditional QTL on 270 chromosome 2, conqADL2, was detected when the ADL at stage V was conditioned on 271 the ADL at stage IV. This major effect QTL could explain 13.41% of the total variance. 272 Moreover, at stage III|II, conqADL10 was identified to be responsible for 4.84% of the 273 total variance. 274 For CEL, five conditional QTL were found at stages III|II, IV|III, and V|IV. At 275 stages IV|III and V|IV, QTL conqCEL1a and conqCEL1b were located at adjacent positions on chromosome 1 and contributed to 3.93% and 10.70% to the total variance, 277 respectively. Another major conditional QTL at the V|IV stage located on chromosome 278 2 was conqCEL2 and contributed to 14.95% of the total variance. When the CEL at 279 stage III was conditioned on the CEL at stage II, conqCEL5 and conqCEL7 were 280 identified on chromosomes 5 and 7. 281 Conditional QTL mapping for IVDMD revealed six QTL on chromosomes 1, 2, 6, 282 9, and 10. These conditional QTL were identified at the III|II and V|IV stages. When 283 the IVDMD at stage III was conditioned on the IVDMD at stage II, QTL 284 conqIVDMD1a, conqIVDMD2, and conqIVDMD6 were responsible for 6.39%, 6.06%, 285 and 5.21% of the total variance, respectively. In addition, QTL conqIVDMD1b, 286 conqIVDMD9, and conqIVDMD10 identified at stage V|IV, could explain 8.68%, 287 9.08%, and 7.67% of the total variance in the IVDMD. 288 Similar to that for the IVDMD, six conditional QTL were identified for NDF. All 289 these QTL were also found at stages III|II and V|IV. Moreover, three conditional QTL 290 for NDF detected at stage III|II overlapped with the conditional QTL for IVDMD 291 detected at the same stage (conqIVDMD1a, conqIVDMD2, and conqIVDMD6). At 292 stage V|IV, three conditional QTL conqNDF1b, conqNDF2b, and conqNDF9 293 contributed to 9.00%, 6.91%, and 10.12% of the total variance in the NDF. conqNDF1b 294 was found to be overlapped with conqIVDMD1b, which was identified at stage V|IV 295 for IVDMD. 296

297
The area of maize in China has been expanding annually. The livestock industry 298 requires large amounts of silage maize owing to the increase in demand for meat and 299 milk. However, common maize is still the main variety promoted in China because of 300 the policy of self-sufficiency in food production. Although silage maize has a large 301 amount of biomass, higher than that of common maize, its starch and dry matter content 302 is lower than that of common maize. Therefore, silage made of grain and forage maize varieties is more popular in China. The critical issue is that the choice of harvesting 304 time is a trade-off among yield, feeding value, and silage quality. 305 Lodging is a critical phenomenon for maize production, causing severe yield 306 reduction during the reproductive stage. Stronger cell walls can provide more powerful 307 mechanical support to avoid lodging. By contrast, cell wall digestibility, which is 308 strongly affected by cell wall lignification, is a vital characteristic for the nutritional 309 value of forage maize [42]. Cell wall lignification is a dynamic and complex process 310 that occurs throughout maize growth. In the present study, we used a maize RIL 311 population to investigate the dynamic changes in cell wall component and digestibility 312 traits at five stages from silking to harvest. For each trait, no apparent differences were 313 observed among the beginning stages. ADF, NDF, ADL, and CEL showed a massive 314 increase at stages III and IV and exhibited a slight decrease at stage V. Conversely, 315 IVDMD showed a large decrease at stages III and IV and exhibited a small increase at 316 stage V. This inverse trend may be associated with significant negative correlations 317 between these traits and IVDMD. These results revealed that the optimal harvest time 318 for grain-forage maize was approximately 50 DAS. At this time, the plant has reached 319 physiological maturity and the grain is in full dent stage. After harvesting for grain, the 320 plant can also be harvested for forage or stover to feed animals. Alternatively, forage 321 maize could be harvested between stages II and III, which is roughly in the dough stage. performed dynamic QTL analysis during five developmental stages after silking. 335 As there was no significant difference among the phenotypic values of stages I and 336 II, no conditional QTL was found for any of the traits when stage II was conditioned on 337 stage I. Moreover, several overlapping unconditional QTL were detected between these 338 two stages. In general, unconditional QTL hotspots were located on chromosomes 1, 2, 339 7, 8, and 9. The co-localized QTL on chromosome 1 was located in bin 1.08, which 340 corresponded to the 215-250 Mb physical region in the maize genome (Version 5.60) 341 [45]. This genomic region was also detected as a meta-QTL hotspot for cell wall content 342 and digestibility [31]. Moreover, this locus was shown to be associated with stalk 343 strength in another study [38] in which QTL mapping was performed for rind 344 penetrometer resistance (RPR) using the same RIL population that was used in the 345 present study. We found that another pleiotropic QTL for RPR was associated with ADF, 346 NDF, and IVDMD. This pleiotropic QTL was localized in bin 8.05 and explained 6.61-347 14.06% of the phenotypic variance for each trait. Therefore, this major effect QTL 348 should be the target for fine mapping and gene cloning. 349 In addition to the pleiotropic QTL mentioned above, we detected several 350 pleiotropic QTL for cell wall and digestibility; however, these were not associated with 351 RPR. On chromosome 2, a co-localized QTL in bin 2.07 was found to be related to CEL 352 and ADL at various stages. Another overlapped QTL in bin 2.04 was repeatedly 353 detected for all four traits except ADL. Thus, we concluded that a cellulose synthesis 354 gene is responsible for the genetic variation at this locus. Another overlapping QTL 355 region was found on chromosome 7. At stage III, this QTL located in bin 7.03 was 356 detected for all five traits. A hotspot QTL on chromosome 9 was associated with ADF, 357 digestibility QTL hotspots should be the focus of attention for further cloning of 360 candidate genes and favorable alleles. The genes underlying these QTL should be 361 responsible for secondary cell wall digestibility but not affect cell wall rigidity. In this 362 context, it is possible to improve the digestibility of silage through molecular marker-363 assisted modification of these QTL without reducing stalk strength. Furthermore, a 364 simultaneous improvement in forage digestibility and lodging resistance can be 365 achieved by pyramiding multiple effective genes underlying all the pleiotropic QTL. 366 Most conditional QTL were mainly detected for stages III|II and V|IV owing to the 367 significant phenotypic differences between the two adjacent stages. Six out of 26 368 conditional QTL explained more than 10% of the phenotypic variation for each trait.

Ethics approval and consent to participate 515
Not applicable 516

Consent for publication 517
Not applicable 518

Availability of data and material 519
The datasets supporting the conclusions of this article are included within the article 520 and its additional files. 521

Competing interests 522
The authors declare that they have no competing interests. 523    Table 1. Phenotypic values of cell wall components and digestibility traits at different developmental stages in a maize recombinant inbred 551 line (RIL) population 552   Figure 1 Changes in cell wall components and digestibility traits after silking in a recombinant inbred line (RIL) population. ADF, acid detergent ber; ADL, acid detergent lignin; CEL, cellulose; NDF, neutral detergent ber; IVDMD, in vitro dry matter digestibility. I-V represent 10-50 days after silking, respectively.