Genomic Prediction of Drought Tolerance during Seedling Stage in Maize using Low-Cost Markers

24 Drought tolerance in maize is a complex and polygenic trait, especially in the seedling stage. 25 In plant breeding, such traits can be improved by genomic selection (GS), which has become 26 a practical and effective tool. In the present study, a natural maize population named 27 Northeast China core population (NCCP) consisting of 379 inbred lines were genotyped with 28 diversity arrays technology (DArT) and genotyping-by-sequencing (GBS) platforms. Target 29 traits of seedling emergence rate (ER), seedling plant height (SPH), and grain yield (GY) 30 were evaluated under two natural drought environments in northeast China. adequate 31 genetic variants have been found for genomic selection, they are not stable enough between 32


23
Abstract 24 Drought tolerance in maize is a complex and polygenic trait, especially in the seedling stage. 25 In plant breeding, such traits can be improved by genomic selection (GS), which has become 26 a practical and effective tool. In the present study, a natural maize population named 27 Northeast China core population (NCCP) consisting of 379 inbred lines were genotyped with 28 diversity arrays technology (DArT) and genotyping-by-sequencing (GBS) platforms. Target  29 traits of seedling emergence rate (ER), seedling plant height (SPH), and grain yield (GY) 30 were evaluated under two natural drought environments in northeast China. adequate 31 genetic variants have been found for genomic selection, they are not stable enough between 32 two years. Similarly, the heritability of the three traits is not stable enough, and the 33 heritabilities in 2019 (0.88, 0.82, 0.85 for ER, SPH, GY) are higher than that in 2020 (0.65, 34 0.53, 0.33) and cross-two-year (0.32, 0.26, 0.33). The current research obtained two kinds of 35 marker sets: the SilicoDArT markers were from DArT-seq, and SNPs were from the GBS and 36 DArT-seq. In total, a number of 11,865 SilicoDArT, 7,837 DArT's SNPs, and 91,003 GBS 37 SNPs were used for analysis after quality control. The results of phylogenetic trees showed 38 that the population was rich in consanguinity. Genomic prediction results showed that the 39 average prediction accuracies estimated using the DArT SNP dataset under the 2-fold 40 cross-validation scheme were 0.27, 0.19, and 0.33, for ER, SPH, and GY, respectively. The 41 result of SilicoDArT is close to the SNPs from DArT-seq, those were 0.26, 0.22, and 0.33. 42 For SPH, the prediction accuracies using SilicoDArT were more than ones using DArT SNP, 43 In some cases, alignment to the reference genome results in a loss to the prediction. The trait 44 with lower heritability can improve the prediction accuracy using filtering of linkage 45 disequilibrium. For the same trait, the prediction accuracy estimated with two types of DArT 46 markers was consistently higher than those estimated with the GBS SNPs under the same 47 genotyping cost. Our results show the prediction accuracy has been improved in some cases 48 of controlling population structure and marker quality, even when the density of the marker is 49 reduced. In the initial maize breeding cycle, Silicodart markers can obtain higher prediction 50 accuracy with a lower cost. However, higher marker density platforms i.e. GBS may play a 51 role in the following breeding cycle for the long term. The natural drought experimental 52 station can reduce the difficulty of phenotypic identification in a water-scarce environment. 53 The accumulation of more yearly data will help to stabilize the heritability and improve 54 predictive accuracy in maize breeding. The experimental design and model for drought 55 resistance also need to be further developed.

59
Maize (Zea mays L.) is an important source of food, animal feed, and energy because of its 60 high yield potential (Jiang et al. 2018). Consequently, poor maize harvests, usually caused by 61 environmental stresses, will have a tremendously detrimental impact on human life. Seed 62 germination and seedling establishments are the initial stages of maize life and influence 63 seedling emergence rate, uniformity, and robustness, which determine maize yield potential 64 (Tian et al. 2014). However, these stages of maize are frequently subjected to severe drought 65 stress worldwide due to increasingly global warmer weather, which decreases the emergence 66 rate and retards seedling growth. Thus, it is necessary to breed more drought-tolerant maize 67

hybrids. 68
Selection of elite drought-resistant maize inbred lines is the basis of obtaining 2017) and grain quality (Guo et al. 2020). Compared with the widely used marker-assisted 84 selection (MAS) technology, the GS strategy is advantageous in that it does not need to 85 identify significant markers, genes, or quantitative trait locis (QTL) related to traits in 86 advance, and can be used to predict the performance of populations without genetic structure 87 analysis (Bernardo 2016   The process from maize seed germination to seedling establishment involves a variety of 129 complicated metabolic transformations, physiological activities, and cellular and tissue 130 differentiation, which are regulated by many genomic loci, including major and micro effect 131 loci under a drought environment. An abundant variation of these loci takes place among 132 different maize germplasms. In this study, we used the combined genotyping data of natural 133 maize population obtained by DArT and GBS sequencing platforms, as well as the phenotypic 134 data of emergence rate (ER), seedling plant height (SPH) for seedling stage, and grain yield 135 (GY) for the mature stage under drought environment, to conduct GS of drought resistance. 136 The main objectives were to: (1) assess the phenotypic variation in the NCCP under the 137 natural drought environment; (2) estimate the genomic prediction accuracy for ER, SPH, and 138 GY traits using DArT and GBS markers; (3) try to improve the prediction accuracies by 139 increasing the marker quality; and (4) control population structure for improve the prediction 140 accuracies. 141 October 2020. In the 2019 maize growing season, the average temperature was 21 °C and the 151 average rainfall was 4 mm. In the 2020 maize growing season in, the average temperature was 152 21°C and the average rainfall was 3mm. During the maize seedling stage, the highest 153 temperature across two years was 35°C, and the precipitation was less than 20 mm, causing 154 serious drought stress in the seedling stage, where three leaveswere maintained in the whole 155 growing season (Fig.1). Trials in our research did not take irrigation measures, but the seeds 156 were planted before precipitation according to the weather forecast for germination. A 157 completely randomized block design with three replications was applied in each evaluated 158 environment. The inbred lines were sown using a seed for each hole in a single plot, 3 m long, 159

Materials and Methods
with plant spacing of 10 cm and row spacing of 60 cm. Target traits of emergence rate (ER), 160 seedling plant height (SPH) and grain yield (GY) were evaluated to represent the drought 161 tolerance of the tested plant materials. The target traits of ER and SPH traits were measured at 162 the seedling stage, i.e. 20 days after planting. While GY was measured at the maturity stage, 163 the moisture content readings were made by LDS-1G Grain Moisture Meter. ER was 164 measured as the proportion of surviving plant number to number planted. SPH was measured 165 as the distance from the plant base to the highest of the seedling plant. The average dry weight 166 of five plants was regarded as the GY score for entries. 167 to estimate the variance components. 173

Phenotypic data analysis and heritability estimation
Where is the trait of interest, μ is the overall mean, , , and × 175 are the effects of the i-th genotype, j-th year, and i-th genotype by j-th year interation, 176 respectively.
is the effect of k-th replication. is the residual effect of the i-th 177 genotype, j-th year, k-th replication. Genotype is considered as the fixed effect, whereas all 178 other terms are declared as the random effects. Years with heritability below 0.05 were 179 excluded from the across location analysis. 180 Broad-sense heritability (H 2 ) based on the entry means within trials was estimated as follows:  (Table S1). 201 Genotyping of all the inbred lines with DArT-seq method was carried out at the SAGA 202 is different from GBS, which does not rely on the reference genome information, but mainly 215 depends on the DArT company's sequencing the sample library data. All reads were aligned 216 by sequence analysis based on maize tags of metagenome representation. Two types of data 217 were generated by DArT-seq: SilicoDArT and SNP. The former is the presence/absence 218 variation of the tag sequences, while the last one, SNPs, have been widely used in genetics 219 and molecular marker-assisted selection breeding fields. The type of SNP data was obtained 220 by B73 RefGen_v4 using the BLAST tool, and some of the fragments in SilicoDArT cannot 221 be aligned to the chromosomes and abandoned. For this reason, the number of SilicoDArT is 222 more than SNP. SilicoDArT of 62,794 were acquired, and 39,659 SNPs were located on 223 maize chromosomes 1-10; however, 547 SNPs were anchored on the zero chromosome (Table  224 S1). 225 In both DArT and GBS datasets, TASSEL 5 (Bradbury et al. 2007) was used to filter out 226 markers with minor allele frequency (MAF) < 0.05 and a missing rate > 20%. The samples 227 with a missing rate of more than 20% were rejected. The number of lines in DArT and GBS is 228 379 and 378, respectively. In total, the markers are 11,865 SilicoDArT and 7,837 SNPs in the 229 DArT platform and 91,003 SNPs in the GBS dataset (Table 1)

Genomic prediction model 243
The RR-BLUP model assumes homogenous variance of all markers and shrinks all marker 244 effects equally to zero. RR-BLUP is equivalent to BLUP and uses the realized relationship 245 matrix estimated from the markers. The genomic prediction analysis was performed with the 246 ridge regression BLUP (rrBLUP) package (Endelman, 2011). The mixed model is described 247 where y is the vector (n × 1) of observations, X is the vector (n × 1) of individuals and β 250 is the fixed effects; ε is the vector (n × 1) of independently random residuals with assumed 251 distribution N (0, Iσε 2 ), Z is the design matrix (n × m) for random effects and u is the vector of 252 random effects with u ∼ N (0, Kσu 2 ) and K being an identity matrix in this case. In addition, n

260
To evaluate the impact of TPS on the accuracy of genomic prediction, the random 261 cross-validation scheme was used to randomly generate different training and prediction sets 262 and assess the prediction accuracy within each TPS. The inbred lines ranged from 10% to 263 90%, with an interval of 10%, were extracted from the total lines and were used to establish 264 the training population to predict the remaining inbred lines. The genomic prediction was 265 repeated 100 times for each TPS. For the same trait, the prediction accuracies estimated from 266 the DArT and GBS datasets were compared. 267 SNP sets in GBS/DArT dataset were filtered by R language, and the SNPs were removed 284 when r 2 was greater than or equal to specific threshold. Genomic prediction accuracy across 285 100 times repetition was calculated using 2-fold cross-validation. 286  We performed a principal component analysis (PCA) to assess the prediction accuracy 305 affected by genetic structures. All the inbred lines used in the present study were divided into 306 two subgroups because it is necessary to ensure that the subsets have enough material for 307 genomic prediction. The prcomp function from the stats package and plot function from the 308 base package were executed to draw PCA plots under the R environment. Genomic prediction 309 accuracy across 100 times repetition was estimated within each subgroup using 2-fold 310 cross-validation. 311

Cross-validation for assessment 312
In this study, the size of the training population with the highest average prediction accuracy 313 and the lowest standard deviation is about 50%, which is equivalent to 2-fold cross-validation. 314 A 2-fold cross-validation analysis was implemented for all traits, i.e., 50% of individuals 315 randomly selected for the training population and 50% in the testing population. Additionally, 316 different training population sizes were used to evaluate the performance of population size. 317 For prediction, regression models were evaluated by the prediction accuracy, a Pearson's 318 correlation between the GEBVs and the phenotypic data in the testing population 319

Phenotypic data and correlation analysis 323
The phenotypic mean and standard deviation of the three traits in 2019, 2020 and BLUEs 324 across two years are shown in Table 2 (Table 2).  Table 2). The phenotypic data analysis based on BLUEs revealed 339 that ER, SPH, and GY showed normal distribution or skewed normal distribution (Fig.2). All 340 three traits showed significant positive correlations, and the correlation coefficients are 0.36, 341 0.35, and 0.25 (P < 0.01) between ER and SPH, ER and GY, and SPH and GY, respectively. 342

Quality control and marker distribution on genotypic data
343 Two kinds of marker sets were obtained in the current research. The SilicoDArT markers 344 were from DArT-seq, and SNPs were from the GBS and DArT-seq. The plots of MAF and 345 missing rate distribution before filtering for three marker sets are shown in Fig. 3

. The 346
SilicoDArT markers had an average MAF of 0.10, and SNPs' average MAF of 0.11 of filtering is 0.11. Before filtering, the highest peak is found at the MAF intervals of 0 to 0.05 349 for three marker sets. However, we are not sure whether the small MAF values are from the 350 actual situation or sequencing errors, so they can not be used for further analysis. The MAF 351 distribution of other intervals (0.05 to 0.5) is relatively uniform. Average minor allele 352 frequency after filtering are 0.24, 0.24, and 0.23 for SilicoDArT, SNPs from DArT-seq, and 353 SNPs from GBS, respectively (Table 1). Consistent MAFs were found on both genotyping 354 platforms (Fig.3 a,b,c), potentially indicating that sequencing quality is reliable. It is worth 355 noting that higher quality marker sets were used in the present study. 356

Results of the distribution of missing rates showed a decreasing trend in SilicoDArT and 357
SNPs from DArT-seq markers, and the missing rate is mainly concentrated in the range of 358 0-10%. On the other hand, the highest peak, accounting for more than 40 percent, was found 359 in the GBS markers for the missing rate, which appears at the interval of 90% to 100%. For 360 each interval of 0 to 10%, 10% to 20%, 20% to 30%, 30% to 40%, and 40% to 50% in GBS 361 markers, the missing ratio accounts for about 10%. There were almost zero accounts in the 362 range of 50% to 90% (Fig.3 d,e,f). The irregular missing rates results of GBS markers were 363 not apparent and thus may not be suitable for alignment to a single reference genome in the 364 SNP calling process due to the diversity of materials. 365 High-quality marker sets were obtained after filtering had been done, and the average 366 missing rates are 0.070, 0.073, and 0.068 for SilicoDArT, SNPs from DArT-seq, and GBS 367 markers (Table 1). Heterozygosity is usually a critical index for marker filtering, especially 368 for inbred lines. In this research, the proportion heterozygous of the three marker datasets was 369 very low before and after filtering, which means our material is homozygous and implies the 370 sequencing is accurate. Specifically, the averages of heterozygosities were less than 0.01 371 before filtering for the three marker sets. And after filtering, the heterozygosity rates were 372 increased, but all were still less than 0.03 (Table S1 and Table 1). 373 When comparing pre-filtering and post-filtering, the number of markers decreased by 374 81%, 80%, and 88% for SilicoDArT and SNPs from DArT-seq and SNPs from GBS, 375 respectively (Table S1 and Table 1). The missing rate is mainly responsible for the reduction 376 of markers, and those were 0.25, 0.24, and 0.56 for the SilicoDArT, SNPs of DArT, and SNPs 377 of GBS shown in Table S1. The high missing rate indicates that the diversity of the population 378 is rich, but a large number of markers will be lost when the marker set is summarized. 379 The SNP marker sets from DArT-seq and GBS cover almost the entire maize genome 380 (Fig.4). Generally, the density of markers at both ends of the chromosomes was high and 381 decreased toward the centromere position. The SNP marker sets from DArT-seq and GBS 382 platforms were enriched at both ends of chromosomes 5 and 8. In addition, markers from 383 GBS had a more uniform distribution than ones in DArT, whereas DArT's SNPs had higher 384 density on the long arms of chromosomes 8 and 9. Uneven distribution may potentially lead 385 to a decrease in prediction accuracy when implementing a genomic prediction. 386

387
The results of the phylogenetic tree are shown in Fig.5 using SNPs of DArT-seq and GBS. were classified as unknown. DArT SNPs produced a better phylogenetic tree, and the markers 395 more clearly separate the Reid, Lvxi, and Huanglvxi from NCCP. However, GBS markers did 396 not clearly separate the materials of heterotic groups in the present study. Some of the 397 materials from local breeders should belong to the mixed group or others, but they were 398 unable to determine sources for a long time, so they were classified into the unknown group. 399 Although some materials have source information, we are not sure that the information is 400 correct when the materials are collected from third parties. Therefore, it isn't easy to obtain a 401 perfect phylogenetic tree in practical breeding work. In the pedigree, we have conflicts with a 402 phylogenetic tree based on markers, making it hard to decide which one we can believe. At 403 least, the source of inbred lines of NCCP is rich based on the result of DArT markers. 404 However, it may cause trouble for the genomic prediction. 405

The genomic prediction accuracy estimated from the DArT and GBS
406 markers 407 Genomic prediction accuracies were estimated for all the three target traits with the 408 BLUE values and the SNP datasets filtered with different parameters (Fig. 6). Using DArT 409 SNP datasets, the prediction accuracy increased continuously as the TPS increased across all 410 three traits evaluated under a natural drought stress environment (Fig. 6). The prediction 411 accuracy slightly increased across all three traits when the TPS increased from 50% to 90%. 412 The smallest standard error was observed in prediction accuracy when 50% to 60% of the 413 training population size was used to predict the target traits of ER. For the target traits of SPH 414 and GY, 40% of the training population size gained the smallest standard error, which 415 indicated that 40 to 60% of the total genotypes assigned as the training set could achieve good 416 prediction accuracy in the DArT SNP markers. For the GBS SNP markers, a similar trend was 417 observed for GY. With the increase of the TPS, the average and median of prediction accuracy 418 showed slightly different for the target traits of ER and SPH (Fig. 6). The prediction 419 accuracies estimated in a single year are similar trend to using BLUEs across two years (Fig.  420 S1, S2). The H2 of the three traits was the highest in 2019, and the highest prediction 421 accuracy was obtained, but the improvement was very limited. SilicoDArT markers also 422 performed this analysis and got a consistent trend with the DArT SNP markers (Fig. S3). 423 According to the above results, 50% of the training population size was selected for further 424

analysis. 425
Genomic prediction accuracies estimated using 2-fold cross-validation (50% TP) for all 426 three traits are shown in Fig.7. The prediction accuracy between genotypic platforms has 427 reached a significant level for three interesting traits. In contrast, the two types of the marker 428 within the DArT-seq significantly differ in the SPH trait, which may mean the marker's 429 alignment to the reference genome results in a loss to the prediction for some traits. For three 430 traits, the prediction accuracy estimated from the DArT-seq markers was significantly higher 431 than markers from the GBS (P ＜0.01). Among the three traits, the prediction accuracy of 432 SPH was the lowest in both the DArT and GBS markers, which consisted of the lowest 433 heritability of SPH. The average prediction accuracies using the SNPs from the DArT-seq 434 were 0.27, 0.19, and 0.33 for ER, SPH, and GY traits, respectively. The result of SilicoDArT 435 is close to the SNPs from DArT-seq, which were 0.26, 0.22, and 0.33. While, the prediction 436 accuracies estimated using the GBS markers were -0.02, -0.06, and 0.20 for the traits of ER, 437 SPH, and GY, respectively. 438

Different prediction accuracy of DArT/GBS marker quality 439
The result of prediction accuracy estimated in all the traits under three marker sets filtered 440 with combinations of MAF and missing rate was presented in Table 3, Table S2 and Table S3. 441 For the ER trait, the average prediction accuracies range from 0.19 to 0.29, and the highest 442 accuracy corresponds to 2,762 markers in SilicoDArT. Interestingly, higher strict of missing 443 values increases the prediction accuracy, while MAF had little influence on the prediction 444 accuracy. This trend has also been found in ER and SPH traits using three markers; for GY 445 traits, the trend weakens. In most cases, the prediction accuracy is not significantly improved 446 when using more stringent filter markers except SPH. SPH is the only trait that responds to 447 high marker quality among the ER, SPH and GY traits using three marker sets. 448 The LD decay distances of DArT SNP and GBS SNP are 6.16Kb and 3.83Kb, 449 respectively (Fig.S4). The result of prediction accuracy estimated in all the traits under the 450 different DArT/GBS markers filtered by LD is presented in Table 4. SilicoDArT did not 451 perform this filter because the location information on each chromosome was missing. 452 Different LD filtering criteria correspond to different marker densities, and with the strictness 453 of filtering, the number of markers dropped dramatically. SPH obtained higher average 454 prediction accuracy with the original number of markers of about 1/10 using SNPs from 455 DArT-seq. For GBS markers, the traits of ER and GY were filtered by LD0.1, and the 456 prediction accuracy of 2,333 SNPs was higher than that of 91,003 SNPs before strictly 457 filtering. A small amount of strict markers improved the prediction accuracy on both 458 platforms, but the traits were not the same. The filtering criteria of marker quality may vary 459 due to the differences of genotyping platforms, which need to be further explored. 460

Prediction accuracy estimated within subgroups 461
In general, the population with close relationships is more likely to obtain high prediction 462 accuracy. We want to control the population structure and improve the prediction. Due to the 463 NCCP not having too many lines, it should be divided into two subgroups, accounting for 464 both population structure and population size. The PCA plots show the distribution of two 465 subgroups using PC1 and PC2 in the DArT and GBS markers, and the first two PC explained 466 35% and 31% of the genetic variation of DArT and GBS markers, respectively (Fig.8) The prediction accuracies for all the traits estimated within subgroups are shown in Table  470 5. Using the DArT SNP markers, the prediction accuracies in Subgroup 1 with 50% of TPS 471 were 0.17, 0.11, and 0.11 for ER, SPH, and GY traits, respectively. For Subgroup 2, the 472 prediction accuracies were 0.29, 0.20, and 0.33 for ER, SPH and GY traits, respectively. The 473 prediction accuracy of Subgroup 1 was significantly lower than when using the whole 474 population for three traits, whereas Subgroup2's prediction accuracies were higher than those 475 using the entire population. The result of the GBS marker is different from that of DArT, and 476 the variation of prediction accuracy is not apparent in the two subgroups. heritabilities of the three traits were more than 0.8 in 2019, but the average prediction 504 accuracies were only slightly higher than that in 2020 and across two years. This means that 505 the prediction of drought resistance needs to accumulate more year data to improve the 506 heritability across years and stable phenotype. 507 In the same trait, the prediction accuracy estimated from the DArT markers was higher 508 than that estimated from the GBS marker dataset, and the difference was significantly. Both 509 DArT and GBS are cost efficient and high-throughput genotyping platforms that can be used genotyping platforms, about $30 per sample, to control the budget within an acceptable range 515 for breeders. This price can obtain DArT markers with normal quality, but a little bit lower 516 quality for GBS markers. That means the genome coverage from GBS has not reached 1×, 517 which is probably the reason why the prediction accuracies are not suitable when using GBS 518 markers in this research. Our research indicates that DArT genotyping platform is more 519 suitable than GBS when working with a lower budget. GS using DArT markers is also being 520 implemented in CIMMYT maize breeding programs for improving grain yield and the other 521 major agronomic traits. 522 This study showed that the average prediction accuracy increased as the TPS increased in 523 all the traits in both genotyping platforms. In contrast, the standard error was first increased 524 and then decreased (Fig. 6). Relatively high prediction accuracies with the slightest standard 525 error were observed in all traits when 50% to 60% of the total genotypes were used as a 526 training set. Results of this study are consistent with previous reports (Guo et al. 2020) that 527 can be applied in the real breeding program. The marker density is an essential effect on 528 prediction accuracy if population sets and heritabilities are the same. Our research found that 529 the prediction accuracy of about 8,000 markers from DArT is significantly higher than that of 530 90,000 markers from GBS, which emphasized the importance of marking quality and 531 indicated 8,000 markers is enough to implicate genomic prediction in maize. Low-density 532 marker populations appeal to GS because they decrease multicollinearity and computational 533 time consumption and allow more individuals to be genotyped for the same cost. However, 534 this prediction was made in a breeding program's initial stage, and higher marker density 535 could play a role in the following breeding cycle for the long term. As a mature genotyping 536 platform, GBS markers have a superb performance on genomic prediction (Liu et al. 2021). 537 The higher marker density of GBS also has a high application prospect in the multi-cycle 538 breeding process. 539 We tried to increase the prediction accuracy by improving the marker quality and controlling 540 the population structure, but at the same time, the marker density was also reduced. Our 541 results show that the prediction accuracy improved in some cases (Table 3, 4, 5), even if the 542 density of the marker was reduced. However, this strategy doesn't always work, and may be 543 due to the fact that some complex traits require markers with specific effects. Also, there is a 544 tradeoff between the number of markers and marker quality or population structure because 545 marker quality decreases as the number of markers increase in a specific marker dataset. This 546 study suggests that genomic prediction can also be performed using high-quality and 547 low-density markers in the initial breeding stage under the drought environment, especially 548 when heritability is low. 549 The low-cost DArT platform increases the possibility of integrating genomic selection 550 strategies into practical breeding programs for small breeding companies and the private 551 sector. The SilicoDArT markers It is possible to capture lost effects using a single reference 552 genome for a diversity-rich population, which is essential for challenging drought prediction. 553 The key to the genomic prediction of drought resistance is to improve heritability and 554 stabilize the phenotype in maize. Natural drought experimental sites can increase the 555 heritability of planting materials by accumulating data for many years without a greenhouse. 556 Fuxin county provides a chance to study genomic prediction for drought tolerance. The 557 experimental design and model for drought resistance also need to be further developed. 558

707
The datasets generated during and/or analysed during the current study are available from the 708 corresponding author on reasonable request. 709 Figure 1 Distribution of temperature and rainfall in maize planting plots of Fuxin Mongolian Autonomous County in 2019 and 2020. The horizontal axis shows the date from sowing to harvest. The panels of a and c are Temperature (°C) information, including maximum temperature, minimum temperature, and average temperature. The panels of b and d are rainfall information, and the unit is mm.