Superoxide dismutase (SOD): from genomics and expression profiling to 2 heritability and phenotypic variability in triticale under drought stress

The main objectives of the current study were to find the possible structural association between 26 the activity of enzymatic antioxidants and grain yield of triticale plants along with identify the 27 genotypic variability of SOD activity among different triticale elite genotypes using novel 28 statistical methods. Additionally, the difference between expression rate of SOD isozymes (Mn- 29 SOD, Cu/Zn-SOD, and Fe-SOD) were tested to find out about the relationship between expression 30 rate and the drought resistance of triticale selected genotypes. The final goal of this study was to 31 introduce source code in SAS language to estimate the genetic parameters based on combined 32 ANOVA model that could be applicable in all field studies. capable of distinguishing elite genotypes and pointing out the proper features for selection criteria. Bioinformatics results indicated that SOD is more important than other enzymatic antioxidant for 53 being considered as selection criteria or candidate gene for transgenic purposes. Based on 54 expressional results, Mn-SOD announced as general isozyme that probably is highly expressed in 55 most of the species, while, Cu/Zn-SOD as was introduced as a genotype specific isozyme that is 56 likely more expressed in tolerant genotypes of triticale.


24
Background 25 The main objectives of the current study were to find the possible structural association between 26 the activity of enzymatic antioxidants and grain yield of triticale plants along with identify the 27 genotypic variability of SOD activity among different triticale elite genotypes using novel 28 statistical methods. Additionally, the difference between expression rate of SOD isozymes (Mn-29 SOD, Cu/Zn-SOD, and Fe-SOD) were tested to find out about the relationship between expression 30 rate and the drought resistance of triticale selected genotypes. The final goal of this study was to 31 introduce source code in SAS language to estimate the genetic parameters based on combined 32 ANOVA model that could be applicable in all field studies. 33 Methods 34 Thirty genotypes of triticale were studied during six years in three different locations to model the 35 relationship of enzymatic antioxidants and grain yield under drought stress condition. In every 36 year, two separate field trials were used to assess the effect of normal irrigation (irrigation trial) 37 and drought stress conditions (stress trials) by withholding irrigation at early heading stage. 38 Accordingly, based on the results of genetic variability, heatmap analysis, biplot graph, and 39 clustering, two genotypes with the highest genetic distance were selected to appraise the 40 differential expression profiling of three SOD isozyme in shoot and root organs. 41

Results 42
Field experiments and bioinformatics results showed that superoxide dismutase (SOD) was the 43 most influential antioxidant in triticale tolerance to drought stress. Additionally, Mn-SOD showed 44 high expression levels for both genotypes but with higher rates in tolerant genotype. Also, Mn-45 SOD expression rate in shoot was higher than root, but the expression level of Fe-SOD was 46 conversely higher in the root. Additionally, Fe-SOD expression levels showed no significant 47 difference in two used genotypes. The expression level of Cu/Zn-SOD was higher in root of the 48 tolerant genotype while lower in root of susceptible genotype, as compared with their shoot. 49

Conclusion 50
Heatmap analysis that is applied for the first time in a screening program, showed to be highly 51 capable of distinguishing elite genotypes and pointing out the proper features for selection criteria. 52 Bioinformatics results indicated that SOD is more important than other enzymatic antioxidant for 53 being considered as selection criteria or candidate gene for transgenic purposes. Based on 54 expressional results, Mn-SOD announced as general isozyme that probably is highly expressed in 55 most of the species, while, Cu/Zn-SOD as was introduced as a genotype specific isozyme that is 56 likely more expressed in tolerant genotypes of triticale. 57 Keywords: biplot, heatmap, genetic variability, Mn-SOD, Cu/Zn-SOD, Fe-SOD, promotor 58 analysis, feature selection. 59

60
High demands for agricultural products as a result of growing population has led researchers to 61 apply and review altered procedures and approaches so as to increase the production and 62 adoptability of the crops. Many efforts have been done in order to combine varied capabilities of 63 different plant species into one unique plant. Subsequently, over 100 years ago scientists had been 64 able to made and release a new plant species by crossing between wheat (Triticum sp.) and rye 65 (Secale cereale) plants . The intention was to increase the capability of wheat, 66 as one of the most significant sources of food amongst cereals in the world, to resist the harsh 67 environmental conditions such as drought stress. Different crossing between rye and wheat with 68 altered polyploidy levels resulted in different types of plants with different properties inherited 69 from different genome sets (Stepochkin, 2019). Triticale (×Triticosecale Wittmack) is a hexaploidy 70 species (AABBRR) resulted from crossing between tetraploid wheat (T. durum Desf; 2n=28: 71 AABB), known as durum wheat, and diploid rye (2n=14: RR). Final release of triticale inherited 72 favorable features and properties from both of its progenitors such as higher resistance to drought 73 stress and higher plant production in comparison with wheat plants (Mergoum et al., 2019). 74 Resistance to drought stress is a complex process involving different types of physiological and 75 molecular networks and responses that has not been properly unveiled yet. Drought stress, similar 76 to the most of other stresses, indirectly causes oxidative stress in plants as a result of disturbing 77 the balance between production of oxygen radicals and their detoxification (Behera et al., 2020;78 Saed-Moucheshi et al., 2014b). In such condition, defense systems and stress resistance 79 mechanisms in plants are activated to increase the amount of radicals' detoxification. Accordingly, 80 the content of antioxidants in plant cells are increases as the results of activating and/or enhancing 81 the expression of related genes. Superoxide dismutase (SOD) is family of genes that its enzymatic 82 products are able to actively dismutase the superoxide radicals (O2 -) in different plants organelles 83 (Elshafei, 2020;Hossain et al., 2015). According to numerus scientific reports, SOD family 84 contains three main members i.e. magnesium-SOD (Mn-SOD), copper/zinc-SOD (Cu/Zn-SOD), 85 and iron-SOD (Fe-SOD The main objectives of the current study were to find the structural association between the activity 95 of enzymatic antioxidants and grain yield of triticale plants along with identify the genotypic 96 variability of SOD activity among different triticale elite genotypes under normal irrigation and 97 drought stress conditions by means of novel statistical methods. Also, the promotor analyzing of 98 enzymatic antioxidants: SOD, glutathione reductase (GR), peroxidase (POD), catalase (CAT), and 99 ascorbic peroxidase (APX) was carried out to evaluate the variation in regulatory network among 100 these antioxidants in line with identify the key regulatory elements in these pathways. In addition, 101 the difference between expression rate of SOD isozymes (Mn-SOD, Cu/Zn-SOD, and Fe-SOD) in 102 selected triticale genotypes, with highest genetic distance, were tested to find out about the 103 relationship between their expression rate and the resistance of triticale plants to drought stress. 104 The In every year, 28 advanced elite genotypes, originated from international CIMMYT center, and 132 two commercial cultivars (Sanabad, and juanilo) of triticale ( were then imported into excel files for further investigation and analysis. and MDH were reported as micro mol per gram (µmol/g) fresh weight (FW), the amount of 207 chlorophyll and protein as milligram per gram (mg/g) FW, and proline content as micromolar per 208 gram (µM/g) FW. The activity of enzymatic antioxidants was reported as micromole of free 209 radicals decomposed in mg of total protein sample (u/mg Prt). In addition, at the harvest time, 210 grain yield of the two central rows in each plot were harvested and their grain yield were recorded 211 according to kilogram per square meter (kg/m 2 ). 212 After sampling the youngest leaves of the plants in greenhouse study at each time span, after 213 induction of the drought stress, the activity of SOD enzyme was measured along with the relative 214 expression of SOD isozymes. For extracting the RNA from leaf samples, Dena Zist isolating kits 215 were used. After the extraction, the quantity of extracted RNA was measured by nanodrop method. 216 Following the quantity measurement of RNA, the extracted samples were treated with DNase in 217 order to remove any DNA contamination, afterward, the positive and negative quality control (QC) 218 of the RNA in each sample were performed; the positive QC was checked by direct electrophoresis 219 of 3 mL extracted samples on gel (1% of 0.5 TBE buffer) and observing two separate and sharp 220 bands on the gel representing each of 16s and 24s RNA in the samples, while, the negative QC 221 was checked by observing no band on electrophoresis gels posterior to polymerase chain reaction 222 (PCR) of samples conducted by using the primer of internal control (housekeeping: elongation 223 factor 1) gene that was used in this study. Then, the RNA of the extracted samples was transformed 224 into complementary DNA (cDNA) by use of Fermentas Co. kits. The produced cDNA's were then 225 used for measuring the relative expression of three isozymes of SOD gene i.e. Mn-SOD, Cu/Zn-226 SOD, and Fe-SOD based on RT-PCR method. The corresponding pair of primers for each isozyme 227 was designed in AlleleID software according to the sequence of the isozymes in wheat genome. In 228 this process, the elongation factor 1 (EF1) was used as the internal control (housekeeping) gene. 229 Each pair of forward (F) and backward (revers: R) primers that were used for expression analysis 230 of The obtained data related to field experiments were subjected to combined analysis of variance 239 (combined ANOVA) and descriptive statistic's, such as mean and standard error. After that, the 240 mean squares (MS) obtained from combined ANOVA for each of normal irrigation and drought 241 stressed trials were used for estimating the residual (res) variance, variance of block (blk) within 242 year (yr), genotype (gn) by year variance, genotypic variance, environmental variance, phenotypic 243 variance, broad sense (general) heritability, genotypic coefficient of variation (GCV), and 244 phenotypic coefficient of variation (PCV) based on the following formulas (Fan,  Factorial ANOVA based on completely randomized design was used for analyzing expression data 254 by using proc GLM in SAS software. The least significant difference (LSD) method was then used 255 for mean comparison of expression data. Excel software package was used for drawing graphs 256 regarding the mean comparisons alongside evaluating bioinformatics data related to promotor 257 analysis. The mean values of triticale genotype across all years for all measured features were used 258 for building a regression model based on stepwise feature selection method in each of normal 259 irrigation and drought stress conditions. Th stepwise regression model was estimated in SAS 260 software by using proc REG. Heatmap by use of "gplots" and "agricolae" (for standardizing data) 261 libraries, biplot graph based on principal component analysis (PCA) by use of "factoextra" library, 262 and cluster analysis based on Ward method and Euclidian distance by use of "Nbclust" library 263 were performed in R statistical software (R 3.5). 264

265
3.1.ANOVA and genetic variability 266 The first three out of six years of the field study were performed in one location (Shiraz), the fourth 267 year in another location (Zarghan), and the last two years in a separate location (Sanandaj); 268 therefore, the effect of location, as a source of variance, in combined ANOVA is actually nested 269 inside the effect of year that makes it impossible to separate these two effects from each other. 270 Subsequently, the effect of year, as a representative of both year and location, was considered 271 along with the effect of condition (stress), normal irrigation or drought stress condition within each 272 year, and the effect of genotype in final ANOVA model as the main source of variation (Table 3). 273 The effect of year was statistically significant (p<0.05) for all measured features in triticale plants 274 consisting of H2O2, MDH, chlorophyll, carotenoid, protein, proline, POD, CAT, APX, GR, SOD, 275 and grain yield which indicates the overall difference between years and locations regarding these 276 features. Similarly, there was significant differences between normal irrigation and drought stress 277 conditions regarding all features in Table 3. However, no significant effect of stress × year 278 interaction was found for chlorophyll content, protein content, POD and APX activities. Save for 279 proline content, the main effect of genotype along with the interaction effects of year × genotype, 280 stress × genotype, and stress × year × genotype were statistically significant for all other measured 281 features (Table 3). 282 Genetic parameters were estimated by taking year and genotype as random effects and calculating 283 the expected values for each source of variation in ANOVA mixed model separately in each stress 284 condition (Table 4). Accordingly, environmental, genotypic, and phenotypic variance were 285 extracted from data and transformed into standard indices i.e. heritability, phenotypic coefficient 286 of variation (PCV), and genotypic coefficient of variation (GCV). The highest heritability under 287 normal irrigation and drought stress conditions were observed in carotenoid content (62%) and 288 SOD activity (61%), respectively. The protein content showed the lowest percentage of heritability 289 in both normal irrigation (6.81%) and drought stress conditions (4.15%). Under normal irrigation 290 condition, proline content had the highest PCV (70%) and GR activity showed the highest GCV 291 (26.6%). The lowest PCV (7.5%) and GCV (1.6%), on the other side, were obtained for H2O2 in 292 normal irrigation condition. MDH and CAT activity were the feature with the same lowest 293 percentage of PCV (33.6%) in drought stress condition. The highest PCV percentage under 294 drought stress condition was observed in APX activity. SOD activity showed the highest GCV 295 (46.4%) in drought stress condition, while the lowest in this condition was estimated for protein 296 content (2.19%). 297 The average values for all triticale genotypes regarding SOD activity, as the most important 298 variable of the study as a result of the highest heritability and genetic variability in drought stress 299 condition, and grain yield across all six years experiments doubled with their standard errors for 300 each stress condition are presented in Table 5. Genotypes 21 (849.5 g/m 2 ) and 17 (625.1 g/m 2 ) 301 under normal irrigation condition and genotypes 19 (380.6 g/m 2 ) and 7 (257 g/m 2 ) under drought 302 stress condition showed the highest and the lowest grain yield, respectively. Genotype 3 showed 303 above average grain yield for normal irrigation condition while lower than average point for stress 304 condition; whereas, the grain yield of genotype 28 was above average point under both normal and 305 stress conditions. Genotype 28 showed the highest SOD activity than all other triticale genotypes 306 in both normal irrigation (17.5 u/mg protein) drought stress conditions (49.9 u/mg protein). The 307 lowest SOD activity for normal and stress conditions were obtained in genotype 1 (4.4 u/mg 308 protein) and genotype 3 (16.4 u/mg protein), respectively (Table 5). 309

3.2.Feature selection and promotor analysis 310
In order to investigate the effects of biochemical features on grain yield of triticale genotypes and 311 finding the most influential features, stepwise regression analysis was performed for normal 312 irrigation and drought stress conditions (Table 6). Accordingly, the results showed that the total 313 protein content was the sole feature selected as the most influential feature on grain yield in normal 314 irrigation condition. The model r-square, known as coefficient of determination, reached to 89% 315 showing a high impact of total protein content on grain yield. Also, the regression coefficient for 316 protein content was negative (b = -27.79) showing a decreasing effect for this feature toward grain 317 yield under normal irrigation condition. On the other side, stepwise regression model for drought 318 stress condition selected the content of free proline, SOD and GR activities as the influential 319 feature on grain yield of triticale in this condition (Table 6). At the first step of the feature selection 320 for drought stress condition, SOD activity was inserted into model showing that this feature is the 321 most influential independent variable toward the grain yield of triticale. After that, proline and GR 322 activity were respectively inserted into the stepwise model as the second and the third influential However, the results achieved by the appropriate statistical methods are clearly indicated by the 360 phenotypic properties and therefore, molecular and genotypic investigations would enable us to 361 test these results. Consequently, promotor analysis of the used enzymatic antioxidants in this study 362 was appraised so as to find any different or similar pattern regarding their regulation mechanism 363 to explain and/or test the actual influence of SOD in drought stress condition. Unfortunately, there 364 is still no bioinformatic database for triticale genome sequence, so we performed the bioinformatic 365 investigation of promotor analysis in wheat genome, its parental line that shares the same AA and 366 BB genomes with triticale. The regulatory elements that were identified in the promotor regions 367 of the enzymatic antioxidants along with their corresponding carrier genome sets are presented in 368 Table 7. The results showed that the greatest number of the enzymatic antioxidant genes were 369 distributed among A and D genome groups. Also, the results revealed a significant difference 370 between the regulatory elements of SOD genes and the regulatory elements of the other antioxidant 371 genes. A high number of regulatory elements in promotor regions of SOD genes were related to 372 circadian regulation, abscisic acid production, auxin response pathways, MYB binding sites, and 373 methyl jasmonate responsive genes that are important compounds of stress responses in plants and 374 were mostly distributed in A chromosome set. In addition, gibberellin and light responding related 375 regulatory elements showed high distributions in SOD promotor regions, especially in A genome 376 set. The total numbers of regulatory elements of SOD promotor regions (1441) were significantly 377 higher than any other considered antioxidant (Table 7). 378

3.3.Genotype selection for expression analysis 379
According to the heritability, genotypic variability, using advanced biostatistics methods and 380 different multivariate models alongside bioinformatics approaches, it was clarified that SOD was 381 the most effective antioxidant under drought stress condition in triticale. Therefore, for achieving 382 higher triticale grain yield under drought stress condition, SOD activity has the capability of being 383 used as an indirect criterion, in order to be used in breeding programs of triticale. Therefore, for 384 differential expression analysis goals regarding the genes involved in the SOD response pathways, 385 all of the used triticale genotypes in this study were assessed to find those with the highest genetic 386 distance. Accordingly, the similarities and distances between all pairs of 30 triticale genotypes 387 were considered based on the Euclidean distance method and then, they were transformed into two 388 separate cluster dendrograms, representing both normal and stress conditions (Figure 3). Under 389 normal irrigation condition, three distinct clusters were distinguished (Figure 3a). This cluster 390 dendrogram indicated that genotypes 28 and 1 were the two cases presenting the highest genetic 391 dissimilarity, or distance, based on all measured featured. However, genotypes 2 and 3, that were 392 grouped in the same clustered as the genotype 1, were two alternatives for being considered as the 393 genotypes revealing the highest genetic distance from genotype 28, rather than genotype 1. Under 394 drought stress, the highest genetic distance was observed between genotype 28 and genotype 3 395 (Figure 3b). Similar to the normal condition result, genotypes 1 and 2 were placed in the same 396 cluster with genotype 3 under drought stress condition. Accordingly, the alternatives for genotype 397 28 could be genotype 17 and genotype 19 that were in same cluster as the genotype 28. As a result, 398 genotype 1, 2, and 3 versus genotypes 17, 19, and 28 were the genotypes with the highest genetic 399 distances that could be considered in the expressional study. These results were clearly verified by 400 the results obtained from the heatmap analyses ( Figure 1) in which genotype 1, 2, and 3 showed 401 the similar colors for almost all measured features, while they showed very dissimilar colors from 402 genotypes 28, 17, and 19 ( Figure 1a and figure 1b). Furthermore, considering the colors saturation 403 in heatmap for these two genotype sets, revealed that they showed the highest dissimilarity 404 regarding the activity of SOD enzyme. However, we were obliged to select just two nominees as 405 the genotypes with the highest dissimilarity for differential expression study. This requirement 406 drew out our attention into biplot and mean comparison for SOD and grain yield. The biplot graphs 407 indicated that, under both normal irrigation and drought stress conditions, genotype 28 was the 408 closest to SOD alongside the grain yield, among the mentioned genotypes ( Figure 2a). Conversely, 409 genotypes 1, 2, and 3 showed similar distance from SOD and grain yield under both normal and 410 stress conditions. As the last option for finding the most proper genotypes with highest genetic 411 variability, mean comparison for SOD and yield were taken into consideration (Table 4). Finally, 412 according to the results of mean comparison for triticale genotypes under both stress conditions, 413 genotype 28 and genotype 3 were chosen to be applied for investigation of SOD isozyme 414 differential expression. 415

3.4.Differential expression analysis of SOD isozymes 416
The relative expression analysis based on real-time PCR method for SOD isozyme containing Mn-417 SOD, Cu/Zn-SOD, and Fe-SOD was carried on two selected genotypes, genotype 28 and 3, with 418 highest genetic variability with each other. The analysis of variance for the SOD isozymes were 419 carried out based on a factorial experiment, in which genotype and time spans after inducing 420 drought stress were the two experimental factors. The main effect of genotype on Mn-SOD and 421 Cu/Zn-SOD gene expression was significant for both shoot and root tissues, but it did not show a 422 significant effect on Fe-SOD gene expression of these tissues (Table 8). Mn-SOD and Cu/Zn-SOD 423 in root tissue were not significantly influenced by the main effect of time after the drought stress, 424 whereas, these two isozymes in the shoot tissue coupled with the Fe-SOD isozyme in both root 425 and shoot tissues showed significant responses to the effect of time spans after drought stress. 426 Except Fe-SOD expression rate in both shoot and root tissues, the two other SOD isozymes showed 427 no significant response to the interaction effect of genotype × time span (Table 8). 428 The Mn-SOD relative expression in both triticale genotypes showed an increasing trend in 429 response to longer duration of time after drought stress in both shoot and root tissues (Figure 4a). 430 The expression rate of Mn-SOD was higher in shoot than expression rate in the root tissue. 431 Additionally, genotype 28 showed significantly higher expression rate regarding Mn-SOD 432 isozyme in both root and shoot in comparison with genotype 3. The expression rate of Cu/Zn-SOD 433 gene showed different patterns regarding two used triticale genotypes in both shoot and root tissues 434 (Figure 4b). In genotype 3, Cu/Zn-SOD expression rate in shoot increased by changing the time 435 span from 12 h to 36 h, but there was no significant different between 36 and 72 h in this tissue. 436 Quite the contrary, 36 h showed the lowest Cu/Zn-SOD expression rate than the two other time 437 spans in root tissue of genotype 3. Quite the contrary, the expression rate of Cu/Zn-SOD isozyme 438 in both root and shoot tissues showed a growing (inducing) pattern in response to the longer time 439 span after the drought stress in genotype 28. In both shoot and root tissues of the both genotypes, 440 Fe-SOD expression rate was reduced in response to longer time span, but the decreasing rate in 441 was slower in genotype 28 than the other genotype (Figure 4c). Comparison the expression rates 442 of the three assessed isozymes of SOD gene with regard to the genotypes, showed significantly 443 higher expression rate of Mn-SOD and Cu/Zn-SOD in genotype 28 than that of genotype 3. The 444 relative expression of Fe-SOD, on the other hand, showed no significant difference between both 445 genotypes. 446 It should be noted that the measurement of SOD activity in shoot and root tissues showed an 447 increasing pattern in response to longer time span after stress induction in both genotypes. In other 448 words, the longer the drought stress lasted, the higher the activity of SOD enzyme reached. 449 Nevertheless, the SOD activity was significantly higher in genotype 28 in comparison with the 450 genotype 3. Moreover, shoot tissue showed a higher activity of SOD enzyme than that of the root 451 tissue (the data are not presented). 452

453
In this study, grain yield of triticale and all biochemical features significantly responded to the 454 effects of environment and genotypes as the main sources of variations in the combined ANOVA. 455 This result is indicating the differences among triticale genotypes in response to altered 456 environments, which is alluding to the proficiency of the screening programs for introducing high 457 potential genotypes. Thus, the estimation of genetic parameters would alter in different conditions. 458 Consequently, heritability, PCV, and GCV of almost all features showed significantly higher 459 values in drought stress condition than that in the normal irrigation condition which indicates that 460 screening and breeding programs, in this set of triticale genotypes, would be more efficient in 461 drought stress condition. On the other hand, the main aim of the current study was to find proper 462 methods to reach out to high potential genotype to tolerate drought stress condition. Accordingly, 463 the assessments were done for both normal and stress condition, but the results of the drought 464 stress condition superseded over the results of the normal irrigation condition, if there was any 465 conflict between the results of these two conditions, such as selecting the influential features on 466 grain yield or finding the highest heritable feature. Estimation of genetic parameters for both 467 conditions showed that SOD activity had the highest heritability and GCV under drought stress 468 condition, while it had almost low values of these parameters under normal condition. These 469 results indicate the importance of this enzyme in triticale genotypes under water shortage 470 condition, because GCV shows high variabilities among genotypes and the higher screening ability 471 for finding appropriate genotypes in this regard, also, its high heritability shows high efficiency of 472 screening based on this feature and higher response to selection (RS), than any other feature, in 473 the following generations. In agreement with the result of genetic variability, heatmap analysis 474 clearly showed higher variations of colors among triticale genotypes in relation to SOD activity 475 under drought stress condition. Furthermore, grain yield showed significantly low heritability and 476 genetic variability in comparison with most of other biochemical features, especially SOD activity. 477 Consequently, direct screening for grain yield would not confidently result in higher grain yield in 478 the next generation of breeding program, but using an indirect selection criterion such as SOD for 479 yield would be highly efficient, on condition that it shows a high collinearity and genetical 480 association, either positive or negative, with grain yield. In concordance with heatmap results that 481 showed similar color patterns and saturations for grain yield and SOD activity, and also the result 482 of biplot in which SOD was the closest feature to the grain yield in drought stress condition, the 483 results of the feature selection pointed out that SOD was the most influential feature, with positive 484 impact, on grain yield of triticale. As a result, the capability of SOD for being used as an important 485 indirect criterion for grain yield is verified by all used biometrical and biostatistical methods. To 486 put it simply, genotypes that show higher SOD activity under drought stress condition probably 487 have higher grain yield. 488 Thus far, there are many studies concerning with the effect of environmental stresses on triticale 489 plants, however, few studies have focused on the impacts of biochemical features on grain yield 490 and their capability of being used as indirect criterion; and there is still no study, to the best of our 491 knowledge, considering the heritability and genetic variability of biochemical features in triticale 492 genotypes. In addition, there is still no related program or source code for estimation of genetic 493 parameters e.g. heritability and genetic variability, and therefore, the SAS programming code that 494 was used in the current study (supplementary materials) is the first program ever written and 495 practically used to estimate genetic parameters. The results achieved by this code was successfully 496 tested by the results of the hand calculations and the both results were in agreement, with just a 497 little difference as a result of different decimal rounding. It also should be noted that, application 498 of heatmap in the current study and using color saturation for identifying the association between 499 yield and biochemical feature in order for finding possible selection criteria, has never been applied 500 by any other study before. Consequently, our results showed the efficiency heatmap analysis in 501 screening and breeding programs and also showing a combination of relationships among features, 502 relationships among genotypes, and cross relationships between features and genotypes. In The results of genetic parameters, feature selection, heatmap, and biplot graph alluded to SOD as 517 the most important antioxidant in triticale under drought stress condition. However, these results 518 are based on phenotype and probable products of gene actions, and they cannot show the molecular 519 changes inside the cells. Therefore, in order to evaluate the antioxidants at the genomic and 520 molecular level and revealing any probable similar or dissimilar pattern regarding their regulation 521 mechanisms explain the aforementioned results, the upstream sequences of antioxidant genes were 522 appraised by bioinformatics tools. This method of appraising upstream sequences of a gene is 523 known as promotor analysis and important region for regulating the gene action. The result of 524 promotor analysis showed clear differences between promotor region of SOD and other 525 antioxidant genes. The number of regulatory elements in SOD promotor region was significantly 526 higher than the other antioxidants, mostly by more than twofold. Also, most of the regulatory 527 elements in SOD promotor region were compounds related to stress responsive pathways. The 528 numerous regulatory elements in promotor region of SOD indicates that this gen can be regulated 529 by various responsive pathways and is most likely more active than other antioxidant in response 530 to changing conditions. In addition, high number of stress related elements in this region shows 531 that this gene is highly responsive to stress condition and it may be highly increased by occurring 532 stress condition. The result also showed a great number of MYB binding sites in the promotor 533 region of SOD genes, that was significantly higher than other antioxidant genes. MYB elements 534 are specifically acting in the responses to water shortage condition, verifying the importance of 535 SOD under drought stress condition. Accordingly, to regulate the expression of SOD genes, a number of cis acting regulatory elements 557 must be present in the promoter region of this gene to regulate its expression based on intracellular 558 changes. Therefore, higher number of regulatory elements that are able to bind with the regulatory 559 motifs in the promotor region bring about higher sensitivity of the corresponding gene to 560 intracellular changes, such as what was identified for SOD promotor in this study. A recent study 561 (Hu et al., 2019) on the expression pattern of SOD genes in grape wine, showed that its regulatory 562 mechanisms could be active at the three levels of transcription, post-transcription, and translation, 563 and their functional roles are focusing on response to stresses during plant development indicating, 564 in concordance with our results, the high importance of SOD expression under environmental 565 stress conditions and the effect of the regulatory elements in this pathway. 566

567
The overall results of biometrical, biostatistical, and bioinformatics analyses indicated that SOD 568 was the most important antioxidant in response to drought stress condition. Also, combining the 569 final results regarding clustering, biplot graph, and heatmap analysis along with the mean 570 comparison for SOD activity and triticale grain yield led to the selection of two genotypes, i.e. 571 genotype 28 and genotype 3, with the highest genetic distance in order to assess the deferential 572 expression profiling of SOD genes in triticale. Considering the activity of SOD enzymes, in the 573 two selected genotypes, showed that the activity of SOD was increased in both genotypes in 574 response to longer time span from initial point of drought stress, but the rate of increase in genotype 575 28 was over threefold higher than genotype 3. Accordingly, the rate of relative expression of Mn-576 SOD was about threefold, and Cu/Zn-SOD was about twofold higher in genotype 28 than those in 577 genotype 3. The expression levels of Fe-SOD showed no significant difference between the two 578 genotypes. These results indicated that higher activity of SOD enzyme under drought stress 579 condition is most likely due to the higher expression level of Mn-SOD and/or the Cu/Zn-SOD 580 isozymes, while Fe-SOD is, seemingly, not significantly involved in response to the drought stress. 581 In addition, quite the contrary to Mn-SOD and Cu/Zn-SOD, the expression pattern of Fe-SOD was 582 decreased in response to longer time span from drought stress starting point. This result indicates 583 that Fe-SOD is probably more responsive to non-stress signaling pathways and its expression level 584 under drought stress is decreased so as to let the two other SOD isozymes to take action and grow 585 their expression level more rapidly. Furthermore, a comparison between Mn-SOD and Cu/Zn-586 SOD indicated that the expression level of Mn-SOD is continuously increased as the time span 587 increased in both genotypes; while, the Cu/Zn-SOD showed altered pattern in response to different 588 time spans in each genotype. Therefore, these results are possibly pointing out that Mn-SOD is a 589 general responsive isozyme that its expression level is, somewhat, quite independent of the genetic 590 differences between the genotypes within a species, on the other hand, Cu/Zn-SOD is more 591 dependent on genotypic differences between the two genotypes. This result regarding the 592 expression levels of Cu/Zn-SOD made us to announce it as a genotype specific isozyme, meaning 593 that its expression level is highly correlated with the genetic contents and background of its carrier 594 individual and is likely more expressed in the genotypes tolerant to drought stress. Furthermore, 595 the expression level of Mn-SOD in shoot and root tissues was roughly equal. Contrarily, the 596 expression rate of Cu/Zn-SOD in genotype 3 was higher in shoot than that in root, while its 597 expression level in genotype 28 was significantly lower in shoot than that in the root tissue. The 598 first organelles that are affected by drought stress are normally those in the shoot (Du et al., 2020), 599 consequently, this leads the expression levels of Mn-SOD to be significantly higher in shoot. 600 Higher duration of drought stress causes the stress to reach out more to the root system and involve 601 its organelles as well, as a result, the expression level of Cu/Zn-SOD isozyme in the root tissues, 602 mostly in tolerant genotypes, is increased. Conversely, the expression levels of Fe-SOD were 603 continuously higher in root in comparison with the shoot tissue. This result increases the likelihood 604 involvements of Fe-SOD isozyme in non-stressful responsive pathways, because, as mentioned 605 before, the root organelles are less likely affected by negative impacts of drought stress than the 606 shoot organelles. 607 SOD has been reported to be a key enzyme in regulating intracellular ROS levels and maintaining 608 normal physiological conditions under oxidative stress caused by effect of water shortage 609 condition (Mao et al., 2018). The findings of a study performed by Zhao et al. (2018) showed that 610 the SOD expression level and malondialdehyde content in rice were significantly correlated, in 611 addition, the activity of SOD expression level in rice cultivars sensitive to heat stress was much 612 lower than the resistant cultivars. Also, in agreement with our results, Mohammadi et al. (2020) 613 showed that Cu/Zn-SOD was significantly expressed in common bean plants under drought stress 614 condition, but the rate of its induction in tolerant genotype was significantly higher than the 615 susceptible one. In concordance with our result regarding the importance of Cu/Zn-SOD as a more 616 specific isozyme that its expression levels are most likely changed by genetic contents of its 617 carriers individuals, Xie et al. (2019b) stated that it is highly feasible that Cu/Zn-SOD play an 618 important role in regulating total SOD activity and ROS detoxification in stressful condition. 619 Moreover, attempts to regulate Cu/Zn-SOD levels in plant chloroplasts through the gene transfer 620 mechanism resulted in transgenic plants with altered expression associated with the aimed gene. 621 Overexpression of Cu/Zn-SOD genes has been reported to increase stress tolerance in tobacco 622 transgenic plants (Tepperman and Dunsmuir, 1990 (2018), suggested that Mn-SOD is an important agent that plays an important role in response to 630 drought tolerance, however, our results showed that Mn-SOD is possibly a general responsive 631 isozyme that its response to drought stress is similar (or the same) to other stress conditions. 632 633

634
The overall results indicated the high possibility of screening favorable and potential genotypes to 635 announce as new triticale cultivars as the result of great variabilities among use genotypes in all 636 six years of the study. Additionally, a new programming code written in SAS language was 637 successfully applied to estimate the genetic parameters based on the results of combined ANOVA 638 for the first time since their introduction. The genetic parameters showed a high heritability and 639 genetic variability for SOD activity under drought stress condition. The results of feature selection, 640 heatmap analysis, and biplot clearly indicated that SOD had a significant effect on triticale 641 production in drought stress condition. In this regard, heatmap graph, which is based on the 642 saturation of different colors that shows the association within variables and cases (genotypes) and 643 also their cross relationships, was successfully used for the first time in this study and it showed 644 to be highly appropriate for screening favorable genotypes and processing the features that are 645 proper for indirect selection in breeding programs. According to the results of promotor analysis 646 for antioxidant enzymes, there was a significant difference between SOD and other antioxidants 647 regarding the numbers and the types of motifs and regulatory elements available in their promotor 648 regions. Bioinformatics study showed that the most of regulatory elements in SOD promotor 649 regions were related to stress responsive pathways, especially drought stress as the MYB binding 650 site. Accordingly, two genotypes with the highest genetic distance were selected to consider the 651 differential profiling of three different isozymes of SOD genes in shoot and root organs of the 652 triticale genotypes. The result showed that tolerant genotype (genotype 28) had significantly higher 653 levels of expression for Mn-SOD and Cu/Zn-SOD than susceptible genotype (genotype 3), while 654 the expression level of Fe-SOD was not significant between these two genotypes. Mn-SOD 655 showed higher expression level in shoot than root, but the expression level of Fe-SOD was 656 conversely higher in the root. Whereas, the expression level of Cu/Zn-SOD in tolerant genotype 657 was higher in root while its level in sensitive genotype was lower in root, in comparison with shoot. 658 Finally, it could be stated that Mn-SOD is a general isozyme that responses positively to stress 659 condition in probably most of the plant species with different genetic backgrounds, while, Cu/Zn-660 SOD is a genotype specific isozyme that is likely more expressed in some genetic content such as 661 the tolerant genotype of triticale in the current study. Furthermore, genotype 28 can be introduced 662 as a favorable genotype for being selected for further assessments regarding its potentials in 663 drought stress condition. 664 665 writing of the final manuscript was carried out by A. Saed-Moucheshi and edited by A. A. 699 Mozafari. All authors read the manuscript and confirmed the final version. 700 7.7.