Superoxide dismutase (SOD) as a selection criterion for triticale grain yield under drought stress: a comprehensive study on genomics and expression profiling, bioinformatics, heritability, and phenotypic variability

Background The main objectives of this study were to find the possible structural association between the activity of enzymatic antioxidants and the grain yield of triticale plants as well as identifying the genotypic variability which might be effective on this association. Accordingly, expression levels of superoxide dismutase (SOD) isozymes (Mn-SOD, Cu/Zn-SOD, and Fe-SOD) were appraised to distinguish any possible relationship between SOD expression and drought resistance of triticale. A novel analytical method for distinguishing elite genotypes based on measured features was proposed. Additionally, a new programing based on SAS-language (IML) was introduced to estimate the genetic parameters rooted from combined ANOVA model (linear mixed model), which is capable of being used in any field study other than the current one. Methods Thirty genotypes of triticale were studied under normal and drought stress conditions during 6 years (three different locations). Accordingly, based on the results of genetic variability, heatmap analysis, biplot graph, and clustering technique, two genotypes with the highest genetic distance were selected to appraise the differential expression profiling of three SOD isozyme in shoot and root organs. Results Field experiments and bioinformatics results showed that superoxide dismutase (SOD) was the most influential antioxidant in resistance of triticale to drought stress; therefore, it could be used as an indirect selection index in early stages to distinguish resistant genotypes to drought stress. Additionally, Mn-SOD and Fe-SOD showed roughly similar expression levels for both genotypes under drought stress. However, Cu/Zn-SOD expression level was higher in root and shoot of the tolerant genotype than the susceptible genotype. Conclusion Heatmap analysis that is applied for the first time to screen suitable genotypes, showed to be highly 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 being considered as selection criteria or candidate gene for transgenic purposes. Based on expressional results, Mn-SOD announced as a general isozyme that is probably highly expressed in most of the species, while, Cu/Zn-SOD was introduced as a genotype specific isozyme that is likely more expressed in tolerant genotypes Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-02919-5.


(Continued from previous page)
Results: Field experiments and bioinformatics results showed that superoxide dismutase (SOD) was the most influential antioxidant in resistance of triticale to drought stress; therefore, it could be used as an indirect selection index in early stages to distinguish resistant genotypes to drought stress. Additionally, Mn-SOD and Fe-SOD showed roughly similar expression levels for both genotypes under drought stress. However, Cu/Zn-SOD expression level was higher in root and shoot of the tolerant genotype than the susceptible genotype. Conclusion: Heatmap analysis that is applied for the first time to screen suitable genotypes, showed to be highly 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 being considered as selection criteria or candidate gene for transgenic purposes. Based on expressional results, Mn-SOD announced as a general isozyme that is probably highly expressed in most of the species, while, Cu/Zn-SOD was introduced as a genotype specific isozyme that is likely more expressed in tolerant genotypes Keywords: Heatmap, Mn-SOD, Cu/Zn-SOD, Fe-SOD, Promotor analysis, Feature selection

Background
High demands for agricultural products as a result of growing population has led researchers to apply and review altered procedures and approaches so as to increase the production and adoptability of the crops [56]. Many efforts have been done in order to combine varied capabilities of different plant species into one unique plant. Subsequently, over 100 years ago scientists had been able to made and release a new plant species by crossing between wheat (Triticum sp.) and rye (Secale cereale) plants [10]. The intention was to increase the capability of wheat, as one of the most significant sources of food amongst cereals in the world, to resist the harsh environmental conditions such as drought stress. Different crossing between rye and wheat with altered polyploidy levels resulted in different types of plants with different properties inherited from different genome sets [53]. Triticale (×Triticosecale Wittmack) is a hexaploidy species (AABBRR) resulted from crossing between tetraploid wheat (T. durum Desf; 2n = 28: AABB), known as durum wheat, and diploid rye (2n = 14: RR). Final release of triticale inherited favorable features and properties from both of its progenitors such as higher resistance to drought stress and higher plant production in comparison to wheat plants [37].
Resistance to drought stress is a complex process involving different types of physiological and molecular networks and responses which has not been properly unveiled and addressed yet [40,57]. Drought stress, similar to the most of other environmental stresses, indirectly causes oxidative stress in plants as a result of disturbing the balance between production and detoxification of oxygen radicals [7,11,48]. In such condition, defense systems and stress resistance mechanisms in plants are activated to induce higher detoxification of oxygen radicals [27,54]. Accordingly, the content of antioxidants in plant cells are increases as the results of activating and/or enhancing the expression profiling of related genes [26,47,50]. Superoxide dismutase (SOD) is family of genes that its enzymatic products are able to actively dismutase the superoxide radicals (O 2 − ) in different plants' organelles [17,23]. According to numerus scientific reports, SOD family contains three main members i.e. magnesium-SOD (Mn-SOD), copper/zinc-SOD (Cu/Zn-SOD), and iron-SOD (Fe-SOD). Many studies have substantiated altered expression rates for each of the SOD isozymes under varies situation and stresses in different plants species [3,20,48] and even different genotypes within a species [12,17,54]. Therefore, delving over the SOD expressional mechanisms that are active in its responsive pathways along with identification of related regulatory elements may provide new insights toward understanding and recognizing the resistance mechanisms in plants. Additionally, finding the true association between SOD expression rate and plant productivity under drought stress condition would help breeders to identify resistant genotypes with more simplicity and increase the efficiency of breeding programs [34].
The main objectives of the current study were to find the structural association between the activity of enzymatic antioxidants and grain yield of triticale plants along with identify the genotypic variability of SOD activity among different triticale elite genotypes under normal irrigation and drought stress conditions by means of advanced and novel statistical methods. Also, the promotor analyzing of enzymatic antioxidants: SOD, glutathione reductase (GR), peroxidase (POD), catalase (CAT), and ascorbic peroxidase (APX) was carried out to evaluate the variation in regulatory network among these antioxidants in line with identify the key regulatory elements in these pathways. In addition, the difference between expression rate of SOD isozymes (Mn-SOD, Cu/Zn-SOD, and Fe-SOD) in selected triticale genotypes, with highest genetic distance and lowest similarity, were tested to find out about the relationship between their expression rate and the resistance of triticale plants to drought stress. The tolerant genotypes of triticale maintaining high production under drought stress condition, were then introduced as genotypes candidates for being used as new and favorable cultivars.

ANOVA and genetic variability
The first three out of six years of the field study were performed in one location (Shiraz), the fourth year in another location (Zarghan), and the last two years in a separate location (Sanandaj); therefore, the effect of location, as a source of variance, in combined ANOVA is actually nested inside the effect of year that makes it impossible to separate these two effects from each other. Subsequently, the effect of year, as a representative of both year and location, was considered along with the effect of condition (stress), normal irrigation or drought stress condition within each year, and the effect of genotype in final ANOVA model as the main source of variation ( Table 3). The effect of year was statistically significant (p < 0.05) for all measured features in triticale plants consisting of H 2 O 2 , MDH, chlorophyll, carotenoid, protein, proline, POD, CAT, APX, GR, SOD, and grain yield which indicates the overall difference between years and locations regarding these features. Similarly, there was significant differences between normal irrigation and drought stress conditions regarding all features in Table 3. However, no significant effect of stress × year interaction was found for chlorophyll content, protein content, POD and APX activities. Save for proline content, the main effect of genotype along with the interaction effects of year × genotype, stress × genotype, and stress × year × genotype were statistically significant for all other measured features (Supplementary Table 3).
Genetic parameters were estimated by taking year and genotype as random effects and calculating the expected values for each source of variation in ANOVA mixed model separately in each stress condition (Table 1). Accordingly, environmental, genotypic, and phenotypic variance were extracted from data and transformed into standard indices i.e. heritability, phenotypic coefficient of variation (PCV), and genotypic coefficient of variation (GCV). The highest heritability under normal irrigation and drought stress conditions were observed in carotenoid content (62%) and SOD activity (61%), respectively. The protein content showed the lowest percentage of heritability in both normal irrigation (6.81%) and drought stress conditions (4.15%). Under normal irrigation condition, proline content had the highest PCV (70%) and GR activity showed the highest GCV (26.6%). The lowest PCV (7.5%) and GCV (1.6%), on the other side, were obtained for H 2 O 2 in normal irrigation condition. MDH and CAT activity were the feature with the same lowest percentage of PCV (33.6%) in drought stress condition. The highest PCV percentage under drought stress condition was observed in APX activity. SOD The average values for all triticale genotypes regarding SOD activity, as the most important variable of the study as a result of the highest heritability and genetic variability in drought stress condition, and grain yield across all 6 years experiments doubled with their standard errors for each stress condition are presented in Table 2. Genotypes 21 (849.5 g/m 2 ) and 17 (625.1 g/m 2 ) under normal irrigation condition and genotypes 19 (380.6 g/m 2 ) and 7 (257 g/m 2 ) under drought stress condition showed the highest and the lowest grain yield, respectively. Genotype 3 showed above average grain yield for normal irrigation condition while lower than average point for stress condition; whereas, the grain yield of genotype 28 was above average point under both normal and stress conditions. Genotype 28 showed the highest SOD activity than all other triticale genotypes in both normal irrigation (17.5 u/mg protein) drought stress conditions (49.9 u/mg protein). The lowest SOD activity for normal and stress conditions were obtained in genotype 1 (4.4 u/mg protein) and genotype 3 (16.4 u/mg protein), respectively ( Table 2).

Feature selection and promotor analysis
In order to investigate the effects of biochemical features on grain yield of triticale genotypes and finding the most influential features, stepwise regression analysis was performed for normal irrigation and drought stress conditions (Table 3). Accordingly, the results showed that the total protein content was the sole feature selected as the most influential feature on grain yield in normal irrigation condition. The model r-square, known as coefficient of determination, reached to 89% showing a high impact of total protein content on grain yield. Also, the regression coefficient for protein content was negative (b = − 27.79) showing a decreasing effect for this feature toward grain yield under normal irrigation condition. On the other side, stepwise regression model for drought stress condition selected the content of free proline, SOD and GR activities as the influential feature on grain yield of triticale in this condition (Table 3). At the first step of the feature selection for drought stress condition, SOD activity was inserted into model showing that this feature is the most influential independent variable toward the grain yield of triticale. After that, proline and GR activity were respectively inserted into the stepwise model as the second and the third influential variables. The final model of this regression analysis explained about 81% of the total variability in grain yield under the drought stress condition.
Since the feature selection method based on regression analysis is somewhat biased to availability of any correlation within in depended variables, known as collinearity, a deeper delve over the influential features on grain yield of triticale genotypes was conducted by plotting heatmap, showing the cluster analysis of the measured features under both normal and stress conditions (Fig. 1). Under normal irrigation condition (Fig. 1a), three cluster of the features were pointed out. Subsequently, APX, H 2 O 2 , SOD, and MDH were fallen into the same cluster, while, the protein content fell into a different cluster from the yield. This result verified the result of feature selection and stepwise model, in which the protein content showed a negative slope toward the grain yield. Similarly, SOD, MDH, and APX shared the same cluster with grain yield under drought stress condition (Fig. 1b). Meanwhile, GR and proline, which both showed negative regression coefficients for grain yield in stress condition, were grouped in the same cluster and they were clearly separate from the grain yield cluster. The results of heatmap analysis for drought stress condition corroborated the results of the feature selection model, in which SOD showed a positive influence on the grain yield. Biplot analysis based on the first two components (PC1 and PC2) of principal component analysis (PCA) were applied to find out about the geometrical and structural relationship among measured features and genotypes at the same time. The overall of roughly 64% of the total variation among all features was explained by PC1 and PC2 in normal irrigation condition (Fig. 2a). This twodimensional (2D) plot showed that, under normal condition, APX, MDH, H 2 O 2 , SOD, POD, and CAT were placed in the same quarter with the grain yield, while the rest of the features were placed in another quarter. Consequently, the arrow of total protein content in normal irrigation condition had the most obtuse (widest) angle with yield arrow indicating negative relationship of this feature as compared with other variables, approving, once again, the result of the feature selection model and heatmap analysis. In a similar pattern to normal irrigation condition, around 61% of total variation among features under drought tress condition was explainable by PC1 and PC2 extracted from PCA. However, the distribution of measured features under drought stress condition was significantly wider than normal irrigation condition; the features were scattered among all four quarter of the biplot (Fig. 2b). Grain yield, MDH, SOD, and APX were distributed in the same quarter of biplot indicating the positive association among these variables. CAT, GR, and proline content were the features in the same quarter directly reverse to the grain yield arrow, showing the negative association of these feature with yield, that once more, verifies the results achieved by feature selection model and the heatmap analysis in drought stress condition.
Among all enzymatic antioxidants considered in this study i.e. SOD, POD, GR, CAT, and POD, all advanced statistical methods selected SOD as the most effective enzyme on grain yield. However, the results achieved by the appropriate statistical methods are clearly indicated by the phenotypic properties and therefore, molecular and genotypic investigations would enable us to test these results. Consequently, promotor analysis of the used enzymatic antioxidants in this study was appraised so as to find any different or similar pattern regarding their regulation mechanism to explain and/or test the actual influence of SOD in drought stress condition. Unfortunately, there is still no bioinformatic database for triticale genome sequence, so we performed the bioinformatic investigation of promotor analysis in wheat genome, its parental line that shares the same AA and BB genomes with triticale. The regulatory elements that were identified in the promotor regions of the enzymatic antioxidants along with their corresponding carrier genome sets are presented in Table 4. The results showed that the greatest number of the enzymatic antioxidant genes were distributed among A and D genome groups. Also, the results revealed a significant difference between the regulatory elements of SOD genes and the regulatory elements of the other antioxidant genes. A high number of regulatory elements in promotor regions of SOD genes were related to circadian regulation, abscisic acid production, auxin response pathways, MYB binding sites, and methyl jasmonate responsive genes that are important compounds of stress responses in plants [4] and were mostly distributed in A chromosome set. In addition, gibberellin and light responding related regulatory elements showed high distributions in SOD promotor regions, especially in A genome set. The total numbers of regulatory elements of SOD promotor regions (1441) were significantly higher than any other considered antioxidant (Table 4).

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

Differential expression analysis of SOD isozymes
The relative expression analysis based on real-time PCR method for SOD isozyme containing Mn-SOD, Cu/Zn-SOD, and Fe-SOD was carried on two selected genotypes, genotype 28 and 3, with highest genetic variability with each other. The analysis of variance for the SOD isozymes were carried out based on a factorial experiment, in which genotype and time spans after inducing drought stress were the two experimental factors. The main effect of genotype on Mn-SOD and Cu/Zn-SOD gene expression was significant for both shoot and root tissues, but it did not show a significant effect on Fe-SOD gene expression of these tissues (Fig. 4). Mn-SOD and Cu/Zn-SOD in root tissue were not significantly influenced by the main effect of time after the drought stress, whereas, these two isozymes in the shoot tissue coupled with the Fe-SOD isozyme in both root and Table 4 Results of the promotor analysis for enzymatic antioxidants in different sets of chromosomes in wheat shoot tissues showed significant responses to the effect of time spans after drought stress. Except Fe-SOD expression rate in both shoot and root tissues, the two other SOD isozymes showed no significant response to the interaction effect of genotype × time span (Fig. 4). The Mn-SOD relative expression in both triticale genotypes showed an increasing trend in response to longer duration of time after drought stress in both shoot and root tissues (Fig. 4a). The expression rate of Mn-SOD was higher in shoot than expression rate in the root tissue. Additionally, genotype 28 showed significantly higher expression rate regarding Mn-SOD isozyme in both root and shoot in comparison with genotype 3. The expression rate of Cu/Zn-SOD gene showed different patterns regarding two used triticale genotypes in both shoot and root tissues (Fig. 4b). In genotype 3, Cu/Zn-SOD expression rate in shoot increased by changing the time span from 12 h to 36 h, but there was no significant different between 36 and 72 h in this tissue. Quite the contrary, 36 h showed the lowest Cu/Zn-SOD expression rate than the two other time spans in root tissue of genotype 3. Quite the contrary, the expression rate of Cu/Zn-SOD isozyme in both root and shoot tissues showed a growing (inducing) pattern in response to the longer time span after the drought stress in genotype 28. In both shoot and root tissues of the both genotypes, Fe-SOD expression rate was reduced in response to longer time span, but the decreasing rate in was slower in genotype 28 than the other genotype (Fig. 4c). Comparison the expression rates of the three assessed isozymes of SOD gene with regard to the genotypes, showed significantly higher expression rate of Mn-SOD and Cu/Zn-SOD in genotype 28 than that of genotype 3. The relative expression of Fe-SOD, on the other hand, showed no significant difference between both genotypes.
It should be noted that the measurement of SOD activity in shoot and root tissues showed an increasing pattern in response to longer time span after stress induction in both genotypes. In other words, the longer the drought stress lasted, the higher the activity of SOD enzyme reached. Nevertheless, the SOD activity was significantly higher in genotype 28 in comparison with the genotype 3. Moreover, shoot tissue showed a higher activity of SOD enzyme than that of the root tissue.

Discussion
In this study, grain yield of triticale and all biochemical features significantly responded to the effects of environment and genotypes as the main sources of variations in the combined ANOVA. This result is indicating the differences among triticale genotypes in response to altered environments, which is alluding to the proficiency of the screening programs for introducing high potential genotypes. Thus, the estimation of genetic parameters would alter in different conditions. Consequently, heritability, PCV, and GCV of almost all features showed significantly higher values in drought stress condition than that in the normal irrigation condition which indicates that screening and breeding programs, in this set of triticale genotypes, would be more efficient in drought stress condition. On the other hand, the main aim of the current study was to find proper methods to reach out to high potential genotype to tolerate drought stress condition. Accordingly, the assessments were done for both normal and stress condition, but the results of the drought stress condition superseded over the results of the normal irrigation condition, if there was any conflict between the results of these two conditions, such as selecting the influential features on grain yield or finding the highest heritable feature. Estimation of genetic parameters for both conditions showed that SOD activity had the highest heritability and GCV under drought stress condition, while it had almost low values of these parameters under normal condition. These results indicate the importance of this enzyme in triticale genotypes under water shortage condition, because GCV shows high variabilities among genotypes and the higher screening ability for finding appropriate genotypes in this regard, also, its high heritability shows high efficiency of screening based on this feature and higher response to selection (RS), than any other feature, in the following generations. In agreement with the result of genetic variability, heatmap analysis clearly showed higher variations of colors among triticale genotypes in relation to SOD activity under drought stress condition. Furthermore, grain yield showed significantly low heritability and genetic variability in comparison with most of other biochemical features, especially SOD activity. Consequently, direct screening for grain yield would not confidently result in higher grain yield in the next generation of breeding program, but using an indirect selection criterion such as SOD for yield would be highly efficient, on condition that it shows a high collinearity and genetical association, either positive or negative, with grain yield. In concordance with heatmap results that showed similar color patterns and saturations for grain yield and SOD activity, and also the result of biplot in which SOD was the closest feature to the grain yield in drought stress condition, the results of the feature selection pointed out that SOD was the most influential feature, with positive impact, on grain yield of triticale. As a result, the capability of SOD for being used as an important indirect criterion for grain yield is verified by all used biometrical and biostatistical methods. To put it simply, genotypes that show higher SOD activity under drought stress condition probably have higher grain yield.
Thus far, there are many studies concerning the effect of environmental stresses on triticale plants, however, few studies have focused on the impacts of biochemical features on grain yield and their capability of being used as indirect criteria; and there is still no study, to the best of our knowledge, considering the heritability and genetic variability of biochemical features in triticale genotypes. In addition, there is still no related program or source code for estimation of genetic parameters e.g. heritability and genetic variability, based on combined ANOVA over different year and location; therefore, the SAS programming code that was used for estimating the biometrical parameters in the current study (supplementary materials) is the first program ever written and practically used to this regard. The results achieved by this code was successfully tested by the results of the hand calculations and the both results were in agreement with each other, except just a little difference as a result of different decimal rounding. It also should be noted that, application of heatmap in the current study and using color saturation for identifying the association between yield and biochemical feature in order for finding possible selection criteria, has never been applied by any other study before. Consequently, our results showed the efficiency of heatmap analysis in screening and breeding programs and also showing a combination of relationships within features, genotypes, and cross relationships between features and genotypes. In accordance with our results regarding the importance of SOD activity, the relationship of biochemical features and grain yield of triticale was modeled in the study of Saed-Moucheshi [45], without estimation of their heritability and other genotypic parameters, and it showed that showed that MDH, proline contents and SOD were the important features with significant impacts on triticale grin yield. Similarly, Riasat et al. [43] nominated SOD activity as the most appropriate trait to obtain higher tolerant triticale genotypes. In addition, there are several studies on different crops that are in concordance with our results regarding the importance of biochemical features as the selection criteria [1,33,44,54,58]. Saed-Moucheshi et al. [47] showed that SOD along with CAT and proline content were the effective biochemical traits toward the grain yield of barley genotypes and they introduced these traits as the candidate features capable of being used as indirect criteria in barely breeding. However, in the study of Vosough et al. [58], SOD was not able to be selected as an important feature toward grain yield of wheat landraces.
The results of genetic parameters estimation, feature selection, heatmap analysis, and biplot graph alluded to SOD as the most important antioxidant in triticale under drought stress condition. However, these results are based on phenotype and indirect possible products of gene actions, which cannot show the molecular changes inside the cells. Therefore, in order to evaluate the antioxidants at the genomic and molecular level and revealing any plausible similar or dissimilar pattern regarding their regulation mechanisms explaining the aforementioned results, the upstream sequences of the antioxidant genes were appraised by bioinformatics tools. This method through which, the upstream sequences of a gene is evaluated, is known as promotor analysis. The result of promotor analysis showed clear differences between promotor region of SOD and other antioxidant genes. The number of regulatory elements in SOD promotor region was significantly higher than the other antioxidants, mostly by twice as many regulatory elements as others. Also, most of the regulatory elements in SOD promotor region were compounds related to stress responsive pathways. The numerous regulatory elements in promotor region of SOD indicates that this gen can be regulated by various responsive pathways and is most likely more active than other antioxidant in response to changing conditions. In addition, high number of stress related elements in this region shows that this gene is highly responsive to stress condition and it may be highly increased by occurring stress condition. The result also showed a great number of MYB binding sites in the promotor region of SOD genes, that was significantly higher than other antioxidant genes. MYB elements are specifically acting in the responses to water shortage condition, verifying the importance of SOD under drought stress condition.
By combining the results obtained from field study regarding the high heritability and influence of SOD in drought stress condition and the results of promotor analysis showing a high number of stress responsive regulatory elements for this gene, the importance of SOD to increase the tolerance of plants under drought stress condition was clearly verified. SOD is one of the active antioxidant enzymes in dismutation pathways of oxygen free radicals. This enzyme orchestrates the breakdowns of the superoxide radicals (O 2 − ) into hydrogen peroxide (H 2 O 2 ) in the different organelles [20,29,46]. Since superoxide radial is the most abundant radicals in some organelles, such as chlorophyll cells [16,25,28,42], it might be said that SOD is the first defense level against occurring oxidative stress in plants [48] which is concordance with the results of the current study regarding high number of stress regulatory elements making SOD to quickly respond to the changing environment. Lin et al. [32] reported that considering the regulatory elements within promotor region of any encoding gene is one of the important strategies to study the molecular responses of plants to different environments. By occurring drought stress in the plant tissues, the plant starts to produce some signal components in accordance with the type of stress, and leads them to reach out to the highly responsive genes [41]. These responsive genes are again producing some other regulatory proteins and compounds and transfer them toward the required genes. These signaling proteins and compounds are known as cis acting regulatory elements [21]. Additionally, most part of the expression regulatory mechanism of plant genes is located within 1000 base pairs of the upstream region to the transcription initiation site [13]. Accordingly, to regulate the expression of SOD genes, a number of cis acting regulatory elements must be present in the promoter region of this gene to regulate its expression based on intracellular changes. Therefore, higher number of regulatory elements that are able to bind with the regulatory motifs in the promotor region bring about higher sensitivity of the corresponding gene to intracellular changes, such as what was identified for SOD promotor in this study. A recent study [24] on the expression pattern of SOD genes in grape wine, showed that its regulatory mechanisms could be active at the three levels of transcription, posttranscription, and translation, and their functional roles are focusing on response to stresses during plant development indicating, in concordance with our results, the high importance of SOD expression under environmental stress conditions and the effect of the regulatory elements in this pathway.
The overall results of biometrical, biostatistical, and bioinformatics analyses indicated that SOD was the most important antioxidant in response to drought stress condition. Also, combining the final results regarding clustering, biplot graph, and heatmap analysis along with the mean comparison for SOD activity and triticale grain yield led to the selection of two genotypes, i.e. genotype 28 and genotype 3, with the highest genetic distance in order to assess the deferential expression profiling of SOD genes in triticale. Considering the activity of SOD enzymes, in the two selected genotypes, showed that the activity of SOD was increased in both genotypes in response to longer time span from initial point of drought stress, but the rate of increase in genotype 28 was over threefold higher than genotype 3. Accordingly, the rate of relative expression of Mn-SOD was about threefold, and Cu/Zn-SOD was about twofold higher in genotype 28 than those in genotype 3. The expression levels of Fe-SOD showed no significant difference between the two genotypes. These results indicated that higher activity of SOD enzyme under drought stress condition is most likely due to the higher expression level of Mn-SOD and/or the Cu/Zn-SOD isozymes, while Fe-SOD is, seemingly, not significantly involved in response to the drought stress. In addition, quite the contrary to Mn-SOD and Cu/Zn-SOD, the expression pattern of Fe-SOD was decreased in response to longer time span from drought stress starting point. This result indicates that Fe-SOD is probably more responsive to non-stress signaling pathways and its expression level under drought stress is decreased so as to let the two other SOD isozymes to act and grow their expression level more rapidly. Furthermore, a comparison between Mn-SOD and Cu/Zn-SOD indicated that the expression level of Mn-SOD is continuously increased as the time span increased in both genotypes; while, the Cu/Zn-SOD showed altered pattern in response to different time spans in each genotype. Therefore, these results are possibly pointing out that Mn-SOD is a general responsive isozyme that its expression level is, somewhat, quite independent of the genetic differences between the genotypes within a species, on the other hand, Cu/Zn-SOD is more dependent on genotypic differences between the two genotypes. In addition, the result related to the expression levels of Cu/Zn-SOD made us to announce it as a genotype specific isozyme, meaning that its expression level is highly correlated with the genetic contents and background of its carrier individual and is likely more expressed in the genotypes tolerant to drought stress. Furthermore, the expression level of Mn-SOD in shoot and root tissues was roughly equal. Contrarily, the expression rate of Cu/Zn-SOD in genotype 3 was higher in shoot than that in root, while its expression level in genotype 28 was significantly lower in shoot than that in the root tissue. The first organelles that are affected by drought stress are normally those in the shoot [15,35], consequently, this leads the expression levels of Mn-SOD to be significantly higher in shoot. Higher duration of drought stress causes the stress to reach out more to the root system and involve its organelles as well, as a result, the expression level of Cu/Zn-SOD isozyme in the root tissues, mostly in tolerant genotypes, is increased. Conversely, the expression levels of Fe-SOD were continuously higher in root in comparison with the shoot tissue. This result increases the likelihood involvements of Fe-SOD isozyme in non-stressful responsive pathways, because, as mentioned before, the root organelles are less likely affected by negative impacts of drought stress than the shoot organelles.
SOD has been reported to be a key enzyme in regulating intracellular ROS levels and maintaining normal physiological conditions under oxidative stress caused by effect of water shortage condition [30,36]. The findings of a study performed by Zhao et al. [61] showed that the SOD expression level and malondialdehyde content in rice were significantly correlated with each other, in addition, the activity of SOD expression level in rice cultivars sensitive to heat stress were much lower than the resistant cultivars. Also, in agreement with our results, Mohammadi et al. [38] showed that Cu/Zn-SOD was significantly expressed in common bean plants under drought stress condition, but the rate of its induction in tolerant genotype was significantly higher than the susceptible one. In concordance with our result regarding the importance of Cu/Zn-SOD as a more specific isozyme that its expression levels are most likely changed by genetic contents of its carrier individuals, Xie et al. [60] stated that it is highly feasible that Cu/Zn-SOD play an important role in regulating total SOD activity and ROS detoxification in stressful condition. Moreover, attempts to regulate Cu/Zn-SOD levels in plant chloroplasts through the gene transfer mechanism resulted in transgenic plants with altered expression of the aimed gene. Overexpression of Cu/Zn-SOD genes has been reported to increase stress tolerance in tobacco transgenic plants [55]. In this study, the transgenic tobacco plants had about 30 times higher SOD activity than that in non-transgenic plants. In the study of Sheoran et al. [51], the expression of Mn-SOD gene in C306 and KAW3717 wheat genotypes was significantly increased in response to drought stress conditions; this study significant consistence between expression rate of Mn-SOD and higher overall enzymatic activity of SOD which is in line with the results of the current study related to the higher expression levels of Mn-SOD in drought stress condition. Consequently, some researchers, such as Awan et al. [4] Sheoran et al. [51] and Wang et al. [59], suggested that Mn-SOD is an important agent that plays an important role in response to drought tolerance, however, our results showed that Mn-SOD is possibly a general responsive isozyme that its response to drought stress is similar (or the same) to other stress conditions.

Conclusion
The overall results indicated the high possibility of screening favorable and potential genotypes to announce as new triticale cultivars as the result of great variabilities among use genotypes in all six years of the study. Additionally, a new programming code written in SAS language was successfully applied to estimate the genetic parameters based on the results of combined ANOVA for the first time since their introduction. The genetic parameters showed a high heritability and genetic variability for SOD activity under drought stress condition. The results of feature selection, heatmap analysis, and biplot clearly indicated that SOD had a significant effect on triticale production in drought stress condition. In this regard, heatmap graph, which is based on the saturation of different colors that shows the association within variables and cases (genotypes) and also their cross relationships, was successfully used for the first time in this study and it showed to be highly appropriate for screening favorable genotypes and processing the features that are proper for indirect selection in breeding programs. According to the results of promotor analysis for antioxidant enzymes, there was a significant difference between SOD and other antioxidants regarding the numbers and the types of motifs and regulatory elements available in their promotor regions. Bioinformatics study showed that the most of regulatory elements in SOD promotor regions were related to stress responsive pathways, especially drought stress as the MYB binding site. Accordingly, two genotypes with the highest genetic distance were selected to consider the differential profiling of three different isozymes of SOD genes in shoot and root organs of the triticale genotypes. The result showed that tolerant genotype (genotype 28) had significantly higher levels of expression for Mn-SOD and Cu/Zn-SOD than susceptible genotype (genotype 3), while the expression level of Fe-SOD was not significant between these two genotypes. Mn-SOD showed higher expression level in shoot than root, but the expression level of Fe-SOD was conversely higher in the root. Whereas, the expression level of Cu/Zn-SOD in tolerant genotype was higher in root while its level in sensitive genotype was lower in root, in comparison with shoot. Finally, it could be stated that Mn-SOD is a general isozyme that responses positively to stress condition in probably most of the plant species with different genetic backgrounds, while, Cu/Zn-SOD is a genotype specific isozyme that is likely more expressed in some genetic content such as the tolerant genotype of triticale in the current study. Furthermore, genotype 28 can be introduced as a favorable genotype for being selected for further assessments regarding its potentials in drought stress condition.

Methods
This study was consisting of overall three different types of practical research: the first type was consisted of field experiments for considering the phenotypic features, the second one was comprised of a greenhouse experiment for considering expression levels of SOD isozymes in triticale plants, and the third research was containing bioinformatic evaluation and promotor analysis of antioxidant genes.

Field study
The field study section of this research was consisted of 12 separate trials that were performed in six con-  Table 1.
In every year, 28 advanced elite genotypes, originated from international CIMMYT center, and two commercial cultivars (Sanabad, and juanilo) of triticale (Supplementary Table 2) were sown in mid-November (10th to 20th) in two separate trials. In each year, two separate sites in the same location were used for two different conditions containing normal irrigation and drought stress conditions. In all experiments, the irrigation schedule was determined according to the time of depleting about 40% of available soil water capacity (SWC). SWC is the difference between wilting point (WP) and field capacity (FC) based on following formula [19]: SWC: soil water capacity, FC: field capacity, WP: wilting point.
The WP and FC of the soil were accurately measured just the day before sowing date in every year in order to calculate the soil available water capacity. For measuring FC and WP the method of Samarah et al. [49] was applied in which the soil samples were weighed before and after oven-drying at 105°C for a period of 24 h. Also, for determining the irrigation schedule, the soil samples of 5 fixed experimental units at the depth of 0-40 cm were collected every day, started in the second day after each irrigation, and transferred to the Soil Science Labs and SWC was regularly determined. The drought stress treatment was implemented by withholding irrigation from starting point of heading stage until the crop harvest. In all experiments, either normal irrigation or drought stress trial, 150 kg urea fertilizer and 100 kg triple superphosphate fertilizer were applied just before the sowing of the triticale seeds. All experiments, in both normal irrigation and drought stress trials (12 overall trials), were arranged based on randomized complete block design (RB) with three replications. Each experimental plot (unit) in all experiments was consisted of four rows with 20 cm between rows and 2 m row length (1 m width by 2 m length). At harvest time, two central rows in each plot were used for measuring grain yield.

Greenhouse study
Following year to the six years of field experiments (2020), a greenhouse experiment was carried out in the greenhouse of Crop Production and Plant Breeding Department of School of Agriculture, Shiraz University, using two selected triticale genotypes consist of genotype 28 and genotype 3, in order to measure the relative expression rate of SOD isozymes (Mn-SOD, Cu/Zn-SOD, and Fe-SOD) under drought stress condition. The triticale seeds (10 seeds per pot) were sown in plastic pots containing 250 g soil. After emergence, the number of plantlets in each pot was reduce to three plantlets. In order to schedule the irrigation times in the greenhouse experiment, RWC and WP of the soil were accurately determined prior to seeds sowing. The content of water for reaching to 100% of FC for plastic pots were measured by weighing the pots (pot + 250 g soil + the content of water for 100% FC). During the experiment, 8 random pots were daily weighed to determine the content water in soil and when, the water content met 80% FC, the pots were irrigated again up to 100% FC. At four leaves stage, water shortage stress was applied by letting the water content of corresponding units meet 60% FC. At this stage, in order to obtain an accurate and precise relative gene expression, the experimental units were divided into two sets: treatment set was consisted of 18 pots (2 genotypes × 3 time-spans × 3 biological repeats) which used for water shortage stress, and control set was consisted of 6 pots (2 genotypes × 3 time-spans) that used as the control for all three replications within each sampling time for each genotype. This division was done due to the destructive manner of the sampling that took place for both shoot (the youngest leaves) and roots of plants in each experimental unit (pot). Sampling for stressed plants were done at 12 h, 36 h, and 72 h after applying the water shortage stress (after reaching 60% FC water content). Sampled tissues were frozen by liquid nitrogen immediately after their separation and transferred to the genetic lab and kept under 80°C.

Bioinformatic study
In order to investigate the sequences of the promoters regions of all genes corresponding the antioxidants activities that were measured in this study i.e. superoxide dismutase (SOD), glutathione reductase (GR), catalase (CAT), ascorbic peroxidase (APX), and peroxidase (POD), the complete sequence (CDS) of the mRNA related to these genes in wheat (because the sequence of these genes in triticale have not been identified yet) were downloaded from NCBI site (https://www.ncbi.nlm.nih. gov/nuccore). The identified accession id number for these sequences were consisting of AF387739.1 for APX, X94352.1 for Cat, D85751.1 for GR, AY857755.1 for POD, and JX398977 for SOD. After that, each sequence was uploaded onto PlantEnsemble site (https://plants. ensembl.org/index.html) and blasted against the whole genome of wheat for finding any sequence that shows over 10E-11 similarity score (e-value less than 10-11). Next, 1000 upstream base pair of each hit of PlantEnsemble site were downloaded as promotor regions for aforementioned genes. Finally, each of the downloaded sequences (promotor regions) were reviewed in Plant-Care database site (https://bioinformatics.psb.ugent.be/ webtools/plantcare/html) for distinguishing any probable regulatory elements. For further investigation and more accurate distinguishing of the regulatory elements, RSAT database (https://rsat.ulb.ac.be/rsat) was used for searching among all sequences for all genes at the same time. All distinguished regulatory elements were then imported into excel files for further investigation and analysis.

Measurements
In all field experiments of this study (6 years by 2 stress levels), 25 days after drought stress induction, the flag leaves were sampled for measuring biochemical and antioxidant-related features. The flag leaves were immediately frozen by liquid nitrogen and transferred to refrigerators having the temperature below 80°C till the time of feature measuring. The biochemical features containing content of hydrogen peroxide (H 2 O 2 ) [2], content of malondialdehyde (MDH) [22], the amount of total chlorophyll (TChl) [31], free proline (PRL) [5] and total protein (Prtn) [8] contents were measured along with the activities of antioxidant enzymes i.e. glutathione reductase (GR) [52], ascorbic peroxidase (APX) [39], catalase (CAT) [14], peroxidase (POD) [9], and superoxide dismutase (SOD) [6]. The content of H 2 O 2 and MDH were reported as micro mol per gram (μmol/g) fresh weight (FW), the amount of chlorophyll and protein as milligram per gram (mg/g) FW, and proline content as micromolar per gram (μM/g) FW. The activity of enzymatic antioxidants was reported as micromole of free radicals decomposed in mg of total protein in the sample (u/mg Prt). In addition, at the harvest time, grain yield of the two central rows in each plot were harvested and their grain yield were recorded according to kilogram per square meter (kg/m 2 ).
After sampling the youngest leaves of the plants in greenhouse study at each time span, after induction of the drought stress, the activity of SOD enzyme was measured along with the relative expression of SOD isozymes. For extracting the RNA from leaf samples, Dena Zist isolating kits were used. After the extraction, the quantity of extracted RNA was measured by nanodrop method. Following the quantity measurement of RNA, the extracted samples were treated with DNase in order to remove any DNA contamination, afterward, the positive and negative quality control (QC) of the RNA in each sample were performed; the positive QC was checked by direct electrophoresis of 3 mL extracted samples on gel (1% of 0.5 TBE buffer) and observing two separate and sharp bands on the gel representing each of 16 s and 24 s RNA in the samples, while, the negative QC was checked by observing no band on electrophoresis gels posterior to polymerase chain reaction (PCR) of samples conducted by using the primer of internal control (housekeeping: elongation factor 1) gene that was used in this study. Then, the RNA of the extracted samples was transformed into complementary DNA (cDNA) by use of Fermentas Co. kits. The produced cDNA's were then used for measuring the relative expression of three isozymes of SOD gene i.e. Mn-SOD, Cu/Zn-SOD, and Fe-SOD based on RT-PCR method. The corresponding pair of primers for each isozyme was designed in AlleleID software according to the sequence of the isozymes in wheat genome. In this process, the elongation factor 1 (EF1) was used as the internal control (housekeeping) gene. Each pair of forward (F) and backward (revers: R) primers that were used for expression analysis of Mn-SOD, Cu/Zn-SOD, Fe-SOD, EF1 are as follow: required for this study. The authors also highly appreciate the scientific and grammatical editing that were done by Dr. M.P.
Authors' contributions AS contributed directly in all trials (six years) of the study. FS helped to conduct the measurements in Shiraz location along with bioinformatics and expression analyses in the University of Kurdistan. The final analyses of data and writing the programming codes in SAS-language software were done by AS and cooperated and tested by EF, and FB. EF and FB helped in the bioinformatics section by providing and collecting the data. The trials in Zarhgan location managed mainly by MR. MR also helped the measurements and data analyses for Zarghan location. AAM managed the experiments in Sanandaj location along with the measurements of the features and writing of the manuscript. The writing of the final manuscript was carried out by AS and properly edited by AAM. All authors have read and approved the final version of the manuscript. Availability of data and materials All data generated or analyzed during this study are included in this published article [and its supplementary information files].

Declarations
Ethics approval and consent to participate