**Transcriptomic Analysis of Female Panicles Reveals Gene Expression Responses to Drought Stress in Maize (***Zea mays* **L.)**

**Shuangjie Jia 1,**†**, Hongwei Li 1,**†**, Yanping Jiang 1, Yulou Tang 1, Guoqiang Zhao 1, Yinglei Zhang 1, Shenjiao Yang 2, Husen Qiu 2, Yongchao Wang 1, Jiameng Guo 1, Qinghua Yang 1,\* and Ruixin Shao 1,\***


Received: 8 January 2020; Accepted: 14 February 2020; Published: 24 February 2020

**Abstract:** Female panicles (FPs) play an important role in the formation of yields in maize. From 40 days after sowing to the tasseling stage for summer maize, FPs are developing and sensitive to drought. However, it remains unclear how FPs respond to drought stress during FP development. In this study, FP differentiation was observed at 20 and 30 days after drought (DAD) and agronomic trait changes of maize ears were determined across three treatments, including well-watered (CK), light drought (LD), and moderate drought (MD) treatments at 20, 25, and 30 DAD. RNA-sequencing was then used to identify differentially expressed genes (DEGs) in FPs at 30 DAD. Spikelets and florets were suppressed in LD and MD treatments, suggesting that drought slows FP development and thus decreases yields. Transcriptome analysis indicated that 40, 876, and 887 DEGs were detected in LD/CK, MD/CK, and MD/LD comparisons. KEGG pathway analysis showed that 'biosynthesis of other secondary metabolites' and 'carbohydrate metabolism' were involved in the LD response, whereas 'starch and sucrose metabolism' and 'plant hormone signal transduction' played important roles in the MD response. In addition, a series of molecular cues related to development and growth were screened for their drought stress responses.

**Keywords:** transcriptome analysis; summer maize; drought; female panicle

#### **1. Introduction**

Under the influence of global warming, changes in climatic conditions are creating unusual weather phenomena worldwide, often imposing drought stress on crops [1,2]. From the agricultural perspective, drought often results in decreased crop productivity and growth [3–5], especially for cereal crops. Maize (*Zea mays* L.) is one of the most important cereal crops and has the most extensive planting area globally [6,7]. One of the most important factors limiting maize growth and development around the world is a lack of water [8–11]. Accordingly, improving tolerance of maize to drought stress is essential for achieving high and stable yields in cereal crops.

As a multidimensional stress, water limitation triggers a wide variety of plant responses; these range from responses at the physiological and biochemical levels to the molecular level [12–16]. When

external drought stimuli are perceived and captured by sensors on cell membranes, the signals are transmitted through multiple signal transduction pathways. Then, plant can regulate the expression of drought-responsive genes to protect themselves from the harmful effects of external stimuli [17,18]. The expressed products of drought-responsive genes are mainly proteins involved in signaling cascades and transcriptional regulation (such as protein kinase, protein phosphatase, and transcription factors) and functional proteins [19]. With the rapid development of high-throughput sequencing technologies, transcriptome analyses have been conducted to identify stress-mediated differences at the level of gene expression. Previous research has shown that many significantly differentially regulated genes that were associated with drought tolerance are induced in different organs of maize plants [20–26]. For example, 249 and 3000 differentially expressed genes (DEGs) were involved in root tissues after 6 h of light and severe drought stress, respectively [23]. In leaves, a total of 619 DEGs and 126 transcripts had their expression levels altered by drought stress at flowering time [24]. In tassels, 1902 DEGs were found after 5–7 days of drought stress [25]. In young ears, a total of 1825 DEGs were identified on the 5th day of drought stress at the V9 stage [26].

The panicle stage (from jointing to flowering) is the key stage for panicle differentiation and development in maize, the number of rows per ear and the number of grains per row are dependent on spikelet and floret differentiation at this time [27]. Therefore, the growth and development of female panicles (FPs) play an important role in the formation of maize yields. Although great advances in understanding differentiation of FPs and how drought stress affects genes transcription in FPs have been achieved in the past few decades [26,28–31], so far, progress in understanding the general molecular basis of FP development in response to long-term drought stress across the panicle stage has not been reported.

Accordingly, in this study, maize inbred line PH6WC (6WC) was used as drought-sensitive experimental material [32], soil water was controlled by means of drip irrigation for 30 days, and the gene expression dynamics of developing FPs at 30 days after drought (DAD) were investigated using transcriptomic analysis. The DEGs were identified and assigned to functional categories to reveal the various metabolic pathways in FPs that are involved in responses to long-term drought stress at different levels. Furthermore, differences in transcription factors between treatments were also analyzed. Overall, the exploration and function prediction of drought-response genes in maize FPs represent an efficient approach to improving the molecular breeding of drought-resistance maize cultivars.

#### **2. Materials and Methods**

#### *2.1. Plant Material and Growth Conditions*

Field experiments were conducted at field experiment stations (34◦31 N, 115◦35 E, 50.7 m above sea level) during the maize growing season (June–October, 2018) in Shangqiu (Henan, China). The maize inbred line 6WC was grown in 9 experimental plots (each plot was 2 m wide, 3.3 m long, 1.8 m deep) which were under a movable awning and filled with luvo-aquic soils, with a 20 cm sand filter layer at the bottom. Maize were planted into four rows, with 40 cm between rows in each plot. Two to three seeds were sown at each acupoint, with subsequent thinning to one seedling conducted at the trifoliate stage (V3). The final stand density was 8 plants m<sup>−</sup>2. During the jointing and tasseling stages, topdressing fertilizer was applied, and weeds, insects, and diseases were controlled throughout the experiment. The top soil (0–40 cm layer) had a pH (water) of 7.3, mean mineral P content of 3.24 g/kg, and inorganic N at sowing of 3.60 g/kg. The average daily maximum and minimum temperatures of the field experiment during the trial were 32.98 ◦C and 20.71 ◦C, respectively.

#### *2.2. Drought Stress Treatments*

During FP development, soil moisture was controlled by means of drip irrigation at 80 ± 5% of the field water capacity (FWC) (well-watered, CK), 60 ± 5% of FWC (light drought, LD), and 45 ± 5% of FWC (moderate drought, MD). The meter was checked every morning and evening throughout

the growth period to guide adjustments of the soil moisture. When the treated soil moisture dropped towards its lower limit, moderate drip irrigation was carried out until the upper limit level was reached, and the irrigation volume was measured by a water meter. After 30 DAD, the drought treatment plots were rehydrated to the CK level. Other field management measures reflected standard field management practices.

#### *2.3. Measurement of Morphology and Microscopic Observation of Female Panicles*

Plant height was measured from the ground to the top of the leaves in their natural growth state at 20, 25, and 30 DAD. The length and width of all the green leaves were measured by ruler in order to calculate leaf area, and the leaf area index (LAI) was determined according to this method [33]. Dry matter accumulation in stalks, leaves, tassels and ears of maize plants were measured at 20, 25, and 30 DAD. After 30 min of defoliation at 105 ◦C, dry weight was determined after being dried at 75 ◦C until a constant weight was reached. The percentage of drought limitation was calculated as (T2−T1)/T1. Here, T2 was shoot dry matter under the MD or LD treatment, while T1 was shoot dry matter under the control or LD treatment [34]. FPs were dissected with a dissecting needle at 20 and 27 DAD and analyzed under a stereomicroscope (Guanpujia, SMZ-B2, Beijing, China) to observe the differentiation of developing female inflorescences. There were three biological replicates in each group.

#### *2.4. RNA Isolation and Illumina Sequencing*

Three plants with consistent growth were selected from each treatment at 30 DAD, FP bracts were sampled, and the upper, middle, and lower parts of the ears were mixed evenly and then frozen at −80 ◦C prior to RNA-sequencing (RNA-seq) analysis. Total RNA was extracted using the mirVana miRNA Isolation Kit (Ambion, Inc., Austin, TX, USA). RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Shanghai OE Biotech (Shanghai, China) conducted RNA-Seq library construction and high-throughput sequencing based on total RNA from the female inflorescence. The libraries were sequenced on the Illumina sequencing platform (HiSeq 2500l, Illumina, San Diego, CA, USA) with 125 bp paired-end reads. The raw RNA-seq data have been uploaded to NCBI SRA (BioProject ID: PRJNA604094).

#### *2.5. Read Mapping and Di*ff*erential Expression Analysis*

Base calling was conducted using the raw image data generated by sequencing to obtain sequence data, and the called raw data (raw reads) were stored in fastq format. Raw data (raw reads) were processed using Trimmomatic [35]. The reads containing poly-N runs and low-quality reads were removed to obtain the clean reads. Then, the clean reads were mapped to the reference (NCBI\_B73\_v4) genome [36] using HISAT2 (version 2.2.1.0) [37]. The Fragments Per kb Per Million Reads (FPKM) values were calculated using cufflinks (version 2.2.1) [38,39], followed by differential expression analysis using DESeq (version1.18.0) [40]. Genes with |fold change| > 2 and *p* < 0.05 were identified as differentially expressed genes with *p* presented as raw *p*-values rather than FDR adjusted *p*-values.

#### *2.6. GO and KEGG Enrichment Analysis*

The significantly expressed Gene Ontology (GO) terms were selected by GO enrichment analysis according to the GO database (http://geneontology.org/). The differences in the frequency of assignment of GO terms in the DEG set were compared with the expressed genes in the CK, LD, and MD samples (*p* < 0.05). Functional groups encompassing DEGs were identified based on GO analysis, and pathway analysis was conducted according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/), with manual reannotation based on several databases and a literature search.

#### *2.7. Di*ff*erential Expression Verification by Quantitative Real-Time PCR (qRT-PCR)*

Transcriptome sequencing data were validated by qRT-PCR. Total RNA was reverse transcribed using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TRANS, Beijing, China). The qRT-PCR experimental methods used HiScript Q RT SuperMix for qPCR (Vazyme, Nanjing, China). The primer sequences were designed using Primer 5 and are listed in Supplementary Materials Table S1. The relative quantification 2−ΔΔCt method was used to calculate the expression level of target genes in different treatments.

#### *2.8. Statistical Analysis*

Data collation and graphic rendering were conducted with SIGMAPLOT 14.0 (Systat Software Inc., San Jose, CA, USA), Microsoft Excel, and Microsoft PowerPoint 2016 software (Microsoft Corp., Redmond, WA, USA). All data are expressed as the mean ± SD value from three independent experiments unless otherwise stated. Data were analyzed by one-way ANOVA using Duncan's multiple range test at a *p* < 0.05 significance threshold in SPSS (IBM Corp., Armonk, NY, USA).

#### **3. Results**

#### *3.1. Female Panicle Development, Phenotypic Change and Yield Components*

Plant height, LAI, dry matter accumulation, and the percentage of drought limitation were measured at 20, 25, and 30 DAD. For example, plant height was reduced by 17.63%, 17.01%, and 16.44% under the LD treatment compared with CK, respectively; furthermore, MD significantly decreased plant height by 25.86%, 26.34%, and 29.73% (Figure 1a), respectively. Leaf area index was significantly decreased under MD at 20, 25, and 30 DAD compared with CK plants, but was not significantly changed at 20 and 30 DAD under LD (Figure 1b). For dry matter accumulation, MD significantly decreased the shoot dry matter at 20, 25, and 30 DAD compared with CK plants, but it was not significantly reduced in the LD vs. CK comparison at 20 and 30 DAD (Figure 1c). In addition, the percentage of drought limitation had the highest absolute values at 20, 25, and 30 DAD after MD treatment. However, the percentage of drought limitation was not significantly affected by LD treatment at 20 and 30 DAD compared with CK (Figure 1d).

From 20 to 30 DAD, according to the book *Corn Growth and Development* [41], FPs may be in the process of differentiation. To explore the responses of FPs to drought stress, maize plants were dissected. The length and diameter of FPs were significantly decreased under the MD treatment at 20 and at 30 DAD; silk was seen in CK and LD plants, but there was no floral differentiation MD plants (Figure 2a). Proportion of dry matter in FP was decreased under LD and MD treatments (Figure 2b). To investigate the effect of drought stress on development and number of mature ears, ear size, and dry matter, yield components were determined (Figure 2c, Supplementary Materials Table S2), which showed that LD and MD treatments significantly decreased ear size and increased the bald tip length. For this reason, the numbers of rows and kernels were reduced by 14.00% and 29.00% under the LD treatment and 19.00% and 43.00%, respectively, under the MD treatment. Therefore, drought resulted in great losses in grain yield of 32.00% and 35.00% under the LD and MD treatments, respectively (Figure 2c, Supplementary Materials Table S2).

**Figure 1.** Agronomic traits changes of maize plants in response to different drought stress. Effect of drought stress on (**a**) plant height; (**b**) leaf area index; (**c**) shoot dry matter; and (**d**) the percentage of drought limitation at 20, 25, and 30 days after drought (DAD). Values are the means of the replicates ± sd. Different lowercase letters and \* symbols indicate statistical significance of differences at a *p* < 0.05 level. There were three treatments, including well-watered (CK), light drought (LD), and moderate drought (MD), and five biological replicates were sampled for each treatment.

**Figure 2.** Development of female inflorescence and agronomic traits changes of ears in response to drought stress. (**a**) Contrasting sizes and differentiation of the control and drought-treated plants at 20 and 30 DAD. (**b**) Proportion of dry matter in immature ears at 30 DAD. (**c**) Radar chart showing changes in yield traits for mature maize plants grown under well-watered (CK), light drought (LD), and moderate drought (MD) conditions. (**b**) and (**c**) were calculated with 15 and 10 biological replicates sampled for each treatment.

#### *3.2. Overview of RNA Sequencing and Mapping*

A total of 49.42 million raw reads were obtained from PH6WC transcriptome libraries (Supplementary Materials Table S3). More than 96.72% of them (47.80 million clean reads) remained after discarding low-quality reads and reads containing adaptor sequences, which were then used for downstream analyses. The clean reads were mapped to the B73 reference genome (ZmB73\_RefGen\_v4). Overall, 94.34–94.52% of clean reads from nine samples were mapped onto the reference genome (Supplementary Materials Table S3). On average, approximately 43.93 (91.29%), 43.73 (91.13%), and 43.70 (91.17%) million reads from the CK, LD, and MD treatments, respectively, were uniquely mapped onto the reference genome.

Compared with the CK treatment, only 40 genes (log2 foldchange > 1 and *p* < 0.05), including nine up-regulated and 31 down-regulated genes, showed significantly differential expression in the LD treatment, and a total of 212 up-regulated and 664 down-regulated genes were identified in the MD treatment. A total of 887 DEGs, including 208 up-regulated and 679 down-regulated genes, were identified in the MD versus LD comparison (Figure 3a). A Venn diagram of the DEGs illustrated that there were 10 common genes that appeared in the LD vs. CK and MD vs. CK comparisons, five genes shared between the LD vs. CK and MD vs. LD comparisons, and 565 genes shared between the MD vs. LD and MD vs. CK comparisons. However, there were no DEGs commonly expressed in all three comparisons (Figure 3b).

**Figure 3.** Identification and characterization of differentially expressed genes (DEGs) between the drought treatment and control plants. (**a**) The number of DEGs in three comparison groups. (**b**) A Venn diagram comparison summarizing overlaps in differentially expressed genes among the three comparisons.

#### *3.3. Gene Expression Validation by qRT-PCR*

To investigate the changes in gene expression at the mRNA level, eight randomly selected genes and three specific genes *cuc2* (LOC103629107), *TE1* (LOC541683), *DLF1* (LOC100037791) were analyzed using quantitative real-time RT-PCR for validation of RNA-seq. The level of expression of the genes amplified is shown in Figure 4.

**Figure 4.** The expression patterns of eleven genes in female panicle tissues under well-watered (CK), light drought (LD), and moderate drought (MD) conditions by qRT-PCR. Values are the mean ± SD of three independent experiments. Maize β-actin expression was used as a control.

Raw data were compared to transcriptomics data (Supplementary Materials Table S4), which closely resembled each other, validating the differential expression of the genes identified as being under drought stress.

#### *3.4. GO Annotation and Enrichment*

A total of 26, 629, and 630 DEGs were assigned by GO analysis conducted based on the genes from the LD vs. CK, MD vs. LD, and MD vs. LD comparisons, respectively. The most significantly regulated 20 terms among biological processes from the MD vs. CK and MD vs. LD comparison genes are shown in Figure 5a,b, but there were no significant terms (i.e., terms with gene number > 2 and *p* < 0.05) resulting from the LD vs. CK comparison. The up-regulated terms from the MD vs. CK comparison are involved in "regulation of timing of plant organ formation," "developmental process," and "regulation of cell proliferation." The down-regulated terms "response to water deprivation" and "post-embryonic plant morphogenesis" were also enriched. The up-regulated terms from the MD vs. LD comparison are "regulation of timing of plant organ formation," "regulation of cell proliferation," and "gibberellin biosynthetic process." The down-regulated terms "reductive pentose-phosphate cycle," "phosphate ion homeostasis," and "ethylene-activated signaling pathway" were enriched among genes from the MD vs. LD comparison. The genes associated with GO terms related to development, growth, and responses to stimulus were also significantly different among the three comparisons. These genes that changed in their levels of transcriptional expression were basically the same in the MD vs. CK and MD vs. LD comparison groups, but lower in expression differences in the LD vs. CK comparison group (Figure 5c,d, Supplementary Materials Table S5).

**Figure 5.** Top 20 biological processes enriched by the up-and down-regulated genes in the (**a**) MD vs. CK and (**b**) MD vs. LD comparisons. Expression pattern of the differentially expressed genes associated with (**c**) development progress and (**d**) growth in the three comparisons. Colors indicate the log2 fold change values. Red indicates up-regulation, and green indicates down-regulation in that comparison.

#### *3.5. Metabolic Pathways Related to Soil Drought Stress*

To further characterize genes affected by drought stress, we performed a KEGG pathway classification analysis to identify functional enrichment of DEGs. Thus, 8, 72, and 74 terms were significantly enriched in the transcriptome profile comparisons of LD vs. CK, MD vs. CK, and MD vs. LD groups (Supplementary Materials Table S6). The significant differences in the top 20 enriched KEGG pathways in the MD vs. CK and MD vs. LD comparisons are shown in Figure 6. In the

MD vs. CK comparison, genes associated with the pathway "Starch and sucrose metabolism" were most enriched followed by those associated with "plant hormone signal transduction" (Figure 6a, Supplementary Materials Table S7). In the MD vs. LD group, "Plant hormone signal transduction," "Starch and sucrose metabolism," and "Glycosphingolipid biosynthesis" were the three most enriched terms (Figure 6a, Supplementary Materials Table S7).

**Figure 6.** Top 20 enriched KEGG pathways in the (**a**) "MD vs. CK" and (**b**) "MD vs. LD" comparisons. Pathway entries with the corresponding number of genes (among those pathways with more than two genes) are shown, and the corresponding -log10 *p*-value of each entry is sorted in descending order. The number of DEGs in each pathway is positively related to the size of plot, and the *p*-values shown in red are more significant.

#### **4. Discussion**

#### *4.1. Responses of Plant Growth and Female Panicle Di*ff*erentiation to Soil Drought Stress*

The panicle stage is the most important productive stage in corn development, and soil drought stress in this stage can affect the plant growth rate, prolong the growth processes of the panicle stage, hinder the normal differentiation and development of ears, and ultimately lead to decreased crop seed setting rates and yields [42–46]. Further, the drought response depends on the time and intensity of water loss as well as the developmental stage [47,48]. In this study, LD and MD treatments compared with the CK treatment decreased green leaf area and significantly suppressed shoot dry matter accumulation over the prolonged drought treatment, and relative to the LD treatment, the MD treatment affected plant growth much more (Figure 1), which is consistent with previous research by Boonjung et al. [49].

FPs are the precursor to maize ears, and the differentiation and development of FPs mainly occurs from the V9 to VT stages, which include growth cone extension (V9), spikelet differentiation (V11), floret differentiation (V12), and formation of the sexual organs (VT). Developing organs are sensitive to drought, especially during their early phases [50]. When drought occurs between the V9 and VT stages, how does the degree of drought affect the formation of ears? In our study, spikelet and floral differentiation, as observed under stereomicroscope, were significantly inhibited and thus delayed by soil drought (Figure 2). Some studies have shown that the number of kernel rows is determined at the spikelet differentiation stage, and the floret differentiation period is the key period that affects grain number [51–53]. Here, mature ears in the MD and LD treatments were much shorter and thinner than those in the CK treatment; in addition, the bald tip length and number of unfilled grains were both increased under drought treatments. As indicated above, drought affected grain yields as well (Figure 2c).

#### *4.2. Genes Involved in Development and Growth in Response to Soil Drought Stress*

Drought treatments affected the expression of genes associated with development and growth of the inflorescence (Figure 4). *Terminal ear 1* (*te1*) maize mutants have shortened internodes, abnormal phyllotaxy, leaf pattern defects, and partial feminization of tassels [54]. Similarly, *cup-shaped cotyledon 2* (*cuc2)* mutants have been reported to have abnormalities in the regulation of the shoot meristem boundary and formation and subsequent development [55,56]. Here, *cuc2* were up-regulated under moderate drought stress (Figure 4, Supplementary Table S5), combined with developmental change (Figure 2a), implying that the differential expression of the gene under MD treatments may be related to the mature delay of the FPs tissue. *DLF1* was also up-regulated under MD stress (Figure 4, Supplementary Materials Table S5), suggesting that the trans-activator protein encoded by, this gene plays an important role in the signal transduction pathway and the regulation of plant growth at the FP development stage [57].

#### *4.3. Genes Involved in Auxin Signaling in Response to Soil Drought Stress*

Auxin is an important phytohormone that is closely related to plant resistance to adverse environmental conditions, and it can induce rapid and transient expression of some genes, including auxin response factor genes (ARF) and primary auxin response genes (Aux/IAA, GH3, SAUR and LBD); the protein products of these genes can specifically bind to ARFs to activate or inhibit downstream gene expression under drought [58–61]. In the current study, auxin signaling genes were involved in the response to drought, as the expression of IAA-conjugating genes (GH3) was up-regulated, and the expression levels of auxin biosynthesis genes were down-regulated after MD stress (Figure 7), leading to the reduction of auxin levels (Supplementary Materials Figure S1a). This implies that drought improves GH3 transcription to help maintain endogenous auxin at an appropriate level in maize [62,63]. However, when the concentration of auxin increases, auxin combines with transport inhibitor response 1 (TIR1), causing Aux/IAA ubiquitination and degradation; then, ARF is released, which further activates the expression of small auxin-up RNA (SAUR) genes [64]. SAUR genes are early auxin-responsive genes involved in plant growth, and the SAUR family regulates a series of cellular, physiological, and developmental processes in response to environmental signals [65–67]. In our study, three SAUR genes were down-regulated under the MD treatment (Figure 7), which may explain MD-induced inhibition of maize growth.

**Figure 7.** Genes involved in auxin plant hormone signal transduction pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG). Differentially expressed genes involved in AUX/IAA, GH3, and SAUR were shown by heat-maps, and the number was calculated with log2foldchange in three comparisons. Color of heat-maps represented different fold change. Yellow box means involved significantly differentially expressed genes were mainly up-regulated, and green box means down-regulated.

#### *4.4. Reactive Oxygen Scavenging System and Ion Channel in Response to Soil Drought Stress*

Limited water supply enhances the production of reactive oxygen species (ROS) [68,69], and plants are protected by glutathione *S*-transferase (GST) and other antioxidant enzymes scavenging excessive ROS from the damage caused by ROS [70]. This is because GST comprises a large superfamily of multifunctional protein and participates in ascorbic acid (AsA)/glutathione (GSH) cycling pathways [68]. Here, probable glutathione *S*-transferase *GST12* was significantly down-regulated under the MD treatment compared with the CK treatment. However, *GSTU6* (LOC103637303) and GST activity was up-regulated under MD stress (Supplementary Materials Table S8 and Figure S1b), which may scavenge ROS and protect both plant cell membrane structure and protein activity [71,72], implying that GST is involved in responding to drought stress [73,74].

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/10/2/313/s1, Table S1: The primer sequences for qRT-PCR, Table S2: Ear and grain yield traits at harvest time under CK, LD and MD treatments, Table S3: Number of reads based on RNA-Seq data of CK, LD and MD treatments, Table S4: RNA-Seq expression levels of the eight genes for qRT-PCR, Table S5: Expression pattern of the differently expressed genes about development progress (A) and growth (B) under CK, LD and MD treatments, Table S6: Number of enriched KEGG pathways terms and DEGs in three comparisons (LD vs.CK, MD vs. CK; MD vs. LD), Table S7: Differentially expressed genes (DEGs) in the top 20 enriched KEGG pathways in two comparisons (MD vs. CK; MD vs. LD), Table S8: Genes significantly enriched in Glutathione metabolism in three comparisons (LD vs.CK, MD vs. CK; MD vs. LD), Figure S1: Effect of drought stress treatments on the content of IAA (a) and the activities of GST (b) in female panicles. CK, well-watered; LD, light drought; MD, moderate drought, five biological replicates were sampled for each treatment.

**Author Contributions:** S.J. did the experiment and wrote the manuscript. H.L. analyzed the data. Y.J., Y.T., G.Z. and Y.Z. did a part of experiment. S.Y., H.Q., Y.W. and J.G. gave some good suggestions on the manuscript. Q.Y. revised the manuscript. R.S. designed this experiment, revised the manuscript and provided the funding. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Open Foundation of CAAS/Key Laboratory of Crop Water Use and Regulation, Ministry of Agriculture (FIRI2019-02-0103), Open Foundation of State Key Laboratory of Crop Biology (2019KF03), Scientific and Technological Innovation Talents in Colleges and Universities in Henan (20HASTIT036) and Central Public-interest Scientific Institution Basal Research Fund (Farmland Irrigation Research Institute, CAAS) (FIRI2017-19).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Comparison of Biochemical, Anatomical, Morphological, and Physiological Responses to Salinity Stress in Wheat and Barley Genotypes Deferring in Salinity Tolerance**

#### **Muhammad Zeeshan 1, Meiqin Lu 3, Shafaque Sehar 1, Paul Holford <sup>4</sup> and Feibo Wu 1,2,\***


Received: 27 November 2019; Accepted: 9 January 2020; Published: 15 January 2020

**Abstract:** A greenhouse hydroponic experiment was performed using salt-tolerant (cv. Suntop) and -sensitive (Sunmate) wheat cultivars and a salt-tolerant barley cv. CM72 to evaluate how cultivar and species differ in response to salinity stress. Results showed that wheat cv. Suntop performed high tolerance to salinity, being similar tolerance to salinity with CM72, compared with cv. Sunmate. Similar to CM72, Suntop recorded less salinity induced increase in malondialdehyde (MDA) accumulation and less reduction in plant height, net photosynthetic rate (Pn), chlorophyll content, and biomass than in sensitive wheat cv. Sunmate. Significant time-course and cultivar-dependent changes were observed in the activities of antioxidant enzymes such as superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX), and glutathione reductase (GR) in roots and leaves after salinity treatment. Higher activities were found in CM72 and Suntop compared to Sunmate. Furthermore, a clear modification was observed in leaf and root ultrastructure after NaCl treatment with more obvious changes in the sensitive wheat cv. Sunmate, rather than in CM72 and Suntop. Although differences were observed between CM72 and Suntop in the growth and biochemical traits assessed and modified by salt stress, the differences were negligible in comparison with the general response to the salt stress of sensitive wheat cv. Sunmate. In addition, salinity stress induced an increase in the Na<sup>+</sup> and Na+/K<sup>+</sup> ratio but a reduction in K<sup>+</sup> concentrations, most prominently in Sunmate and followed by Suntop and CM72.

**Keywords:** antioxidants; ultrastructure; osmotic stress; salinity; wheat; barley

#### **1. Introduction**

Saline soils are a major problem in many countries with the Environment Program of United Nations estimating that of the 9–34% of the world's irrigated land is adversely affected by salinity [1]. Salinity can kill plants and other soil organisms and is referred to as a "silent killer" in some regions or as "white death" in others as it invokes images of a lifeless, shining land studded with dead trees. Approximately 32 million ha of dry lands [2] and 60 million ha of irrigated land [3] are affected by human-induced soil salinization, and it is well documented that salinity is one of the most severe environmental stresses hampering crop production [4,5]. At high electrical conductivity (EC) resulting

from salinization, crop yields can decline drastically rendering crop cultivation no longer profitable and making soil amendments inevitable [6]. World agriculture needs to feed about 2.3 billion people globally by 2050 [7]; thus, it is imperative to understand the mechanisms associated with tolerance to salinity so that breeding programs and agronomic practices can be put in place that will allow production to meet this increasing demand [8].

Saline soils limit plant growth due to osmotic stress, ionic toxicity, and a reduced ability to take up essential minerals [9,10]. In extreme cases, root cells may lose water instead of absorbing it due to the hyperosmotic pressure of the soil solution. Water deficits affect a cascade of physical, signaling, gene expression, biochemical, and physiological pathways and processes, resulting in decreased cell elongation, wilting, and, ultimately, plant death; these harmful effects of salinity can be considered as water-deficit effects [3,11,12]. In saline soils, NaCl comprises 50–80% of the total soluble salts [13] causing elevated, and potentially toxic, concentrations of Na<sup>+</sup> and/or Cl<sup>−</sup> in the plant. These ions affect many enzymes or cellular functions such as photosynthesis signaling systems [14–16]. In addition, because of their physicochemical similarities and a shared transport system, the Na<sup>+</sup> in the soil solution of saline soils competes for uptake with K<sup>+</sup> [17] and can lead to K<sup>+</sup> deficiency [18,19]. The induced K<sup>+</sup> deficiency inhibits growth because it plays a critical role in maintaining cell turgor, membrane potentials, and enzyme activities.

As a consequence of the primary effects of salinity described above, secondary stresses such as oxidative stress often occurs due to an overproduction of reactive oxygen species (ROS) [20]. These ROS cause lipid peroxidation leading to increased membrane fluidity and permeability [21,22], the denaturation of functional and structural proteins [23], and can affect nucleic acids through base modifications, induce inter- and intra-strand crosslinks, crosslinks with proteins as well as creating strand breaks [24]. However, plants have developed comprehensive internal resistance systems to combat the outcomes of ROS that are comprised of enzymatic as well as non-enzymatic antioxidants [25]. ROS-scavenging enzymes include those that are playing a direct role in the processing of ROS such as superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT), and those glutathione reductase (GR) and ascorbate peroxidases (APX) that mediate in the reaction cycle of antioxidant chemicals such as glutathione (GSH) and ascorbic acid (AsA) [26–30]. The other half of the antioxidant machinery includes nonenzymatic antioxidants comprising of ascorbic acid, phenolic compounds (flavonoids, anthocyanins), α-tocopherol, carotenoids, and amino acid cum osmolyte proline. Besides the synthesis and modulation of osmolytes some phytohormones and regulatory molecules also play prominent role in triggering salinity tolerance effector molecules [31].

Barley and wheat have different salt tolerances capacities and are grown as major grain crops in both saline and non-saline soils [32]. Previous studies have focused on salinity stress in either barley or wheat alone, with little inter-specific comparison. Thus, this study is the first to compare the mechanisms that confer salinity tolerance in these two species. We aimed to explore the similarities or differences in their physiological mechanisms upon exposure to salt stress. We also hypothesized that there may be species-specific mechanisms that can be co-related with the salt sensitivity of wheat or the tolerance of barley. Thus, this research can enhance our understanding of holistic salinity tolerance mechanisms and will aid in the breeding of salt-tolerant crops.

#### **2. Material and Methods**

#### *2.1. Plant Material and Growth Conditions*

A shade house hydroponic experiment was carried out on the Zijingang Campus, Zhejiang University, China. Two wheat (*Triticum aestivum* L.) cv. Suntop (salt-tolerant) and Sunmate (salt-sensitive) and a salt-tolerant barley (*Hordeum vulgare* L.) cv. CM72, were used in the experiment. Suntop and Sunmate are two high yielding Australian Prime Hard varieties bred by Australian Grain Technologies (AGT, Narrabri, Australia. Although they were derived from the same cross, Suntop and Sunmate differ significantly in salt tolerance. Seeds of each cultivar were disinfected in 2% (*v*/*v*)

H2O2 then washed thoroughly using double distilled water (ddH2O). The seeds were germinated on filter paper in germination boxes in a plant growth chamber (23/18 ◦C, day/night) in darkness for 3 days and incubated for further 4 days in the light. The uniform 7-day-old seedlings of each cultivar were selected and transplanted into 5 L pots in a hydroponic solution containing 4.5 L basic nutrient solution (BNS) as described by [33,34], with continuous aeration using air pumps. Each container was covered with a polystyrene plate with 7 evenly spaced holes (2 plants per hole) and placed in a greenhouse with natural light and a temperature of 20 ± 2 ◦C/day and 15 ± 2 ◦C/night. At the two-leaf stage, plants were treated with 100 mM NaCl; non-NaCl treated plants were used as controls (BNS). The solution pH was adjusted to 5.8 with NaOH or HCl, as required. The experiment was arranged in a randomized block design with four replications. Plants were sampled at 1, 5, 10, and 15 days after treatment (DAT) for time-course analysis of the salt treatments. For morphological and physiological analyses, plants were harvested 25 DAT and either analyzed immediately or frozen in liquid nitrogen and stored at −80 ◦C for further analysis.

#### *2.2. Measurement of Growth Traits and Mineral Concentrations*

Shoot height, root length, and fresh weights were determined 25 DAT (days after treatment) and then the samples were separated into shoot and root and perfectly washed with ddH2O to eliminate any foreign material. Samples from each treatment with 4 biological replicates were oven dried at 75 ◦C for three days and subsequently, the dry weight of the roots and shoots were determined in gram. Later each dried sample was weighed (about 0.2 g), ground, and made into ashes by heating the samples at 550 ◦C for half a day. Before dilution with ddH2O, the ashes were digested in 30% HNO3 and then Na<sup>+</sup> and K<sup>+</sup> concentration were quantified using flame atomic absorption spectrometry (SHIMADZU AA-6300, Columbia, Maryland, USA) [35].

#### *2.3. Measurement of Photosynthesis Parameters, Chlorophyll Contents, and Chlorophyll Fluorescence*

Intact, second fully-expanded leaves from the apex were used to measure relative chlorophyll content with the help of a handheld chlorophyll meter (Minolta SPAD-502, Tokyo, Japan) according to Wu et al. [36]. Three measurements were recorded from each leaf and averaged. The gas exchange parameters (i.e., photosynthetic rates (Pn), intercellular CO2 concentrations (Ci), stomatal conductance (Gs), and transpiration rates (Tr)) were measured on a bright sunny day between 9 a.m. to 11 a.m. using a Li-Cor-6400 portable photosynthesis system (Li-Cor, Lincoln, NE, USA).

Chlorophyll fluorescence (*Fv*/*Fm*) was measured at 25 DAT according to Genty et al. [37]. Both treated and control plants were shifted to an experimental room, kept in the dark for 25 min, flag leaf were cut for the determination of chlorophyll fluorescence using a pulse-modulated chlorophyll fluorimeter using IMAGING-PAM (Walz; Effeltrich, Germany) image processing software. Fluorescence values observed comprised of *Fo*, initial/minimal fluorescence, *Fm*, the maximal fluorescence value, and *Fv*/*Fm*, the maximum quantum yield of PSII photochemistry. The data were noted at five different points at 40, 70, 120, 150, and 180 mm from leaf tips.

#### *2.4. Lipid Peroxidation and Antioxidative Enzyme Activity Assay*

The roots and second leaf from the apex were sampled at 1, 5, 10, and 15 DAT. Lipid peroxidation was measured in the tissues and expressed as malondialdehyde (MDA) content using the TBA (thiobarbituric acid) method according to the Wu et al. [34]. The activity-specific and non-specific absorbance was determined at 532 and 600 nm, respectively. Enzymatic antioxidants activity was determined as described by Leul and Zhou [38]. Briefly, 0.2 g of frozen leaf and root plant tissue were ground in a pre-chilled mortar and pestle and homogenized in 2 mL of 1 M Tris buffer (pH 8). Later, the samples were briefly centrifuged at 10,000× g at 4 ◦C for about 15 min and the supernatants were used for the following assays. The activity assay of superoxide dismutase (SOD, EC 1.15.1.1), peroxidase (POD, EC 1.11.1.7), and catalysis (CAT, EC 1.11.1.6) were recorded according to Wu et al. [34]. Ascorbate peroxidase (EC 1.11.1.11) activity was determined at 290 nm using ascorbate (AsA)

as a substrate and 2.8 mM cm−<sup>1</sup> as an extinction co-efficient [39], while Jiang and Zhang [40] methods were used to determine the activity of glutathione reductase (GR, EC. 1.6.4.2).

#### *2.5. Cell Ultrastructure*

For transmission electron-microscopy, fresh roots (about 2–3 mm in length) and leaf pieces (about 1 mm2) without veins were hand sectioned and treated with 100 mM sodium phosphate buffer (PBS; pH 7.0) containing 2.5% (*v*/*v*) glutaraldehyde and placed overnight at 4 ◦C, then briefly washed; this step was performed 3 times with the same buffer. Each sample was treated for 60 min with 1% osmium tetroxide OsO4 (*v*/*v*) followed by washing with PBS (sodium phosphate buffer) further 3 times. Thereafter, the leaf and root samples were dehydrated with a diluted series of ethanol (50%, 60%, 70%, 80%, 90%, 95%, and 100%) for about 20 min in each solution, later all the samples were dried for 20 min in concentrated acetone. Finally, ultrathin sections (80 nm) were cut and affixed to copper grids for study using transmission electron microscopy (TEM 1230EX, JEOL, Akishima, Tokyo, Japan) at 60 kV.

#### *2.6. Statistical Analysis*

Statistical analysis was performed with the Data Processing System statistical software package [41] using ANOVA followed by Duncan's Multiple Range Tests (DMRT) to evaluate significant treatment effect at a significance level of *p* < 0.05. Origin Pro (Version 8.0, Origin lab corporation, Wellesley Hills, Wellesley, MA, USA) was used to prepare graphs.

#### **3. Results**

#### *3.1. Plant Growth Parameters*

Salt inhibited the growth of the barley and wheat plants, with the treated plants showing wilting, necrosis and chlorosis (Figure 1A). Salt damage was most severe in the wheat cv., Sunmate, while in the other cultivars, the damage was less pronounced. Salinity stress significantly (*p* < 0.05) reduced plant biomass in the wheat and barley cultivars (Figure 1B). In comparison to the other two cultivars, the effects of salt stress on plant growth was much more noticeable in Sunmate; it had the least effect on shoot height and the biggest effect on shoot weight. Shoot height was reduced under salinity stress by 29%, 12%, and 13% in the Sunmate, Suntop and CM72 cultivars, respectively. The fresh shoot weight was reduced by 68%, 55%, and 59%, while shoot dry weight was reduced 68%, 53%, and 49% in Sunmate, Suntop and CM72 in salinity stress plants respectively. Similarly, compared to the control plants, the root length was reduced under salinity stress by 37%, 8%, and 24% in Sunmate, Suntop, and CM72, respectively, while the reductions in fresh root weight were 42%, 33%, and 11% and dry root weight were 65%, 39%, and 30% in Sunmate, Suntop, and CM72, respectively in salinity treated plants (Table S1).

**Figure 1.** Morphology (**A**) and growth parameters (**B**) of seedlings of two wheat cv., Suntop and Sunmate, and one barley cv., CM72, 25 days after treatment with NaCl. Control and NaCl represent 0 and 100 mM NaCl, respectively. Values are means + SE (*n* = 4). For each parameter, means annotated with the same letter are not significantly different from each other according to Duncan's Multiple Range Tests at *p* ≤ 0.05.

#### *3.2. Chlorophyll and Photosynthetic Parameters*

Gas exchange parameters were recorded 25 DAT (days after treatment) and significant (*p* < 0.05) decreases in net photosynthetic rate (Pn), stomatal conductance (Gs), intercellular CO2 concentration (Ci), and transpiration rate (Tr) were detected in both wheat and barley in comparison to their respective controls (Figure 2A–F). No significant changes were observed in Gs, chlorophyll contents and *Fv* to *Fm* ratios, however, significant differences were observed in Pn, Ci, and Tr among all the cultivars in the salinity treated plants. Interestingly, the two-salt tolerant cultivars, Suntop and CM72, showed no significant difference in regard to the photosynthetic parameters, but differences were noted in Sunmate, which is salt-sensitive.

**Figure 2.** Effect of salinity stress on photosynthetic traits in two wheat cv., Suntop and Sunmate, and the barley cv., CM72, 25 days after treatment with 100 mM NaCl. Pn (**A**), Gs (**B**), Ci (**C**), Tr (**D**) and *Fv*/*Fm* (**F**), represent net photosynthetic rate, stomatal conductance, intercellular CO2 concentration, transpiration rate, and a maximum quantum yield of photosystem II photochemistry of the second fully expanded leaves, respectively. The chlorophyll content was measured as SPAD (**E**) (Soil Plant analysis Development). Values are means + SE (*n* = 4). For each parameter, means annotated with the same letter are not significantly different from each other according to Duncan's Multiple Range Tests at *p* ≤ 0.05.

#### *3.3. Shoots and Roots Na*+*, K*<sup>+</sup> *Concentration, and Na*+*:K*<sup>+</sup> *Ratio*

Internal Na<sup>+</sup> and K<sup>+</sup> concentrations were determined, salinity significantly (*p* < 0.05) increased Na+, decreased K<sup>+</sup> content, and increased Na+:K<sup>+</sup> ratio in both shoots and roots of all cultivars in the saline-treated plants relative to the control plants (Figure 3A–F). In general, roots contained more Na<sup>+</sup> and K<sup>+</sup> compared to shoots, regardless of the cultivar or treatment. With regard to the plants given the salt treatment, in both shoots and roots, the increase in the Na<sup>+</sup> content followed the trend CM72 < Suntop < Sunmate in both organs, while the K<sup>+</sup> content decreased in the following trend CM72 > Suntop > Sunmate. Interestingly, the increase in Na<sup>+</sup> content among the cultivars was inversely proportional to the decrease in K<sup>+</sup> content. As a consequence of the changes in both minerals, the Na+:K<sup>+</sup> ratio increased under salt stress. The greatest Na+:K<sup>+</sup> ratios were observed in Sunmate (0.339) and Suntop (0.127), while the smallest were observed in CM72 (0.075) in shoots. The same trend in Na+:K<sup>+</sup> values were also observed in the roots (Table S2).

**Figure 3.** Effect of salinity stress on Na<sup>+</sup> and K<sup>+</sup> concentrations and Na+:K<sup>+</sup> ratio in shoots (**A**,**C**,**E**) and roots (**B**,**D**,**F**) of two wheat cv., Suntop and Sunmate, and the barley cv., CM72, 25 days after treatment with 100 mM NaCl. Error bars represent SE (*n* = 4). Different letters indicate significant differences (*p* ≤ 0.05) among the 3 cultivars.

#### *3.4. Lipid Peroxidation Assay and Antioxidative Enzyme Activities*

Lipid peroxidation measurements at 1, 5, 10, and 15 DAT showed that salt stress induced significant changes among the cultivars and treatments (Figure 4A,B and Supplementary Tables S3 and S4). Regardless of the cultivar, MDA contents were significantly increased by the salt treatment in both leaves and roots, indicating enhanced lipid peroxidation. In general, in plants given salinity treatments, the MDA content was highest in Sunmate followed by Suntop then CM72, with the highest increase observed 15 and 10 d after treatment in leaves and roots, respectively.

Significant differences (*p* < 0.05) were detected among the cultivars and between the treatment in both roots and leaves for all measured antioxidant enzymes (Supplementary Tables S3 and S4, and expression relative to control activities in Figure 5A–J). In general, in the leaves, the relative activities of all enzymes were highest in Suntop, followed by CM72 and then Sunmate; however, in the roots, there was little difference between CM72 and Suntop. For SOD in the leaves, the activities of this enzyme were similar on Days 1 and 5, rose on Day 10, and then dropped below those measured on Day 1. In the roots, SOD activities rose on Day 10 and remained high on Day 15. For POD in the leaves, activities rose on Day 5 and then declined during the remainder of the assessment period, whilst in the roots POD activities did not rise until Day 10 and then declined. For APX in the leaves, activities started to rise on Day 10 and were highest on Day 15 whilst in the root's activities rose on Day 10 and remained high. For CAT in leaves, activities were highest on Day 5 and then declined, except in Suntop, where a decline was observed on Day 10. In the roots, CAT activities tended to stay on similar levels throughout the treatment period. Finally, for GR for the two wheat cv., there was little change in activities in the leaves during the assessment period; however, for CM72, GR activities were highest on Days 5 and 10 and then declined. In the roots, for Suntop, GR activities increased on Day 5 and then remained high, for Sunmate, activities increased on Day 10, and then reduced and for CM72, GR activities were high throughout the experiment. Peaks of antioxidant enzyme activity were observed generally earlier in shoots than in roots, while the earliest peak was observed for CAT and the latest for APX.

**Figure 4.** Effect of salinity stress on malondialdehyde contents (MDA, nmol-1 FW) in leaves (**A**) and roots (**B**) of Suntop (Aqua), Sunmate (Red), and CM72 (Orange) 1, 5, 10, and 15 days after treatment with 100 mM NaCl. The data are expressed as a percentage of control content. Different letters indicate significant differences (*p* ≤ 0.05) among the 3 cultivars within each sampling day. Error bars represent SE (*n* = 4).

**Figure 5.** Effect of salinity stress on (**A**,**B**) superoxide dismutase (SOD) (U g−<sup>1</sup> FW); (**C**,**D**) peroxidase (POD) (OD 470 g−<sup>1</sup> min<sup>−</sup>1); (**E**,**F**) ascorbate peroxidases (APX) (mmol g−<sup>1</sup> FW min<sup>−</sup>1); (**G**,**H**) catalase (CAT) (mmol-1 FW min-1; (**I**,**J**) glutathione reductase (GR) (mmol g-1 FW min<sup>−</sup>1) activities in leaves (**A**,**C**,**E**,**G**,**I**) and roots (**B**,**D**,**F**,**H**,**J**) of Suntop (Aqua), Sunmate (Red), and CM72 (Orange) 1, 5, 10, and 15 days after treatment with 100 mM NaCl. The data are expressed as a percentage of control activities. Different letters indicate significant differences (*p* ≤ 0.05) among the 3 cultivars within each sampling day. Error bars represent SE (*n* = 4).

#### *3.5. Leaf and Root Ultrastructure*

The chloroplast ultrastructure of Sunmate was more severely affected by salt stress relative to controls and also to Suntop and CM72. Under control conditions, the chloroplasts of Suntop mesophyll cells usually had normal morphology with distinct grana and stroma lamellae, large starch grains with numerous plastoglobuli and well-organized, round mitochondria (Figure 6A); after the salt treatment, there were fewer plastoglobuli, no starch grains were apparent, and the grana and stroma lamellae were diffuse (Figure 6B). In contrast, chloroplasts of Sunmate were severely damaged by salinity stress, i.e., the chloroplast envelope showed disintegration with reduced grana stacks and less distinct thylakoids membranes, swollen oval-shaped mitochondria and larger osmophilic plastoglobuli (Figure 6C,D). As with Suntop, the chloroplasts of CM72 remained relatively normal in response to the salt treatment except for the disappearance of starch grains and thinner lamellae (Figure 6E,F).

When viewed using transmission electron microscopy, the root cells of all cultivars grown without salt treatment (control) had dense cytoplasm and organelles, and organized and large nuclei and nucleoli (Suntop, Figure 7A; Sunmate, Figure 7C; CM72, Figure 7E). Treatment with salt induced a number of ultrastructural changes from mild to severe, with the most visible alteration being the disappearance of nucleoli and vacuoles in Sunmate (Suntop, Figure 7B; Sunmate, Figure 7D; CM72, Figure 7F). Suntop and CM72 had clear nucleoli and larger and several vacuoles in comparison with Sunmate. However, the size of the nucleoli increased in CM72, and to a lesser extent in Suntop, upon exposure to salt.

**Figure 6.** Transmission electron micrographs of chloroplasts of leaves of Suntop (**A**,**B**), Sunmate (**C**,**D**), and CM72 (**E**,**F**) under control (top panel) and 100 mM NaCl (bottom panel). CW, cell wall; G, grana; MTC, mitochondria; PG, plastoglobuli; SG, starch grains.

**Figure 7.** Electron micrographs of roots of Suntop (**A**,**B**), Sunmate (**C**,**D**), and CM72 (**E**,**F**) under control (top panel) and 100 mM NaCl (bottom panel). CW, cell wall; Nu, nucleolus; N, nucleus; Vac, vacuole.

#### **4. Discussion**

The effects of the treatments differed for different plant organs; therefore, the effects on shoots and roots of both species are considered separately.

Reduced biomass, a marked perturbation in photosynthetic parameters along with reduced chlorophyll contents resulting from salinity stress were observed in the wheat and barley cultivars. These effects are possibly due to either single or combined effects of reduced stomatal conductance, inhibition of metabolic phenomena, and increased ROS generation which can increase oxygen-induced cellular damage [42]. The reductions in stomatal conductance (Gs), photosynthesis rates (Pn), and leaf chlorophyll contents due to salinity were greater in Sunmate than in Suntop and CM72 (Figure 2). In a study conducted using sorghum, Netondo et al. [43] found that changes in stomatal conductance (Gs) and intercellular CO2 concentration (Ci) were positively correlated under salt stress, concluding that stomatal conductance (Gs) was the key factor arresting net photosynthesis rates (Pn) under saline stress. The lower stomatal conductance (Gs) accompanied by low chlorophyll contents in Sunmate could contribute to the higher inhibition of net photosynthesis rates (Pn). Usually, plants close their stomata upon the onset of stressful conditions to save water, consequently reducing stomatal conductance (Gs) and photosynthesis [44]. The effect of salinity might be a secondary influence, arbitrated by the lower partial pressure of CO2 in the green parts of the plant due to the stomata closure on the photosynthesis-related enzyme activities [45,46].

The *Fv*/*Fm* ratio reflects the photochemical efficiency of PSII [47]. Results in this study show that even a small but significant reduction in *Fv*/*Fm* with the greatest decline in Sunmate followed by CM72 then Suntop was recorded (Figure 2); these results are consistent with that presented by Ahmad et al. [48] and Ibrahim et al. [47]. NaCl stress can disturb the photosynthesis biochemistry, limiting the efficiency of two photosystems due to the disordering of chloroplast integrity [47]. In our study, salinity altered leaf chloroplast ultrastructure causing swelling of thylakoids, diffuse granular and stroma lamellae, a larger number of large-sized plastoglobuli and a reduction in leaf chlorophyll content in the sensitive wheat cultivar, Sunmate; these changes were not seen to the same extent in Suntop and CM72 (Figure 6). There may be several reasons for the disruptions to thylakoids and the chloroplast envelope. These include higher accumulation of lipids in chloroplasts, ion toxicity, or

imbalance [49], and osmotic imbalance between chloroplast and stroma [50], which in turn cause a reduction of photosynthetic efficiency and the electron transport activity of chloroplasts [51].

Additionally, upon exposure to salinization, severe disruption of nuclei and nuclear membranes of roots were detected in Sunmate but to a lesser extent in Suntop and CM72 (Figure 7). Salinity largely affects roots because of their direct contact with the soil. Therefore, to protect the whole plant from the adverse effects of salinity, roots should better tolerate salinity stress [47]. Zhang and Blumwald [52] noted that in tolerant plants, Na<sup>+</sup> is kept away from the cytosol by compartmentalizing it into the vacuole, and due to the lack of this ability in sensitive plants, dehydration and ionic imbalance disturbed the metabolic process in sensitive plants [53]. We observed that the nucleolus disappeared in Sunmate. A common consequence of this type of alteration inside the nucleus would be a loss of function and/or even cell death [53].

It is important to determine Na<sup>+</sup> and K<sup>+</sup> concentrations and Na+:K<sup>+</sup> ratios in shoots and roots to understand mechanisms of salinity tolerance [54]. In our study, under the salinity stress, significant differences were found in shoot Na<sup>+</sup> and K<sup>+</sup> along with Na+:K<sup>+</sup> ratios in both species relative to controls, with the most severe effect in Sunmate (Figure 3A–F, Table S2). In general, CM72 accumulated less Na<sup>+</sup> and more K<sup>+</sup> in shoots followed by Suntop and then Sunmate. Hence, the low Na+:K<sup>+</sup> ratios in CM72 and Suntop may explain the tolerance of these cultivars. Root to shoot Na<sup>+</sup> and K<sup>+</sup> translocation is limited, as all genotypes accumulated more Na<sup>+</sup> and K<sup>+</sup> in roots than in shoots. These results are consistent with the idea of differences in translocation restricting Na<sup>+</sup> movement to the shoot being one of the mechanisms of salinity tolerance. Na<sup>+</sup> and K<sup>+</sup> are interdependent under salinity stress. Previous studies have found a decrease in K<sup>+</sup> content in several plant species resulting from high salinity [35,55]. Increased Na<sup>+</sup> concentrations in root zones have an antagonistic effect on K<sup>+</sup> uptake. Consequently, a deficiency of K<sup>+</sup> has created stunting growth and reducing yields [56].

The most general consequence of salinization is the accumulation of hazardous substances in plant cells especially ROS such as singlet oxygen (O2), superoxide radicals (O2 −), and hydrogen peroxide (H2O2); these species cause damage to proteins, lipids, and nucleic acids thereby promoting rapid plant death [57]. Malondialdehyde (MDA), a product of polyunsaturated fatty acid peroxidation [58], is commonly considered as a sign of the extent of oxidation damage under stress [27,59]. The hostile influences of NaCl stress on lipid peroxidation have been reported in other plants, for example in *Brassica juncea* [60] and *Vicia faba* [61], and MDA has been widely recognized as a good salinity tolerance marker in plant species [62]. We found significantly lower MDA contents in the shoots and roots of CM72 and Suntop compared to the Sunmate (Figure 4A,B). These data suggest that CM72 and Suntop were better protected than Sunmate against oxidative damage under salinity stress.

After salt treatment, tolerant plants eventually develop an enhanced antioxidant enzyme system to handle the effects of ROS. In our study, significantly increased SOD, POD, APX, CAT, and GR activities were found in roots and leaf tissues of both species in the NaCl treated plants (Figure 5A–J). However, the relative activities of these enzymes were recorded higher in Suntop and CM72 than in cultivar Sunmate in both tissues. SOD provides the first line of defense against ROS and protects plants from severe damage generated by O2 <sup>−</sup> and H2O2 in the presence of metal ions [63]. Many studies have found that salinity positively promotes SOD activity in tolerant cultivars in both roots [64] and leaves [65]. Subsequent reactions are required to convert the H2O2 produced by SOD to H2O because H2O2 is still toxic to plants and reactions involving POD, CAT, and APX are important. Our research corroborates previous studies [47,48] had measured an enhanced activity of SOD, POD, and CAT in plants treated with a high NaCl dose and the activities of these enzymes were again higher in the two tolerant genotypes. Feki et al. [66] and Koca et al. [67] also demonstrated that tolerance to salinity in wheat and sesame genotypes was associated with lower MDA contents and higher activities of antioxidant enzymes. Thus, it is evident from our results and the results of others [10,68] that the higher POD, CAT, and APX activities coordinate with SOD activity to deal with the undesirable effect of O2 − and H2O2 and the activities of these enzymes are strongly correlated with tolerance to salt-induced oxidative stress in wheat and barley.

The activity of GR in the leaves and roots was higher in CM72 compared to Suntop and Sunmate (Figure 5I,J). Other studies working with salt-sensitive and tolerant genotypes suggested that higher GR activities relate to salt tolerance [46,65]. The higher GR activity might be able to elevate NADP<sup>+</sup> concentrations to gain electrons from the photosynthetic electron transport chain thereby reducing the production of ROS [69]. Our results also suggest that the salt-tolerant cultivar may exhibit a more active ascorbate-glutathione cycle.

#### **5. Conclusions**

Although differences were observed between CM72 and Suntop in the growth and biochemical traits assessed and modified by salt stress, the differences are negligible in comparison with the response to the salt stress of sensitive wheat cv. Sunmate. The distinct differences between wheat and barley were lower MDA content, lower Na+/K<sup>+</sup> ratio and a higher level of APX and GR content in the roots of barley cultivar CM72. These results lead us to infer that differences in response to salinity may be just as great within a species as between species. The most obvious mechanisms for salt tolerance in the tolerant barley and wheat cultivars are the increased activities of ROS-scavenging enzymes and a more balanced Na+:K<sup>+</sup> ratio. Our results indicated that Suntop is highly tolerant against salinity, which is quite similar to barley CM72. Novel salt-tolerant related genes may be identified in Suntop for improving the salt tolerance of wheat cultivars, except for commercial application in saline-alkali soils.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/10/1/127/s1, Table S1: Effect of salinity on plant growth and biomass of Suntop, Sunmate and CM72, 25 days after treatment with 100 mM NaCl, Table S2: Shoot and root Na<sup>+</sup> and K<sup>+</sup> concentration and Na+/K<sup>+</sup> ratio of two wheat cv. Suntop and Sunmate, and one barley cv. CM72, 25 days after treatment with 100 mM NaCl, Table S3: Effect of salinity stress on SOD, POD, CAT, APX and GR activities and MDA contents in the shoots of Suntop, Sunmate and CM72, after 1, 5, 10 and 15 days 100 mM NaCl treatment, Table S4: Effect of salinity stress on SOD, POD, CAT, APX, and GR activities and MDA contents in the roots of Suntop, Sunmate, and CM72, after 1, 5, 10, and 15 days 100 mM NaCl treatment.

**Author Contributions:** Conceptualization, F.W. and M.L.; data curation, M.Z. and F.W.; formal analysis, M.Z.; funding acquisition, F.W.; investigation, M.Z. and S.S.; methodology, F.W. and M.Z.; project administration, F.W.; supervision, F.W.; validation, M.Z., S.S., P.H., and F.W.B.; writing—original draft, M.Z.; writing—review and editing, M.Z., F.W., and P.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was supported by the Key Research Foundation of Science and Technology Department of Zhejiang Province of China (2016C02050-9-7).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Transcriptomic Profiling of Pomegranate Provides Insights into Salt Tolerance**

**Cuiyu Liu 1,2, Yujie Zhao 1,2, Xueqing Zhao 1,2, Jinping Wang 1,2, Mengmeng Gu 3,\* and Zhaohe Yuan 1,2,\***


Received: 17 October 2019; Accepted: 24 December 2019; Published: 27 December 2019

**Abstract:** Pomegranate (*Punica granatum* L.) is widely grown in arid and semi-arid soils, with constant soil salinization. To elucidate its molecular responses to salt stress on mRNA levels, we constructed 18 cDNA libraries of pomegranate roots and leaves from 0 (controls), 3, and 6 days after 200 mM NaCl treatment. In total, we obtained 34,047 genes by mapping to genome, and then identified 2255 DEGs (differentially expressed genes), including 1080 up-regulated and 1175 down-regulated genes. We found that the expression pattern of most DEGs were tissue-specific and time-specific. Among root DEGs, genes associated with cell wall organization and transmembrane transport were suppressed, and most of metabolism-related genes were over-represented. In leaves, 41.29% of DEGs were first suppressed and then recovered, including ions/metal ions binding-related genes. Also, ion transport and oxidation-reduction process were restricted. We found many DEGs involved in ABA, Ca2+-related and MAPK signal transduction pathways, such as ABA-receptors, Ca2<sup>+</sup>-sensors, MAPK cascades, TFs, and downstream functional genes coding for HSPs, LEAs, AQPs and PODs. Fifteen genes were selected to confirm the RNA-seq data using qRT-PCR. Our study not only illuminated pomegranate molecular responses to salinity, but also provided references for selecting salt-tolerant genes in pomegranate breeding processes.

**Keywords:** pomegranate; salt stress; transcriptome; tissue-specific; signaling transduction pathways; transcription factors

#### **1. Introduction**

Soil salinization is defined as the excess or deposition of salt ions in land, which may interfere with plant growth. With the aggravation of soil salinization, it has become a considerable threat to healthy and sustainable development of worldwide agriculture [1]. Approximate 20% of the global cultivated lands and 50% of the irrigated lands are affected by salinity [2]. Plants exposed to saline conditions mainly suffer from osmotic stress, ion toxicity, and nutrient deficiency [3,4]. Consequently, all of the significant processes involved in plant growth and development, such as photosynthesis, protein synthesis, energy conversion and ion balance, could be affected by salinity [5]. The effects of salinity depend, not only on species, genotypes, and the age of plants, but also on the duration and intensity of stress [6]. Meanwhile, plants can adapt to saline environment with following strategies: (1) Efficiently controlling the uptake, transport, and compartmentalization of toxic ions;

(2) synthesis of osmoregulation substances and activation of antioxidant enzymes; (3) formation of unique morphological structures, such as succulent leaves, salt glands and bladders [7,8].

Plant salt-tolerance is a quantitative trait controlled by multiple genes [4]. Using transcriptomic analysis, researchers can identify the salt-related genes with differential and temporal expression patterns in plants, and provide a clear picture of transcripts in response to salt stress. These salt-related genes are identified and categorized into two types according to the functions of proteins [4]. The first type is the effector code the functional proteins, which are directly involved in the physiological and biochemical responses to salt stress in plants. These proteins include the superoxide dismutase (SOD) [9], ascorbate peroxidase (APX) [10], high-affinity potassium transporter (HKT) [11], Na+/H<sup>+</sup> antiporter (NHX) [12], aquaporin (AQP) [13], late embryogenesis abundant (LEA) [14], and H+-ATPase (VHA) [15], etc. The second type is regulator that involved in regulating the expressions of genes and the signal transduction pathways, such as transcription factors (TFs) and various kinases [16]. These well-characterized TFs include Apetala2/ethylene response factor (AP2/ERF), dehydration-responsive element-binding protein/C-repeat binding factor (DREB/CBF), WRKY, NAC, basic leucine zipper (bZIP), MYB, and basic helix-loop-helix (bHLH) family members [4]. These genes regulate the expressions of downstream genes via various ways, then may influence the plants salt-tolerance eventually [17]. Under salt stress, many signaling transduction pathways are stimulated in plants, such as Ca2<sup>+</sup>-related, abscisic acid (ABA), and mitogen-activated protein kinase (MAPKs), as well as the crosstalk networks among them, which play crucial roles in responses to salt stress [4,18].

Pomegranate (*Punica granatum* L.) is an emerging commercial fruit tree of the Lythraceae family [19]. In recent years, pomegranate is increasingly popular with extensive usage of its fruits and products, and it is considered a 'super fruit' with high nutritional and medicinal values [20]. The species is widely grown in arid and semiarid regions where the availability and irrigation of saline water are significant issues [21]. Thus, it is important to explore pomegranate potential salt tolerance. Considerable research has been done on pomegranate physiological and biochemical responses to salinity, especially on its growth, ion balance, osmoregulation, and the scavenging of reactive oxygen species (ROS) [22–24]. In this study, we used roots and leaves of a pomegranate cultivar 'Taishanhong' with whole genomic sequence released to perform deep transcriptome sequencing and then demonstrate a global representation of potential candidate genes under salt stress. We aim to unravel the fundamental molecular mechanisms in pomegranate underlying the responses to salt stress.

#### **2. Materials and Methods**

#### *2.1. Plant Growth and Stress Treatments*

Uniform rooted cuttings of pomegranate ('Taishanhong') were obtained from the Pomegranate Repository of Nanjing Forestry University, China. They were grown in plastic pots (2.5 L) containing medium (1:1 by volume of perlite and peat, 1.5 ± 0.1 Kg) in a climate controlled chamber for six months (14 h light 26 ◦C/10 h dark 22 ◦C), and fertigated weekly with <sup>1</sup> <sup>2</sup> Hoagland's solution [25].

At the beginning of experiment, forty-five plants were randomly selected and divided into 3 groups, fifteen plants per group. The roots and leaves of the first group of plants were harvested as control samples (T0R and T0L). The The other two groups of plants were fertigated once with 500 mL of <sup>1</sup> <sup>2</sup> Hoagland's solution containing 200 mM NaCl. A saucer was placed under each container to collect leachate solution during the experiment period. The roots and leaves of the second group of plants on day 3 (T1R and T1L) and the third group of plants on day 6 (T2R and T2L) were harvested later. There was not visible difference between the treated and untreated plant (Figure S1). Samples from five plants in the same group were pooled together as a replicate due to the small amount of biomass. All samples were immediately frozen in liquid nitrogen and then stored at −80 ◦C.

#### *2.2. RNA Preparation, cDNA Library Construction and Sequencing*

Total RNA was extracted from pomegranate roots and leaves using the Total RNA Kit (Tiangen, Beijing, China) according to the manufacturer's instructions. RNA purity, concentration and integrity were checked using the NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA), Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA), respectively.

A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEB Next® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. The RNA samples were concentrated using magnetic oligo (dT) beads and then broken into short fragments using an RNA Fragmentation buffer (NEB Next First Strand Synthesis Reaction Buffer (5X), Ambion, Austin, TX, USA). The first-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second-strand cDNA was synthesized subsequently using DNA Polymerase I and RNase H. After adenylation of 3 ends of DNA fragments, NEB Next Adaptor with hairpin loop structure were ligated to prepare for hybridization. The cDNA fragments with suitable lengths (150~200 bp) were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA) to construct the final cDNA libraries. Then 3 μL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 ◦C for 15 min followed by 5 min at 95 ◦C before PCR. Next, the selected cDNA fragments were enriched via PCR with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. PCR products were purified with AMPure XP system and the cDNA libraries were assessed on the Agilent Bioanalyzer 2100 system. Finally, the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, Inc., San Diego, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq 4000 platform (Illumina, Inc., San Diego, CA, USA) and 150 bp paired-end reads were generated.

#### *2.3. Sequence Assembly*

Low quality reads and reads containing adapter and ploy-N were filtered by NGS QC ToolKit [26] from raw reads. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. The clean data from all 18 libraries were separately mapped to the pomegranate genome assembly (ASM220158v1) using HISAT2 software [27]. StringTie [28] was used to construct and identify both known and novel transcripts from HISAT2 alignment results. StringTie was also used to count the reads numbers mapped to each gene. After that, FPKM (Fragments Per Kilobase of exon per million fragments Mapped) were calculated based on the length of the gene and reads count mapped to this gene [27].

#### *2.4. Identification and Functional Annotation of DEGs*

All genes were compared against various protein databases by BLASTX, including the Nr (NCBI non-redundant protein sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins), and Swiss-Prot (A manually annotated and reviewed protein sequence database) with an E-value cut-off of 10−5. Then genes with the best BLAST hit (the highest score) were chosen along with their protein functional annotations. DESeq and Q-value were employed to evaluate the differential expression genes (DEGs) between controls and treatments. The false discovery rate (FDR) and log2FC (log of fold change) were calculated for all genes, and only transcripts with FDR < 0.05 and|log2(fold change)|≥1 were considered as DEGs. To annotate the DEGs with gene ontology (GO) terms, the Nr BLAST results were imported into the Blast2GO program [29]. To examine the expression

patterns of DEGs, the expression data from roots and leaves (T0, T1 and T2) were normalized to 0, log2(T1/T0), log2(T2/T0), and then clustered by Short Time-series Expression Miner software (STEM) separately, using a FDR correction method and *p*-value ≤ 0.05 as the cutoff [30]. The results of GO annotations in each pattern were enriched and refined using TBtools v0.6652 (Toolbox for Biologists, https://github.com/CJ-Chen/TBtools). We used KOBAS software [31] to test the enrichment of DEGs in KEGG pathways.

#### *2.5. Validation of RNA-Seq by qRT-PCR*

We selected 15 DEGs as experimental validation by quantitative real-time PCR (qRT-PCR), which were performed with three biological and three technical replicates for each cDNA sample. The primers for these genes were designed with NCBI primer-BLAST (Table S1). Reverse transcription was conducted with the GoScriptTM Reverse Transcription System (Promega, Madison, WI, USA). We conducted qRT-PCR in a 7500 fast Real-Time PCR system (Applied Biosystems, CA, MA, USA), and analyzed results with the ΔΔCt method, and *Pgactin* (F: ATCCTCCGTCTTGACCTTG, R: TGTCCGTCAGGCAACTCAT) gene was used as a reference gene. Each reaction was carried out in a final volume of 20 μL, containing 7.5 μL of ddH2O, 10 μL of SYBR Green PCR master mix, 0.5 μL of each gene-specific primer and 2 μL of diluted cDNA. The PCR thermal cycling conditions were as follows: 95 ◦C for 10 min; 40 cycles of 95 ◦C for 5 s, 60 ◦C for 30 s and 72 ◦C for 30 s. Data were collected during the extension step: 95 ◦C for 15 s, 60 ◦C for 1 min, 95 ◦C for 30 s and 60 ◦C for 15 s.

#### **3. Results**

#### *3.1. Overview of Sequencing and Mapping*

To obtain a global overview of the salt-induced changes at whole-transcriptomic scale in pomegranate, we constructed 18 cDNA libraries from the roots and leaves of 0 (controls), 3 and 6 days after 200 mM NaCl treatment. A high correlation (R2 > 0.80) between biological replicates was observed for all treatments (Figure S2), which indicated that the biological replicates were reliable in this study. In totally, sequencing libraries yielded 519.47 million reads with 150 bp for both paired ends. After adapter removal and refining, we obtained 155.84 Gb of clean data, and Q30 percentage were all above 90.08% (Table 1). The ratios of reads mapping to the pomegranate genome were high, with values ranging from 93.61% to 95.33%. Then, the unique-mapped reads were 91.89–93.79% and the mapped sequences in genome exon regions were 62.88–69.32% (Table 1). Finally, we assembled 29,226 transcripts from the 18 cDNA libraries of pomegranate roots and leaves under salt stress. To splice the genes more completely and accurately, the StringTie software was employed to reconstruct the transcripts, and then retrieved 34,047 assemblies, including 5396 novel transcripts. Finally, 26,444 genes (including 2421 new genes) were assembled from annotated genes using 7 different databases (see methods), and can be used as a reference of pomegranate genome annotation (Table S2).


**Table 1.** Summary of RNA-Seq results and the alignments in pomegranate genome.

T09 26,542,339 7,962,701,700 50.51 91.32 94.63 92.88 67.69


**Table 1.** *Cont*.

Reads Number: total Number of paid-end reads in Clean Data; Base Number: total number of bases in Clean Data; Q30%: the ratio of nucleotides with quality value ≥ 30; GC content: The ratio of guanidine and cytosine nucleotides; Mapped ratio: The percentage of Mapped Reads in proportion of Clean reads.

#### *3.2. Identification and Annotation of DEGs*

A total of 2255 DEGs were identified from salt treatments with FDR < 0.05 and|log2(fold change)|≥1, including 1080 up-regulated and 1175 down-regulated genes, and more DEGs were found in roots (1623 genes) than in leaves (632 genes) (Table S2). To provide insights into the underlying functions of pomegranate transcripts under salt stress, DEGs were annotated after they were compared to well-studied sequences in the GO, Nr, Swiss-Prot, KEGG, and KOG databases (Table S2). In this study, most of DEGs in T1R were distributed in T subgroup (Signal transduction mechanisms), E subgroup (Amino acid transport and metabolism) and G subgroup (Carbohydrate transport and metabolism) (Figure 1a). DEGs of T2R and T1L samples were mostly concentrated in R subgroup (General function prediction only), Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms), and K subgroup (Transcription) (Figure 1b,c). The DEGs in T2L were distributed in O subgroup (Posttranslational modification, protein turnover and chaperones), P (Inorganic ion transport and metabolisms) and Q subgroup (Secondary metabolites biosynthesis, transport and cabalisms) (Figure 1d). These subgroups were close to the metabolism, such as DNA transcription, cell division, protein modification and energy conversion.

The DEGs of roots and leaves were also analyzed by KEGG pathways enrichment. These pathways were mainly classified into 5 categories (Figure S3) and showed the top 20 pathways (Figure 2). The majority of pathways in pomegranate roots and leaves were associated with metabolisms, such as global and overview maps, amino acid metabolism, carbohydrate metabolism, lipid metabolism, and energy metabolism, etc. (Figure S3). The metabolism-related genes were significantly enriched in metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), cysteine and methionine metabolism (ko00270), phenylpropanoid biosynthesis (ko00940), starch and sucrose metabolism (ko00500), amino sugar and nucleotide sugar metabolism (ko005200) and protein processing in endoplasmic reticulum (ko04141) pathways. Also, there were some DEGs involved in environmental adaptation (plant-pathogen interaction, ko04626) and signal transduction (plant hormone signal transduction, ko04075) pathways (Figure 2; Figure S3). Genes clustered in genetic information processing were enriched in folding, sorting and degradation, replication and repair, transcription and translation, including protein processing in endoplasmic reticulum (ko04141), DNA replication (ko03030), RNA polymerase (ko03020), RNA transport (ko03013) and RNA degradation (ko03018). In addition, genes involved in the photosynthetic pathways, such as carotenoid biosynthesis (ko00906), porphyrin and chlorophyll metabolism (ko00860) and photosynthesis-antenna proteins (ko00196) were suppressed by salinity.

**Figure 1.** COG classification of differentially expressed genes (DEGs) in pomegranate roots and leaves. DEGs enrichments in T1R (**a**), T2R (**b**), T1L (**c**) and T2L (**d**).

**Figure 2.** The top 20 KEGG pathway of DEGs in pomegranate root and leaf. The horizontal axis showed an enrichment factor, the smaller the enrichment factor is the more significant enrichment of DEGs in this pathway. While the vertical axis illustrated the KEGG pathway, and the color illustrated the -log(Q-value), red color is more reliable and significance enrichment in this pathway. The size of black spots showed the number of DEGs enriched into each pathway. DEGs enrichments in T1R (**a**), T2R (**b**), T1L (**c**), and T2L (**d**).

#### *3.3. The Expressional Patterns of DEGs*

Numbers of DEGs increased in roots with the duration of salt-stress, which were opposite in leaves. The genes were both with little overlap in roots and leaves at two points of stress time (Figure 3a). Only a small portion (2.0%) of DEGs (20 up-regulated and 24 down-regulated) shared the common expression tendency between roots and leaves, and even less (1.2%) DEGs showed an utterly opposite tendency between two tissues (3) genes were up-regulated in leaves and down-regulated in roots, 23 genes were down-regulated in leaves but up-regulated in roots) (Figure 3b). The remaining majority of DEGs were exclusively up-regulated or down-regulated in either tissue. Almost 85.9% of up-regulated and 53.4% of down-regulated DEGs were activated in roots after 6-days salt stress, most of down-regulated genes in T1L recovered in T2L (Figure 3c,d), and only 22 (2.0%) genes were suppressed in leaves at two points of time.

**Figure 3.** The numbers of up-regulated and down-regulated genes in treatments (T1R, T2R, T1L, and T2L), comparing to the controls (T0R and T0L). (**a**) All DEGs in T1R, T2R, T1L, and T2L; (**b**) Up- or down-regulated genes in leaves and roots; (**c**) Up-regulated genes in T1R, T2R, T1L, and T2L; and (**d**) Down-regulated genes in T1R, T2R, T1L, and T2L.

The STEM was employed to analyze the expressional patterns of DEGs. All the 1,623 DEGs in roots were clustered into 8 profiles, whereby 1,340 DEGs were significantly clustered into 4 profiles with *<sup>p</sup>*-value <sup>≤</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> (Figure 4a). The two up-regulated patterns were Profile 4 (31.55%, 512 DEGs) and Profile 7 (15.16%, 246 DEGs), and two down-regulated patterns were Profile 3 (28.59%, 464 DEGs) and Profiles 0 (7.27%, 118 DEGs). All the 632 DEGs in leaf samples were also clustered into 8 profiles and then 2 enriched profiles, including one first down-regulated and then the up-regulated pattern of Profile 2 (41.29%, 261 DEGs, *<sup>p</sup>*-value <sup>=</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>32), and one down-regulated pattern of Profile 1 (16.61%, 105 DEGs, *<sup>p</sup>*-value <sup>=</sup> <sup>6</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>17) (Figure 4b).

**Figure 4.** The expression patterns of DEGs in pomegranate roots and leaves under NaCl stress using STEM. The P-value in roots; (**a**) and leaves (**b**) are shown. The trends of each pattern showed in black lines and the clustered genes are shown in red lines. GO enrichments of Profile 4 (**c**), Profile 7 (**d**) Profile 3 (**e**) and Profile 0 (**f**) in root libraries (T0R, T1R and T2R); and Profile 2 (**g**) and Profile 1 (**h**) in leaf libraries (T0L, T1Land T2L).

Next, DEGs within these profiles were subjected to GO-term enrichment analysis. They were classified into three main categories, including cellular component, biological process, and molecular function, and the top 30 sub-categories (if more than 30, *p*-value < 0.05) of each category were shown (Figure 4). In roots, genes in up-regulated pattern of Profile 4 were just inducted in T2R. They were enriched in chitinase activity, chitin binding, oxidoreductase activity, metal ion, and cation binding under molecular function; chitin, glucosamine-containing compound, aminoglycan, amino sugar and cell wall macromolecule metabolic process under biological process (Figure 4c; Table S3). Gens were continuously suppressed with the extension of stress in Profile 7, and they were enriched in endopeptidase activity, peptidase activity, hydrolase activity, and catalytic activity under molecular function and proteolysis and metabolic process under biological process (Figure 4d; Table S3). The top subcategories of down-regulated pattern (Profile 3) were extracellular region (GO:0005576), cell wall (GO:0005618) and external encapsulating structure (GO:0030312) under cellular component, and lipid biosynthetic process (GO:0008610), cell wall organization (GO:0071555), and external encapsulating structure organization (GO:0045229) under biological process. (Figure 4e; Table S4). Many DEGs coding cell wall components were mostly suppressed in T2R samples, including four extensins, two pectin acetylesterase, three polygalacturonase, and three xyloglucan endotransglucosylase/hydrolase proteins. At the same time, transmembrane transport (GO:0055085) and oxidation-reduction process (GO:0055114) in Profile 0 were inhibited under salt stress. Many genes were identified in these biological process, such as adenine/guanine permease, aquaporin, cationic amino acid transporter, sugar carrier, and zinc transporter (Figure 4f; Table S4). Profile 2 showed a first down- and then up-regulated pattern in leaves, there were many subcategories, such as proteolysis and oxidation-reduction process under biological process, oxidoreductase activity, metal ion, and cation binding under molecular function. The top subcategories in down-regulated pattern of Profile 1 were coenzyme binding, cofactor binding and oxidoreductase activity under molecular function, and oxidation-reduction process, cation and ion transport under biological process. We identified genes coding ion transporters, such as cation/calcium exchanger, zinc transporter, ammonium transporter 1 member, copper-transporting ATPase, potassium transporter (Figure 4h; Table S5).

#### *3.4. ABA Signaling Pathway*

The ABA signaling pathway mainly consists of three protein classes: ABA receptors (PYR/PYL/RCAR), type 2C protein phosphatase (PP2Cs) and sucrose non-fermenting1-related protein kinase 2 (SnRK2s). In our study, three ABA receptors, *PYLs*, were significantly down-regulated (Figure 5), while 13 *PP2Cs* have different expressional patterns in pomegranate roots and leaves under salt stress. Among these genes, seven *PP2Cs* were up-regulated in roots or leaves, five *PP2Cs* down-regulated, and another gene up-regulated in roots and down-regulated in leaves (Figure 5). In contrast, just one *SnRK2* was identified among the DEGs, and it was up-regulated in leaves, which allowed the accumulation of phosphorylated downstream ABFs and activation of ABA-response genes.

**Figure 5.** DEGs involved in ABA signal transduction pathway in pomegranate under salt stress. (**a**) ABA signal transduction pathway; (**b**) DEGs annotations and expressional levels. The color represents with log2(fold change), the DEGs with log2FC ≥ 1 were up-regulated (red) and log2FC ≤ 1 were down-regulated (green) (FDR < 0.05).

#### *3.5. Ca2*+*-Related Signaling Pathways*

Thirty DEGs involved in Ca2<sup>+</sup>-related signaling pathway were identified, these transcripts coded function proteins included three Ca2<sup>+</sup>-ATPases (ACAs), four cation/H<sup>+</sup> antiporters (CAXs), two cation/calcium exchangers (CCXs), four glutamate receptor (GLRs), ten calcium-binding proteins (CaM/CMLs), two CBL-interacting protein kinases (CIPKs) and five calcium-dependent protein kinases (CDPKs) (Figure 6). Expectedly, two *ACAs* were up-regulated in roots (T2R), and one was down-regulated in leaves (Figure 6). Ten *CAXs* expressed differently, including two *CAXs* located on vacuoles were both up-regulated in T2R. Two *CIPKs* were significantly up-regulated in leaves, but no obvious change was observed in roots. The majority of *CaM*/*CMLs* were up-regulated in roots and down-regulated in leaves. The *CDPKs* were up- or down-regulated in roots and/or leaves. Also, two DEGs, that were involved in the MAPK signaling pathway were identified, and both were down-regulated in leaves (Figure 6).

**Figure 6.** The DEGs involved in the Ca2+-related signal transduction pathways. (**a**) Ca2<sup>+</sup>-related signal transduction pathway; (**b**) DEGs annotations and expressional levels. The color represents with log2(fold change), the DEGs with log2FC ≥ 1 were up-regulated (red) and log2FC ≤ 1 were down-regulated (green) (FDR < 0.05).

#### *3.6. The Transcription Factors (TFs)*

For identification of the TFs in pomegranate transcriptome, all the mapped genes were analyzed by BLAST against the Plant Transcription Factor Database (PlantTFDB, http://planttfdb.cbi.pku.edu.cn). A total of 1346 TFs that were classified into 55 putative families (Table S6). Among these TFs, twenty-seven TF families, including 151 genes expressed differently under salt stress compared to controls. The most abundant of differential expression TFs included NAC, ERF, MYB-related, C2H2, MYB, bHLH, GRAS, WRKY, LBD, B3 and bZIP family genes (Table 2; Table S6).

Under NaCl stress, we found that 12 of 19 *NAC* genes were up-regulated, and 10 of 19 genes were down-regulated in leaves or roots of pomegranate plants. Interestingly, we also observed that 18 *PgNACs* significantly changed in T2R when compared to controls, and 4 genes were up-regulated or down-regulated in roots, but reversed in leaves (Table 2; Table S6). These results suggested that many *NACs* were involved in salt stress, but there were different potential responding mechanisms of NAC domain genes to salinity. Thirteen *MYB* genes were induced by salt treatment in pomegranate (Table 2). Among these genes, six in roots and four in leaves were up-regulated, and four in roots and eight in leaves were down-regulated. Moreover, sixteen *MYB-related* genes were detected in our RNA-Seq analyses, thirteen genes were down-regulated, and seven genes were up-regulated by salt stress (Table 2; Table S6). But so far, little was reported that the MYB-related type proteins are related with the responses to salt stress.

The AP2/ERF superfamily is divided into three families: the AP2 family proteins containing two repeated AP2/ERF domains, the ERF family proteins, containing a single AP2/ERF domain, and the RAV family proteins containing a B3 domain. There were 19 *PgERFs*, and 3 *PgAP2s* significantly expressed in treatments when compared to controls (Table 2). Fifteen *ERFs* were down-regulated, and all of *AP2* were repressed by salinity. There were 9 up-regulated and 8 down-regulated *C2H2*, three up-regulated and nine down-regulated *bHLH*, five up-regulated and three down-regulated *WRKY*

identified in pomegranate roots and leaves. Remarkably, nine *GRASs* were up-regulated in the root, but the no-significant change in leaves (Table 2; Table S6).


**Table 2.** The different expressional TFs in pomegranate under NaCl stress.

#### *3.7. qRT-PCR Validation*

To confirm the reliability of the expression levels obtained from the RNA-Seq, fifteen DEGs, including five *PgERFs*, three *PgMYBs* and seven *PgNACs* were selected for qRT-PCR assays. The results showed that the expression level of each transcript closely corresponded to the transcript level estimated from the sequence data with R<sup>2</sup> <sup>≥</sup> 0.85 (Figure 7), which implies reproducibility and accuracy of the RNA-Seq results.

**Figure 7.** The relationship between qRT-PCR and RNA-Seq date. The color circle is gene expression level. The lines represent co-relationship with dependent linear equations.

#### **4. Discussion**

RNA-sequencing techniques have proven to be beneficial and economical for scanning the transcribed genes, both in model and non-model plant species. In this study, we used the RNA-Seq approach as a powerful tool to elucidate the molecular responses to salt stress in pomegranate. Finally, we reconstructed the transcripts and identified 5396 novel genes from the 18 cDNA libraries. There were significant differences between the percentage of total mapped reads and the mapped reads in exonic regions. Theoretically, the ratio of sequencing reads that produced from mature mRNA mapped into exonic regions could be 100%. However, a portion of 30.68–37.12% reads were also mapped into intronic and intergenic regions for the following reasons: (1) The coding and non-coding sequences in the reference genome might have been annotated incorrectly; (2) the new alternative splicing (SNPs) and novel genes resulted from the alternative mRNA splicing and transcripts reconstructing; (3) variations exist between sequencing and reference genome genotypes, which are attributed to the gaps and diversity of their sequences [32,33].

Under salt stress, plants firstly experience osmotic stress due to a disorder of water uptake. Then a process of gradual recovery due to the partially or completely re-established uptake of water occurs within days after the salt treatment [34,35]. We identified 2255 DEGs through the 18 mRNA libraries. From the third to sixth day after treatment, there were significant differences between pomegranate roots and leaves, with little overlap of DEGs responding salinity (Figure 3). More salt-related genes were up-regulated in roots, while the majority of suppressed genes recovered with salt-treating process in leaves. In contrast, the osmotic adjustment of leaves started only after roots had reached new water equilibrium [36]. These tissue-specific and time-specific differences between roots and leaves have also been reported in *Arabidopsis thaliana* (L.) [37], *Populus euphratica* (Oliv.) [36], and *Millettia pinnata* (L.) [38]. Many DEGs were involved in down-regulated patterns in roots, such as cell wall organization, transmembrane transport and oxidation-reduction process, but the proteolysis and metabolic process were over-expressed. Otherwise, most of DEGs involved in oxidation-reduction process and ion transport in leaves were suppressed (Figure 4). The cell wall plays an important role in protecting plant from salt toxicity [39]. In our study, extensin, pectin acetylesterase, polygalacturonase, and xyloglucan endotransglucosylase/hydrolase proteins were identified in roots (Table S4). These proteins are important components of the cell wall, which affect plant growth by mediating cell enlargement and expansion [39]. Many DEGs involved in transmembrane transport code for various channels, carriers and pumps, and ion transporters, such as aquaporin, cationic amino acid transporter, sugar carrier, and zinc transporter (Table S4). They worked together to transport coenzymes, amino acids, carbohydrates, and ions under salt stress. The suppressed genes coding proteins, such as dehydrogenase, cytochrome P450s and peroxidase involved in oxidation-reduction process, play essential roles in responding to salinity [40,41] (Tables S4 and S5). These results indicated that salinity affected pomegranate growth and development by inhibiting the cell division and transmembrane transport, as well as slowing down the redox reactions.

On the other hand, over-represented genes in roots were enriched in metabolic process, carbohydrate metabolic process, and catabolic process, etc. These metabolism-related genes in roots were also enriched in KEGG pathways, such as metabolic pathways (ko01100), biosynthesis of secondary metabolites (ko01110), phenylpropanoid biosynthesis (ko00940), starch and sucrose metabolism (ko00500), amino sugar and nucleotide sugar metabolism (ko005200) (Figure 2; Table S3). The carbohydrates such as glycan, starch, sucrose, amino sugar have been demonstrated as osmolytes and energy resouces in plant responses to salt stress, especially in halophyte species [42]. The accelerated metabolisms in pomegranate might contribute to its adaption responses to salt stress. Similar expression patterns of metabolism-related genes were observed in *Oryza sativa* (L.) [43], *Thellungiella halophila* (C. A. Mey.) [42] and *Helianthus tuberosus* (L.) [44] under salt stress. Metal ions/cations, such as K+, Ca2+, Mg2+, Fe2+, Cu2<sup>+</sup>, and Zn2<sup>+</sup>, act as cofactors participating in the catalytic activity of the enzymes they bind to [44]. The up-regulated DEGs involved in metal ion/cation binding might accelerate the catalytic activity. At the same time, the ion transport process was inhibited in roots and leaves (Tables

S3 and S5). The restriction of ion transport indicated that the excess of Na<sup>+</sup> inhibited the uptake of mineral ions [4]. The results correspond to the physiological responses in other pomegranate cultivars under salt stress [45]. We suspected that pomegranate could cope with the toxic ions to some extent via decreasing uptake of ions from soil and increasing utilization of ions in cells. The toxic ions access into the leaf cells with the transpiration, but the detrimental effect is a slow process taking weeks or months, eventually cause salt toxicity in leaves [46]. Plants were only exposed to salinity stress just for six days in our research, so further study is needed to reveal the long-term effect of salinity.

ABA is a critical hormone, which regulates plant growth, development, as well as responses to environmental stresses such as salinity, drought, heat, cold, wound, and pathogen [47,48]. Once ABA is produced, ABA-bound receptors bind ABA, inhibit *PP2Cs*, such as *ABI1* and *ABI2*, thereby activate *SnRK2s* [49]. The *SnRK2s* are involved in plant responses to abiotic stresses, which can phosphorylate the ABA-responsive element binding factor (ABF) and trigger the expression of ABA-responsive genes [50]. In our study, three ABA receptors, *PYLs*, were significantly down-regulated in roots or leaves (Figure 5), which is consistent with results in tea (*Camellia sinensis* L.) [51] and grape (*Vitis vinifera* L.) [52]. The down-regulation of *PYLs* can reduce plant sensitivity to ABA in response to salt stress and help plants adapt to salinity. The expressional differences of *PP2Cs* between roots and leaves in pomegranate are similar to other plants, which indicate that *PP2Cs* may participate in response to salinity via different strategies [52,53]. *PYR*/*PYLs* may regulate *SnRK2s* directly or indirectly, however, whether *SnRK2s* responding to various abiotic stresses in an ABA-independent or ABA-dependent pathway needs further investigation.

Ca2<sup>+</sup> referred as a second messenger plays a very important role in many stress-related signaling transduction pathways. Under various abiotic stresses, over-accumulation of ions induce a temporary fluctuations in plant cytosolic ([Ca2+]cyt) levels [54]. The genes involved in CIPKs and CDPKs pathways, such as *ACAs*, *CAXs*, *GLRs*, *CaM*/*CMLs*, *CBLs* etc., participate in transporting and binding Ca2+, sensing and relaying signals in plant cell under salinity stress [55,56]. In our study, we found two *ACAs* and two *CAXs* located on vacuoles were up-regulated in T2R, and the majority of *CaM*/*CMLs* were up-regulated in roots but down-regulated in leaves (Figure 6). Previous reports revealed a central role for CaMs in the regulation of Ca2<sup>+</sup> channels and pumps, like CNGCs and ACAs [57,58]. The negative effect of CaM on CNGC activity provides a direct feedback pathway to restrict the influx of Ca2<sup>+</sup> into plant cells [58]. The CaM stimulates the activity of these *ACAs* by preventing their auto-inhibition [58]. Collectively, under salt stress, pomegranate plants restrain the excess influx of Ca2<sup>+</sup> into root cells under salinity. Meanwhile, it coped with the excess Ca2<sup>+</sup> via removing Ca2<sup>+</sup> from the cytosol by *ACAs* and *CAXs* in plasma membrane, including efflux of excess Ca2<sup>+</sup> into the outer rhizosphere and/or the influx into vacuoles [59,60].

Specifically, transcription factors, such as NAC, MYB, AP2/ERF, bHLH, WRKY, GRAS family genes, regulated the expression of downstream genes to cope with various stresses [17]. Our study indicated that the TFs of pomegranate participated in the response to salinity with different patterns. For example, most of *NAC* and *C2H2* genes were up-regulated in roots, while most of *MYB* and *MYB-related* genes were down-regulated in leaves (Table 2; Table S6). Interestingly, nine *GRASs* were all up-regulated in T2R, which is mainly due to that these proteins play roles in plant development, including root development, axillary shoot development, and maintenance of the shoot apical meristem [61]. We also found 57 DEGs coding for salt-inducted effectors, including two *SODs*, two *APXs*, nineteen *PODs*, five *LEAs*, five *AQPs*, ten *HSFs*, three *AKTs,* and one *HKT* (Table S7). Two *APXs* were up-regulated and two *SODs* were down-regulated in leaves. Fourteen and nine *PODs* were suppressed in roots and leaves, respectively. Most of *AQPs* (4 of 5) were down-regulated, while all of *LEAs* (5 of 5) and majority of *HSFs* (9 of 10) were up-regulated in pomegranate tissues under salt stress. *LEAs* and *HSPs* were reported to prevent protein denaturation and maintain cell membrane fluidity under stress [62]. Then the DEGs of *LEA* and *HSF* in pomegranate contributed to mitigating salt stress. The expressions of three *AKTs* and one *HKT* were significantly down-regulated, which suggested that the influxes of K<sup>+</sup> and Na<sup>+</sup> in the cell were restricted [63].

Briefly, when pomegranate plants are exposed to high salinity, Na<sup>+</sup> enters the cell and simultaneously generates stress signals. The messengers Ca2+, ROSs (reactive oxygen species), and ABA transduced the signals, and then the cascades, including Ca2<sup>+</sup>-sensors, MAPK cascades and TFs, involved in secondary and tertiary regulatory networks were activated (Figure S4). TFs regulated the expression of downstream functional genes, such as *HSPs*, *LEAs*, *AQPs*, *SODs,* and *APXs*, to cope with salt stress.

#### **5. Conclusions**

In summary, our research sheds light on the pomegranate molecular response mechanisms to salinity. The differentially expressed genes could also provide references for selecting salt-tolerant breeding materials in pomegranate breeding processes.

**Accession Codes:** The raw data were deposited in NCBI Sequence Read Archive (SRA) database (http://www. ncbi.nlm.nih.gov/Traces/sra) under accession number SRP148507.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/10/1/44/s1, Figure S1: The images of untreated and treated plants. Figure S2: Pearson correlation of all samples, Figure S3: Pathways classification into cellular process, environmental information processing, genetic processing, metabolism and organismal system based on the KEGG analysis, Figure S4: Signaling networks in pomegranate under salt stress, Table S1: The primers used for qPCR validation of 15 DEGs, Table S2: Total numbers of DEGs assigned to databases, Table S3: DEGs involved in up-regulated patterns of pomegranate roots, Table S4: DEGs involved in down-regulated patterns of pomegranate roots, Table S5: DEGs involved in expressional patterns of pomegranate leaves, Table S6: The numbers of differently expressed TFs in pomegranate under NaCl stress, Table S7: The salt-inducted DEGs in pomegranate under NaCl stress.

**Author Contributions:** Conceptualization, C.L. and Z.Y.; methodology, C.L. and Y.Z.; formal analysis, C.L., Y.Z., J.W., and X.Z.; investigation, C.L., Y.Z., X.Z., and J.W.; writing—original draft preparation, C.L.; writing—review and editing, C.L., X.Z., M.G., and Z.Y.; supervision, M.G. and Z.Y.; Funding acquisition, Z.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Initiative Project for Talents of Nanjing Forestry University (GXL2014070, GXL2018032), the Natural Science Foundation of Jiangsu Province (BK20180768), the Doctorate Fellowship Foundation of Nanjing Forestry University, and the Priority Academic Program Development of Jiangsu High Education Institutions (PAPD).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Treatment of Sweet Pepper with Stress Tolerance-Inducing Compounds Alleviates Salinity Stress Oxidative Damage by Mediating the Physio-Biochemical Activities and Antioxidant Systems**

### **Khaled A. Abdelaal 1, Lamiaa M. EL-Maghraby 2, Hosam Elansary 3,4, Yaser M. Hafez 1, Eid I. Ibrahim 5, Mostafa El-Banna 6, Mohamed El-Esawi 7,8 and Amr Elkelish 9,\***


Received: 28 November 2019; Accepted: 19 December 2019; Published: 23 December 2019

**Abstract:** Salinity stress occurs due to the accumulation of high levels of salts in soil, which ultimately leads to the impairment of plant growth and crop loss. Stress tolerance-inducing compounds have a remarkable ability to improve growth and minimize the effects of salinity stress without negatively affecting the environment by controlling the physiological and molecular activities in plants. Two pot experiments were carried out in 2017 and 2018 to study the influence of salicylic acid (1 mM), yeast extract (6 g L−1), and proline (10 mM) on the physiological and biochemical parameters of sweet pepper plants under saline conditions (2000 and 4000 ppm). The results showed that salt stress led to decreasing the chlorophyll content, relative water content, and fruit yields, whereas electrolyte leakage, malondialdehyde (MDA), proline concentration, reactive oxygen species (ROS), and the activities of antioxidant enzymes increased in salt-stressed plants. The application of salicylic acid (1 mM), yeast extract (6 g L<sup>−</sup>1), and proline (10 mM) markedly improved the physiological characteristics and fruit yields of salt-stressed plants compared with untreated stressed plants. A significant reduction in electrolyte leakage, MDA, and ROS was also recorded for all treatments. In conclusion, our results reveal the important role of proline, SA, and yeast extracts in enhancing sweet pepper growth and tolerance to salinity stress via modulation of the physiological parameters and antioxidants machinery. Interestingly, proline proved to be the best treatment.

**Keywords:** *Capsicum annuum* L.; salt stress; salicylic acid; yeast; proline

#### **1. Introduction**

Global food safety is seriously dependent on crops and their supplies, which require considerable increases for servicing the gap between production and demand [1]. The necessity of improving crop production has been much more emergent in the last few years due to the expanding population, which will exceed to 9.7 billion by 2050. Undoubtedly, increases in the population will exert pressure on crops and food resources [1]. Simultaneously, global warming, as well as various biotic and abiotic stresses, hinder the growth and yields of agricultural crops [2]. Among abiotic stresses, salinity is recognized as one of the main restricting factors affecting the growth and productivity of agricultural crops, especially in arid and semiarid regions [3]. Salinity stress causes a reduction in growth and biomass, chlorophyll degradation, water status modification, malfunctions in stomatal functions, modifications in transpiration and respiration, and disequilibria in ion ratios [4,5]. Furthermore, plants develop cytotoxic-activated oxygen under saline conditions, which might seriously interfere with healthy metabolisms as a result of the oxidative damage of lipids, proteins, and nucleic acids [6,7]. Salinization may additionally lead to the excessive intracellular generation of reactive oxygen species (ROS) such as hydroxyl radicals (OH) and superoxide radicals (O2 −) [8]. Plants confront these sorts of oxidants by developing several defensive mechanisms, including antioxidant enzymes and molecules that eliminate potentially cytotoxic types of activated oxygen [9,10].

Sweet pepper (*Capsicum annuum* L.) is an important vegetable crop that is grown for local consumption, and which has a high economic value in the Egyptian agricultural market. Farmers started to utilize saline water to partially fulfil crop water demands. The pepper plant is not a salt-tolerant vegetable, and about 14% of fruit yield loss occurs as a result of each increase in salt level of 1.0 dS/m [11]. Previous investigations have been conducted to mitigate the harmful impact of salt stress on sweet pepper, but most have not been sufficient or broadly applicable. As a result, the search for cheaper, ecologically-friendly strategies for salinity amelioration which enhance the growth and productivity of sweet pepper has been very important to the agriculture sector [12].

Numerous studies have found that implementing exogenous chemicals improves salt stress tolerance in plants [13]; examples of such chemicals are phytohormones such as salicylic acid, sterols, and methyl jasmonate [2,14]. Other chemicals such as polyamines, melatonin, and sodium nitroprusside have also been used to enhance the tolerance of various crop plants to saline conditions [15].

Salicylic acid is an essential phenolic compound that regulates plant growth processes and responses to different environmental factors [16]. It is a stress tolerance inducer and an important signal in many physiological processes, such as proline metabolism and photosynthesis. It reduces oxidative stress in plants under environmental stress and enhances plant growth and productivity under salt- [17] and drought-stress conditions [18]. Foliar application of SA-enhanced growth characteristics of sweet pepper plants [6] has increased the chlorophyll concentrations and enzyme activitives in barley plants, as well as counteracting the deleterious impacts of salinity on faba beans [19]. Yeast extracts are the main source of various important compounds, such as amino acids, phytohormones, and vitamins [20,21]. The use of active yeast extracts has been shown to decrease the damaging impact of drought conditions on pea plants, and enhanced the growth performance and yield of stressed plants [22]. Yeast extract applications have led to improvements in the growth characteristics of bean and corn plants, such as the dry weight of leaves, the leaf area, and the number of leaves under drought conditions [23]. The application of yeast and NPK fertilizers has significantly enhanced chlorophyll concentrations and root yields in sugar beet plants [22]. Seaweed extracts have also improved plant tolerance to abiotic stresses. For example, the application of *Ecklonia maxima* seaweed extract has been shown to enhance the tolerance of zucchini squash plants to salinity stress by improving plant performance, shoot biomass yield, fruit quality, leaf gas exchange rate, SPAD index, and leaf nutritional status under saline conditions [24]. Furthermore, proline has a positive impact on the activity of enzymes and osmotic adjustment under stress conditions, while protecting enzyme denaturation and modulating osmoregulation [25]. The application of proline-modulated antioxidant enzymes such as peroxidase (POX) and catalase (CAT) in tobacco plants under salinity conditions plays a significant role in protein

synthesis and accumulation in plants under stress conditions like drought and salinity in order to enhance the growth characteristics and yield [26–30].

Considering the variable effectiveness levels of salicylic acid, proline, and yeast extract on plants, as well as the harmful impact of salinity stress on the growth and productivity of important crops, the present study aims to evaluate and compare the levels of effectiveness of the three stress tolerance inducers, i.e., salicylic acid (1 mM), yeast extract (6 g L−1), and proline (10 mM), on the growth characteristics, antioxidants, physiological and biochemical parameters, and yield of sweet pepper plants (*Capsicum annuum* L.) grown under the same saline conditions in order to determine which stress tolerance inducer should be recommended for further enhancements of crop performance and tolerance.

#### **2. Materials and Methods**

#### *2.1. Experiments Design and Treatments*

Pot experiments were performed at Agricultural Botany Department, Faculty of Agriculture, Kafrelsheikh University, Egypt during the growing seasons of 2017 and 2018. Laboratory analyses were carried out at the Plant Pathology & Biotechnology Lab, and the EPECRS Excellence Center Kafrelsheikh University, Egypt. This research was conducted to study the impacts of salicylic acid (1 mM), yeast extract (6 g L−1), and proline (10 mM) on the growth characteristics and biochemical and yield parameters of salt-stressed sweet pepper plants (*Capsicum annuum* L.). Irrigation water was artificially salinized by applying NaCl at concentrations of 2000 and 4000 ppm. The seeds of sweet pepper cv. California Wonder were obtained from Sun Seed Company in USA. Ten seeds were sown in the nursery using foam trays. Forty-two days after sowing, seedlings were transplanted into pots (30 cm diameter); each pot contained 8 kg soil and 2 plants. The physical and chemical soil characteristics were recorded, according to the methods described by Abdelaal et al. [21], as follows. pH: 8.2; N: 32.4 ppm; P: 10.5 ppm; K: 289 ppm; electrical conductivity: 1.8 dS m−1, soil organic matter: 1.9%; sand: 17.3%; silt: 35.5%; and clay: 47.2%. Fertilizers were added in two equal doses as recommended (NPK, 135:40:35 kg/ha), plus essential micronutrients, whereas the first dose was added 15 days after transplanting and the second at the beginning of flowering stage [31]. The plants were treated twice (20 and 40 days after transplanting) with salicylic acid (1 mM), yeast (6 g L<sup>−</sup>1) and proline (10 mM). The experiment was done in a completely randomized design with five replicates (five pots with two plants each), and the following measurements were recorded after collecting the plant samples.

#### *2.2. Physiological and Biochemical Analysis*

For physiological and biochemical analyses, the samples were collected at 90 days after transplantation for use in the following assays.

#### 2.2.1. Chlorophyll a and b Determination

For chlorophyll a and b determination, 5 mL N-N Dimethyl formamid was added to 1 g sweet pepper fresh leaves and placed in a refrigerator for 24 h. Following the centrifugation at 4000 *g* for 15 min, the optical density was calculated using spectrophotometer at 647 and 664 nm, according to Moran [32].

2.2.2. Calculation of Leaves Relative Water Content (RWC %) and Electrolyte Leakage (EL %)

The relative water content (RWC) in leaves was recorded according to the formula of Sanchez et al. [33] as follows: RWC = (FW − DW)/ (TW − DW) × 100, where FW is fresh weight, DW is dry weight, and TW is turgid weight. Electrolyte leakage (EL %) was estimated using the formula of Dionisio-Sese and Tobita [34] as follows: EL (%) = Initial electrical conductivity/final electrical conductivity × 100.

#### 2.2.3. Proline Content Determination

Proline was assayed according to the method described by Bates et al. [35] with minor modifications. In brief, a plant sample (0.6 g) was extracted in sulfosalicylic acid (5%) followed by centrifugation at 10000 *g* for 7 min. The supernatants were diluted with water, mixed with 2% ninhydrin, heated at 94 ◦C for 30 min, and then cooled. Toluene was then added to the mixture, and the upper aqueous phase was spectrophotometrically assayed at 520 nm.

#### 2.2.4. Calculation of Lipid Peroxidation and Reactive Oxygen Species (Superoxide and Hydrogen Peroxide)

The lipid peroxidation as malondialdehyde (MDA) in plant samples was calculated according to the method described by Heath and Packer [36] with minor modifications. In brief, 0.6 g of plant sample was extracted in TCA (0.1%), followed by centrifugation at 13,000 *g* for 8 min. The supernatants were mixed with thiobarbituric acid (0.5%) and TCA, and heated at 92 ◦C for 35 min, followed by cooling and centrifugation at 12,000 *g* for 8 min. Next, the supernatants' absorbance was measured at 532 and 660 nm. Superoxide and hydrogen peroxide levels were also determined according to the method described by Badiani et al. [37].

#### 2.2.5. Antioxidant Enzymes Activity (CAT and POX)

Plant samples (1.5 g) were extracted in Tris-HCl (100 mM, pH 7.5) containing Dithiothreitol (5 mM), MgCl2 (10 mM), EDTA (1 mM), magnesium acetate (5 mm), PVP-40 (1.6%), aphenylmethanesulfonyl fluoride (1 mM), and aproptinin (1 μg mL<sup>−</sup>1). The mixed solutions were filtered and centrifuged for 8 min at 13,000 rpm. The supernatants were utilized to record enzymes activities. The activity of CAT and POX of leafy samples was determined according to the method described by Aebi [38] and Hammerschmidt et al. [39]. The supernatant absorbance was shown spectrophotometrically to be 470 nm.

#### 2.2.6. Fruit yields

At 120 days after transplanting, the number of fruits per plant, the fruit fresh weight per plant (g), and the total fruit yield (ton hectare<sup>−</sup>1) were recorded.

#### *2.3. Statistical Analysis*

Data represent the mean ± SD (standard deviation). Two-way analysis of variance was performed using SPSS ver. 19 (SPSS Inc., Chicago, IL, USA). A Tukey's test was also carried out to determine whether a significant difference (*p* < 0.05) existed between mean values.

#### **3. Results**

#### *3.1. Chlorophyll a and b Concentrations*

According to our results in Figure 1, the concentrations of chlorophyll a and b were significantly decreased in sweet pepper plants under salt-stress conditions; the lowest values were recorded with 4000 ppm compared with 2000 ppm and control plants in the two growing seasons. However, the salt stressed plants treated with salicylic acid, yeast extract, and proline showed significant increases in chlorophyll a and chlorophyll b concentrations compared with stressed untreated plants in both seasons. Under salt stresses of 2000 and 4000 ppm, the maximum concentrations of chlorophyll a and b were recorded with proline treatment in both seasons.

**Figure 1.** Effect of salinity stress (2000 and 4000 ppm NaCl) and supplementation of SA, yeast, and proline on the contents of (**A**) chlorophyll a, (**B**) chlorophyll b, (**C**) relative water content (RWC) in sweet pepper in the seasons of 2017 and 2018. Data is mean (±SE) of five replicates. Different letters in each Figure represent significant differences at *p* < 0.05.

#### *3.2. Relative Water Content (RWC %)*

Data obtained in Figure 1 showed that RWC decreased considerably in salt stressed plants; the greatest reduction was recorded in the plants exposed to salinity at 4000 ppm compared with control plants. The exogenous application of salicylic acid (1 mM), yeast (6 g L<sup>−</sup>1), and proline (10 mM) caused a significant increase in RWC in salt stressed plants (2000 and 4000 ppm) compared with salt stressed untreated plants. Furthermore, the best treatments under salinity of 2000 ppm were salicylic acid and proline. Under salt treatment at 4000 ppm, the application of yeast extract (6 g L−1) and proline (10 mM) showed the highest RWC in sweet pepper plants compared with SA treatment in stressed untreated plants in both seasons.

#### *3.3. Electrolyte Leakage (EL %)*

It may be noted from Figure 2 that salt stress at 2000 and 4000 ppm caused a significant increase in electrolyte leakage (EL); the maximum increase was recorded with a salinity level of 4000 ppm in both seasons. Interestingly, electrolyte leakage was significantly decreased upon the foliar application of salicylic acid (1 mM), yeast extract (6 g L<sup>−</sup>1), and proline (10 mM) compared with control plants in both seasons. The best treatment was proline under a salt stress of 2000 ppm in both seasons (Figure 2).

**Figure 2.** Effect of salinity stress (2000 and 4000 ppm NaCl) and supplementation of SA, yeast, and proline on the contents of (**A**) electrolyte leakage, (**B**) proline in sweet pepper in the seasons of 2017 and 2018. Data is mean (±SE) of five replicates. Different letters in each Figure represent significant differences at *p* < 0.05.

#### *3.4. Proline Concentration*

It is evident that proline had markedly accumulated in sweet pepper plants; the highest concentration was recorded with a salinity at 4000 ppm in comparison to the control plants (Figure 2). Intriguingly, the application of salicylic acid, yeast extract, and proline resulted in enhanced proline concentration under all salinity levels; the greatest result was observed with proline (10 mM).

#### *3.5. Lipid Peroxidation (MDA) and Reactive Oxygen Species (Superoxide and Hydrogen Peroxide).*

The results showed that lipid peroxidation (i.e., malondialdehyde or MDA), superoxide, and hydrogen peroxide were significantly increased under salt conditions compared with control plants in both seasons (Figure 3). The maximum levels of MDA, superoxide, and hydrogen peroxide were recorded at a salinity level of 4000 ppm, followed by 2000 ppm, in both seasons. On the other hand, the application of salicylic acid, yeast extract, and proline significantly reduced MDA, O2 −, and H2O2 concentrations under all salinity levels compared to the stressed untreated plants. The best results were obtained with SA and proline.

**Figure 3.** Effect of salinity stress (2000 and 4000 ppm NaCl) and supplementation of SA, yeast, and proline on the contents of (**A**) lipid peroxidation, (**B**) superoxide, (**C**) hydrogen peroxide in sweet pepper in the seasons of 2017 and 2018. Data is mean (±SE) of five replicates. Different letters in each Figure represent significant differences at *p* < 0.05.

#### *3.6. Antioxidant Enzymes Activity*

Antioxidant enzyme activities (CAT and POX) were assayed. The data presented in Figure 4 shows that the plants exposed to salt stress (2000 and 4000 ppm) had higher CAT and POX activity compared with control plants in both seasons. On the other hand, applications of salicylic acid (1 mM), yeast extract (6 g L<sup>−</sup>1), and proline (10 mM) led to reductions in the activities of CAT and POX in the salt-stressed plants.

**Figure 4.** Effect of salinity stress (2000 and 4000 ppm NaCl) and supplementation of SA, yeast, and proline on the activity of (**A**) Catalase (CAT), (**B**) peroxides (POX) in sweet pepper in the seasons of 2017 and 2018. Data is mean (±SE) of five replicates. Different letters in each Figure represent significant differences at *p* < 0.05.

#### *3.7. Number of Fruits per Plant, Fruit Fresh Weight, and Total Fruit Yield (Ton Hectare*<sup>−</sup>*1).*

According to our findings in Figure 5, salt stress at 2000 and 4000 ppm caused significant decreases in fruit number per plant, fruit fresh weight, and the total fruit yield (ton hectare−1) in both seasons. The lowest values of these traits were recorded with salt stressed plants at 4000 ppm concentration, followed by 2000 ppm. Nevertheless, the exogenous application of salicylic acid, yeast, and proline significantly improved the number of fruits per plant, fruit fresh weight, and total fruit yield (ton hectare<sup>−</sup>1) in the stressed treated plants compared with stressed untreated plants.

Interestingly, SA and proline treatments gave the maximum values of the three studied characteristics at a salinity concentration of 2000 ppm in the two seasons (Figure 5). Under salinity stress of 4000 ppm, the best results of fruit number per plant, fruit fresh weight, and total fruit yield (ton hectare<sup>−</sup>1) were recorded with proline, followed by SA and yeast extract in both seasons.

**Figure 5.** Effect of salinity stress (2000 and 4000 ppm NaCl) and supplementation of SA, yeast, and proline on (**A**) number of fruits plant<sup>−</sup>1, (**B**) fresh weight plant<sup>−</sup>1, (**C**) total fruit yield (ton hectare<sup>−</sup>1) in sweet pepper in the seasons of 2017 and 2018. Data is mean (±SE) of five replicates. Different letters in each Figure represent significant differences at *p* < 0.05.

#### **4. Discussion**

The exogenous application of salicylic acid, proline, and yeast extract previously exhibited variable effectiveness levels on plant performance and tolerance to the harmful impact of salinity stress; therefore, the present study assessed and compared the effectiveness levels of these stress tolerance inducers on the growth characteristics, antioxidant levels, physiological and biochemical parameters, and yield of sweet pepper plants (*Capsicum annuum* L.) grown under the same saline conditions in order to determine which stress tolerance inducer should be recommended for the enhancement of crop performance and tolerance. In the present study, salt stress significantly decreased the aforementioned physiological parameters of sweet pepper plants. This reduction in chlorophyll a and b concentrations could be due to the effect of salinity on chlorophyll-degrading enzyme (chlorophyllase) activity, which reduces the chlorophyll synthesis level or negatively affects the structure and number of chloroplasts [40–42]. The chloroplast is one of the most vital organelles for photosynthesis and plant production, and is

dramatically affected by abiotic stresses [43,44]. The obtained results indicated that foliar applications of SA (1 mM), yeast extract (6 g L−1), and proline (10 mM) led to increased chlorophyll a and b contents in salt-stressed plants. These findings are in harmony with those obtained by Saleh et al. [25], Abdelaal et al. [21], and Soliman et al. [2]. These results might be due to the antioxidant scavenging influence of SA on chlorophyll degradation under saline conditions [45,46]. Correspondingly, the role of yeast extract in chlorophyll concentration enhancement might be due to the fact that yeast is rich in many essential elements, vitamins, and amino acids, which improve chlorophyll concentrations under stress conditions [21]. Moreover, our results showed that proline application minimized the harmful effects of all salinity levels on chlorophyll a and b concentrations due to its ability to function as a scavenger for ROS. Thus, proline plays a pivotal role in enzyme activation and protects chlorophyll from degradation under salt-stress conditions [47,48].

Additionally, relative water content (RWC) was decreased under salt stress. This decrease may be due to the reduction in water uptake [49] and/or its harmful effect on cell wall structure [50,51]. In contrast, RWC was significantly increased in stressed plants treated with SA, yeast, and proline. The ameliorative effects of these treatments on RWC could be due to the increase in osmoregulators, as well as to osmotic adjustment in plant cells [23,52,53].

Salt stress causes adverse effects on sweet pepper plants, including increased electrolyte leakage percentage. This increase may be due to the damaging effects on plasma membrane and selective permeability resulting in an increase in electrolyte leakage. This result is similar to that obtained in [54,55]. Conversely, the foliar application of SA, yeast, and proline led to decreased electrolyte leakage levels in all treatments. This beneficial effect could be due to the protective role of SA, yeast, and proline in plasma membrane stability and increasing soluble metabolite accumulation. A similar result was indicated by Ishikawa and Evans [56] and Huang et al. [57], who reported that osmoregulators improve plant growth and yield under various stress conditions.

Proline concentration was significantly increased in response to salt-stress conditions. This increment represents an important mechanism to minimize the deleterious impact of salinity stress and enhance plant growth [58]. The foliar application of SA, yeast, and proline under salt conditions may minimize the destructive effect of salinity on plant growth and improve proline accumulation. Similarly, SA application led to improved plant growth characteristics in maize plants under salt conditions [59]. Our results are in agreement with those of Huang et al. [60], Li et al. [61], and Gharsallah et al. [62]. Lipid peroxidation as MDA is an important factors indicating oxidative damage induced by salt stress. Lipid peroxidation was significantly boosted in salt-stressed (2000 and 4000 ppm) sweet pepper plants. Nonetheless, lipid peroxidation content was significantly decreased upon the foliar application of SA, yeast, and proline. These results may be attributed to the pivotal role of these treatments in decreasing oxidative stress damage, and consequently, in causing MDA reduction [25,48,53].

In the current study, superoxide and hydrogen peroxide, which are indicators of oxidative stress, were significantly produced in sweet pepper plants treated with NaCl at 2000 and 4000 ppm. This increase in superoxide and hydrogen peroxide production may be due to the fact that reactive oxygen species have a critical role under stress conditions in adjusting development, differentiation, redox levels, and stress signaling in the chloroplasts, mitochondria, and peroxisomes of plant cells [63,64]. Moreover, the high levels of hydrogen peroxide and superoxide are the main reasons for oxidative stress in the plant cells exposed to various stresses. Our results are supported by the findings of previous studies [65–67]. The application of SA, yeast, and proline on salt-stressed sweet pepper plants led to reductions in the formation of superoxide and hydrogen peroxide. This effect may be due to the role of these treatments in stabilizing protein structures and maintaining the redox states of plant cells, as well as stimulating antioxidant enzymes system [6,21,48]. Under salt stress (2000 and 4000 ppm), antioxidant enzyme activities were significantly increased in sweet pepper plants in order to combat the harmful impact of salt by adjusting osmotic balance. In agreement with our findings, similar results were noted in various plants under saline and drought conditions [68,69]. The activation of CAT and POX enzymes under salt conditions plays a key role in the improvement of plant defense systems. In the current study, the exogenous foliar application of SA, yeast, and proline led to improved antioxidant enzymes activity, as well as guarding the plant cells against oxidative stress and dehydration of the plasma membrane under salt-stress conditions. These results were supported by the findings reported in various plants [70–72].

The reductions of fruit numbers per plant, fruit fresh weight, and total fruit yield (ton hectare<sup>−</sup>1) under salt conditions are possibly due to the adverse impacts of salinity on the growth characteristics and physiological processes such as water uptake, photosynthesis, flowering, and fruit formation, which led to diminished yields. Accordingly, the highest level of salt (4000 ppm) was adversely more effective than the lowest one (2000 ppm). The same trends of salt stress were previously described in faba bean [73] and strawberry plants [74]. Our results indicate that proline treatment was the best, followed by SA and yeast treatments. This useful effect of proline may be due to its pivotal role in osmotic regulation, enzyme activation, and protein synthesis, which consequently enhances the growth and yield characteristics of stressed plants [47,75,76]. Also, SA plays an essential role as a stress tolerance inducer via reducing the oxidative damage and enhancing plant productivity under salt stress. These results are in harmony with previous findings of Gupta and Huang [77], Ahanger et al. [78], and Husen et al. [79].

#### **5. Conclusions**

According to our findings, salt stress caused significant decreases in chlorophyll concentrations, relative water content, and fruit yields. However, lipid peroxidation, proline, electrolyte leakage, and reactive oxygen species were increased. Based on the results, the foliar application of salicylic acid (1 mM), yeast extract (6 g L<sup>−</sup>1), and proline (10 mM) was an effective method by which to overcome the injurious effects of salt stress on sweet pepper plants. It may be concluded that relative water content and chlorophyll concentration, as well as antioxidant enzyme activity, were significantly modulated in the stressed treated sweet pepper plants. In contrast, electrolyte leakage and lipid peroxidation were decreased in treated sweet pepper plants under salt conditions. Thus, the application of salicylic acid, yeast, and proline led to a decrease in the harmful effects of salt stress by regulating osmolytes and antioxidants, which ultimately enhances the growth characteristics and fruit yields of sweet pepper plants. Interestingly, proline proved to be the best treatment for the further enhancement of plant performance and tolerance to salinity stress.

**Author Contributions:** K.A.A., L.M.E.-M., H.E., Y.M.H., E.I.I., M.E.-B., M.E.-E. and A.E. performed the experiments, analyzed the data, and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by king Saud University funding through Researchers Supporting Project number (RSP-2019/118).

**Acknowledgments:** The authors extend their appreciation to king Saud University funding through Researchers Supporting Project number (RSP-2019/118) and also to all members of PPBL and EPCRS excellence center, Fac. of Agric., Kafrelsheikh University.

**Conflicts of Interest:** The authors have declared no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Screening of Provitamin-A Maize Inbred Lines for Drought Tolerance: Beta-Carotene Content and Secondary Traits**

#### **Aleck Kondwakwenda \*, Julia Sibiya, Rebecca Zengeni, Cousin Musvosvi and Samson Tesfay**

School of Agriculture, Earth and Environmental Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville 3209, Pietermaritzburg, South Africa; sibiyaj@ukzn.ac.za (J.S.); zengeni@ukzn.ac.za (R.Z.); cousinmusvosvi@gmail.com (C.M.); tesfay@ukzn.ac.za (S.T.)

**\*** Correspondence: alekondwa@gmail.com; Tel.: +27-332-606-246

Received: 17 September 2019; Accepted: 21 October 2019; Published: 29 October 2019

**Abstract:** Provitamin A maize (*Zea mays* L.) biofortification is an ideal complementary means of combating vitamin A deficiency (VAD) in sub-Saharan Africa where maize consumption is high coupled by high VAD incidences. However, drought remains a major abiotic constraint to maize productivity in this region. Comprehensive drought screening of initial breeding materials before advancing them is important to achieve genetic gain. In this study, 46 provitamin-A inbred lines were screened for drought tolerance in the greenhouse and field under drought and optimum conditions using β-carotene content (BCC), grain yield (GY), and selected morphophysiological and biochemical traits. The results revealed that BCC, morphophysiological and biochemical traits were effective in discriminating among genotypes. Number of ears per plant (EPP), stomatal conductance (Gs), delayed leaf senescence (SEN), leaf rolling (RL), chlorophyll content (CC) and free proline content (PC) proved to be ideal traits to use when indirectly selecting for GY by virtue of having relative efficiency of indirect selection values that are greater than unity and considerable genetic variances under either or both conditions. The findings of this study form the basis of initial germplasm selection when improving provitamin A maize for drought tolerance.

**Keywords:** Provitamin A; maize; drought; morphological; physiological; biochemical; β-carotene

#### **1. Introduction**

Provitamin A maize has orange or yellow endosperm, which contains precursors of vitamin A in the form of carotenoids, hence the name "provitamin A maize". Provitamin A carotenoids include β-carotene, α-carotene and β-cryptoxanthin. Beta-carotene is the most important of the three carotenoids because it has higher provitamin A activity owing to its unique double ring molecular structure [1]. However, the ordinary (not improved) yellow maize grown and consumed throughout the world has β-carotene content of less than 1.5 μg g−<sup>1</sup> [2], which is too low given that the required target of total provitamin A carotenoid is 15 μg g−<sup>1</sup> [3]. Developing provitamin A maize cultivars with higher levels of provitamin A carotenoids through biofortification is a sustainable, cheap and effective complementary solution to VAD challenges faced by many developing nations [4].

An online Biofortification Priority Index (BPI) tool, developed and managed by HarvestPlus, shows that maize provitamin A biofortification as VAD intervention is most suitable for maize consuming developing countries, particularly the Southern Africa region (www.harvestplus.org/knowledgemarket/BPI). However, maize production in this region is vulnerable to drought due to recurring low annual precipitation coupled by poor coping capacity of most farmers. For instance, in 1992 and 2002, most of the southern African countries experienced the worst droughts resulting in over 60% maize yield loss in the whole region [5]. In 2013, 770 million people in the Southern African Development

Community (SADC) were at risk of food insecurity due to severe mid-season dry spells [6]. In 2016, eight of South Africa's nine provinces were declared food-insecure due to drought [7].

Drought stress affects maize at almost all growth stages, but the flowering and grain filling stages are the most susceptible, with yield losses of over 90% reported when drought coincides with these growth stages [8]. Genetic improvement of maize for drought tolerance through breeding is a sustainable solution to reduce the impacts of drought. However, breeding for drought tolerance is a complex task because the trait is controlled by many genes, and is highly affected by genotype and environment interaction. Furthermore, grain yield, which is the trait of interest, has very low variation and heritability under drought conditions, which makes selection difficult [9]. Comprehensive screening forms the foundation of any successful drought tolerance breeding program [10]. Additionally, screening materials for drought tolerance at the initial stages of breeding often imparts tolerance to other related stresses such as low nitrogen stress [11].

To increase the chances of selecting the appropriate genotypes, breeders should meticulously consider all the available information at screening phase. Screening maize for drought tolerance entails the selection of high yielding genotypes under water deficit stress and/or optimum conditions. Indirect selection for grain yield via related secondary traits helps to circumvent the challenge of poor grain yield variation and low heritability under drought conditions [12]. Stress-tolerant indices and multivariate statistics are often useful when selecting best genotypes.

Indirect selection for drought tolerance involves selecting for multi-secondary traits that are highly correlated to grain yield and have high heritability values. Furthermore, an ideal secondary trait should have efficiency of indirect selection relative to direct selection of greater than a unity [13]. Morphological and physiological (morphophysiological) traits that are associated with maize drought tolerance include anthesis-silking interval (ASI), leaf rolling (RL), chlorophyll content (CC), leaf senescence (SEN), and number of ears per plant (EPP) [12,14]. Stomatal conductance (Gs) analysis as a physiological response to drought stress has not been widely applied as a screening criterion for drought tolerance in maize and therefore, information about its correlation with grain yield and heritability under drought stress and optimum conditions is not well established.

Biochemical changes that are drought-induced in plants include increase in stress signaling hormones and proteins regulators such as proline and abscisic acid (ABA) among others [15,16]. Proline is an amino acid, which plays an osmoregulatory role in plants under drought conditions [17]. Despite the presence of genetic variation of proline content in plants under drought stress and wide application of proline analysis in understanding drought tolerance of other crops such as wheat [18], cowpea [19], and peanut [20], it has not been applied in large-scale maize drought tolerance screening.

Given the importance of maize provitamin A biofortification in maize-consuming developing countries and the prevailing devastating impacts of drought to maize productivity, it is important to investigate the effectiveness of integrated application of morphophysiological and biochemical traits in screening provitamin A maize inbred lines for drought tolerance. The objectives of the study were to: (i) determine the level of genotypic variation for drought tolerance among tropical provitamin A maize inbred lines in regards to secondary traits, (ii) screen candidate lines for drought tolerance based on grain yield and β-carotene content, and (iii) identify secondary traits that can be effectively used for indirect selection of grain yield under drought stressed and non-stressed conditions. The study forms the basis of germplasm selection for use in provitamin A maize drought tolerance breeding programs.

#### **2. Materials and Methods**

#### *2.1. Plant Materials and Study Sites*

Fifty inbred lines were screened for drought stress tolerance under managed drought conditions. Inbred lines consisted of 46 provitamin A (orange endosperm) and four drought tolerant white endosperm (non-provitamin A) checks. The 46 provitamin-A inbred lines were sourced from the provitamin A biofortification nurseries of International Maize and Wheat Improvement Center

(CIMMYT) (33) and International Institute of Tropical Agriculture (IITA) (13). The four drought tolerant checks were obtained from the Agricultural Research Council (ARC), Grain Crops Institute, Potchefstroom, South Africa. Names, codes and other information about the inbred lines are given in supplementary Table S1. The respective institutions could not provide pedigree information. The study was carried out across four environments (Env), which were two greenhouse and two field in KwaZulu-Natal Province of South Africa. Greenhouse trials were carried out at the University of KwaZulu-Natal (UKZN), Agriculture Pietermaritzburg campus (29◦46 S, 30◦58 E) from January to April 2017 (Env1) and May to September 2017 (Env2). The field trials were carried out at Ukulinga Research Farm in Pietermaritzburg (29◦40 S, 30◦24 E) from April to August 2017 (Env3) and at Makhathini Research station, Jozini, South Africa (27◦39 S, 32◦17 E) from May to September 2018 (Env4). Supplementary Table S2 shows monthly weather data for the four environments during the respective growing periods.

#### *2.2. Experimental Design and Crop Establishment*

An alpha lattice design was used to screen the 50 genotypes with two replications containing five incomplete blocks with ten genotypes each and two water regimes (water stress, S and optimum conditions, W) across all the four environments. In the field, the plot size was two rows of 5 m with 0.75 m between the rows and intra row spacing of 0.30 m. Plots were planted with two seeds per station and thinned to one plant 2 weeks after crop emergence. In the greenhouse, a plot was made of four 5 L perforated plastic pots with two plants in each pot, which were thinned to one plant per pot 2 weeks after crop emergence. Pine-bark growing media mixed with loam soil at a ratio of 3:1 was used in the greenhouse. In the field, the soil was predominantly black clay loam soil at both sites. The water stress treatment (S) for all the experiments was implemented in accordance with CIMMYT protocols of withholding irrigation at two weeks prior to expected anthesis date [21]. The water stress condition was maintained until 5 weeks after the flowering of 50% of the genotypes then a single irrigation was applied at grain filling stage. In the field, the optimum treatment (W) involved a 10-day interval sprinkler irrigation throughout the growing period. In the greenhouse, W consisted of drip irrigation for 3 min, four times per day. Across all the environments, compound fertilizer was applied at the rate of 150 kg N, 65 kg P and 65 kg K ha−<sup>1</sup> at the time of planting and top-dressing fertilizer was applied at five weeks after emergence at a rate of 60 kg N ha−1. In the field, weeds were controlled using Gramoxone® (Syngenta, Greensboro, NC, USA) at a rate of 5 L ha−<sup>1</sup> and manual weeding whilst hand weeding was done in the greenhouse. Coragen® (Dupont, Washington, DC, USA) and Karate® (Syngenta, Greensboro, NC, USA) insecticides were used to control insects at a rate of 1 L ha−<sup>1</sup> across all environments when it was necessary.

#### *2.3. Plant Characteristics*

#### 2.3.1. Morphophysiological Traits

Data for the following morphophysiological and biochemical traits were collected under both water regimes across all the four environments. Most traits were measured following a CIMMYT protocol described by [21]. Both field and greenhouse grain yield (GY) was estimated per plot area bases and converted to tonnes per hectare. The plot area in the greenhouse was determined by multiplying the cylinder-shaped pot area by 4 since a plot was made up of four pots. Similarly, the plot area for the field experiments was determined by calculating the area of a rectangular shaped plot. Number of ears per plant (EPP) was computed as the number of ears with at least one fully developed grain divided by the number of harvested plants. Days to anthesis (DA) is the number of days after planting to when 50% of the plants in a plot start to shed pollen. Days to silking (DS) is the number of days after planting when 50% of plant in a plot produces silks. Anthesis-silking interval (ASI) was calculated as DS minus DA. Leaf rolling (RL) was scored using a scale of 0 to 10 where 0 = unrolled leaf and 10 = leaf rolled like an onion and scores were converted to percentage (%) with measurements

taken twice after imposing drought treatment. Leaf senescence (SEN) was scored using a scale from 1 to 10 (1 = 10%; 2 = 20%; 3 = 30%; 4 = 40%; 5 = 50%; 6 = 60%; 7 = 70%; 8 = 80%; 9 = 90%; and 10 = 100% dead leaf area) at 3, 5 and 7 weeks after 50% of the plants reached anthesis. Chlorophyll content (CC) was measured from the adaxial surface of the second top fully expanded leaf of five plants per plot at 3, 5 and 7 weeks after 50% of the plants reached anthesis using SPAD-502-Plus chlorophyll meter (Konica Minolta, Osaka, Japan). Stomatal conductance (Gs) was measured from the abaxial surface of the second top fully expanded leaf using a SC–1 leaf porometer (Decagon Devices®, Pullman, WA, USA) at 3, 5 and 7 weeks after 50% of the plants reached anthesis. Chlorophyll content (CC) and stomatal conductance were measured at midday periods (1200–1400 h). Beta carotene content (BCC) and proline (PC) content were determined and measured as described in the following section.

#### 2.3.2. Biochemical Traits

Beta-carotene content was measured from kernels harvested from the self pollinated plants of each of the 46 provitamin inbred lines. A sample of 20 g, which contained about 30 to 50 kernels, was randomly collected from each of the 46 provitamin-A inbred lines and dispatched to Agricultural Research Council (ARC), Science Analytical Laboratory, Pretoria, South Africa (http://www.arc.agric.za) for β-carotene analysis. High-performance liquid chromatography (HPLC) was used for analysis following a protocol for dried maize kernels as described by [22]. The β-carotene analysis was done three times per sample, giving three data points per each genotype.

Proline analysis was performed at UKZN, Crop Science laboratory following a protocol by [23]. Fresh leaf samples were collected from the second top fully expanded leaves for the S and W treatments of both field and greenhouse experiments at 3 weeks after imposing the S treatment. The leaf samples were freeze-dried at very low temperature (−74 ◦C) using liquid nitrogen before grinding them into fine powder. A 0.5 g ground leaf sample was homogenized in 10 mL of 3% aqueous sulfosalicylic acid and the homogenate was filtered. Two ml of filtrate was mixed with 2 mL acid-ninhydrin and 2 mL of glacial acetic acid for 1 h in a water bath at 100 ◦C. After cooling, 4 mL of toluene were added and then mixed vigorously using a rotor. The top mixture containing proline within toluene was decanted from the aqueous phase then taken for UV visible spectrophotometer analysis for the absorbance measurement at a wavelength of 520 nm using a model UV-1800 spectrophotometer (Shimadzu Corporation, Kyoto, Japan). The proline concentration was calculated using the formula shown in Equation (1) [23].

$$\begin{array}{l} \text{Proline Content} \quad (\text{\(\mu\) per gram of dry leaf tissue})\\ \qquad = [(\text{\(\mu\) produce}/\text{mL}) \times \text{mL toture}) / 115.5 \text{ \(\mu\)}\\ \qquad / \text{\(\mu\) mole} / [(\text{\(g\) sample}) / 5] \end{array} \tag{1}$$

where 115.5 is the molecular weight of proline.

#### *2.4. Data Analysis*

#### 2.4.1. Analysis of Variance, Mean Performance and Stress-Tolerant Index

Analysis of variance (ANOVA) for all morphophysiological and biochemical traits was done after carrying out a test of homogeneity of variances. A lattice procedure of R software version 3.5.1 [24] was used to carry out the ANOVA following a mixed model (Equation (2)). Genotypes were considered as fixed effects and environments as random:

$$\mathcal{Y}\_{i\bar{j}k\bar{l}m} = \mu + R\varepsilon\_{l} + B(R\varepsilon)\_{ij} + G\_{k} + E\_{l} + \mathcal{W}\_{m} + GE\_{kl} + G\mathcal{W}\_{km} + E\mathcal{W}\_{lm} + GE\mathcal{W}\_{k\bar{l}m} + \mathcal{E}\_{i\bar{j}k\bar{l}m} \tag{2}$$

where *Yijklm* is the trait of interest, μ is the mean effect, *Rei* is the effect of the *i*th replicate, *B*(*Re*)*ij* is the effect of the *j*th incomplete block within the *i*th replicate, *Gk* is the effect of the *k*th genotype, *El* is the effect of the *l*th environment, *Wm* is effect of the *m*th water regime while *GEkl*, *GWkm*, *EWlm* and *GEWklm* are the respective interactions and E*ijklm* is the residual error term.

Beta-carotene was analysed separately using the general linear model given in Equation (3).

$$Yij = \mu \, +Gj + S(G)ij \,\tag{3}$$

where *Yij*—the performance of *i*th sample of the *j*th genotype; *Gj*—the effect of *j*th genotype; *S*(*G*)*ij*—sample within genotype, which is the error term.

Least significant difference (LSD) test was carried out at 0.5 α level to separate the means.

Stress tolerance index (STI) described by [25] was used to select high yielding genotypes under both conditions as shown in Equation (4):

$$STI = \frac{\text{Y} \dot{p}i \times \text{Y} si}{\text{Y}p^2} \tag{4}$$

where *Ys*—grain yield of a genotype *i* under drought-stressed condition; *Yp*—grain yield of *i* genotype under non-stressed condition, and *Xp*—mean yield of genotypes under non-stressed condition. Genotypes with high STI value and β-carotene content greater than 1.5 μg g−<sup>1</sup> were selected.

#### 2.4.2. Variance Components and Heritability

Variance components and broad sense heritability (H) were computed for morphophysiological and biochemical traits in R software following a procedure described by [26] as shown in Equations (5) and (6).

$$
\delta^2 p = \delta^2 \mathcal{g} + \frac{\delta^2 \mathcal{g} \mathcal{e}}{\mathcal{e}} + \frac{\delta^2}{m} \tag{5}
$$

$$H = \frac{\delta^2 g}{\delta^2 p} \tag{6}$$

where δ2*p*—phenotypic variance, δ<sup>2</sup> *g*—genotypic variances, δ<sup>2</sup> *ge*—genotype by environment interaction variance, δ2—error variance, *r*—number of replications, *e*—number of environments and *H*—broad sense heritability.

Heritability classifications guidelines described by [27] were used to describe the heritability levels exhibited by the measured traits in this study in which values from 0 to 0.3 was low, 0.3 to 0.6 was moderate and >0.6 was high.

#### 2.4.3. Principal Component Biplot, Trait Correlations and Relative Selection Efficiency

Phenotypic and genotypic correlations among morphophysiological and biochemical traits were computed separately for each water regime using the META-R software version 6.04 [28]. Correlation classifications guidelines described by [29] were used to explain the correlations among the traits in this study in which correlation coefficients values from ±0.9 to ±1.00 were considered very high correlations, ±0.7 to ±0.9 were high, ±0.5 to ±0.7 were moderate, ±0.3 to ±0.5 were low and ±0.00 to ±0.3 were negligible. A combined principal component PCA biplot was computed to graphically show the traits associated with each water regime. Equation (7) was used to test the efficiency of indirect selection of grain yield via secondary traits relative to direct selection as outlined by [13].

$$\text{Relative efficiency of indirect selection} = \frac{|r\_{\mathcal{S}}|h\_{\mathcal{X}}}{h\_{\mathcal{G}\mathcal{Y}}} \tag{7}$$

where *rg* is the value of the genotypic correlation between GY and a secondary trait, *hX* is the square root of the broad sense heritability of trait, and *hGY* is the square root of the broad sense heritability of grain yield. According to [13] the most desirable secondary traits should have an efficiency of indirect selection relative to direct selection of greater than unity.

#### **3. Results**

#### *3.1. Analysis of Variance, Mean Performances and Stress-Tolerant Index*

A combined ANOVA was carried out after separate ANOVA had shown significant (*p* < 0.05) effects of genotype, water regime and their interaction for the studied traits (Table 1). The combined ANOVA revealed that the genotype, water regime, environment and their respective interactions had significant (*p* < 0.001 and *p* < 0.05) effects on GY and other traits except DA, EPP and PC, which were not significantly affected by the interaction of the environment and water regime. Genotypes exhibited significantly (*p* < 0.001) different mean BCC (Table 1).

Mean performance of the top ten and bottom five genotypes with >1.5 μg g−<sup>1</sup> BCC and ranked in descending order of STI values are presented in Table 2. Thirty-one genotypes had >1.5 μg g−<sup>1</sup> BCC. Mean BCC was 2.05 μg g−<sup>1</sup> with genotype 24 (CLHP0022) ranked first in the STI ranking of genotypes with >1.5 μg g−<sup>1</sup> BCC. Genotypes 1 (CLHP00306) and 11 (CLHP00378) had the highest BCC of 4.22 μg g−<sup>1</sup> whilst 14 (CLHP0343) had the lowest value of 0.64 μg g<sup>−</sup>1. The mean performance of all the 50 genotypes under the two water regimes are given in supplementary Table S3. Mean STI was 0.53. Fifty percent of genotypes surpassed the mean STI value. The highest STI value was 0.94 exhibited by entry 27 (CLHP0005) whilst entry 20 (CLHP0364) had the least STI value of 0.23. Genotype 50 (CML569) had the highest STI value of 0.71 among the four drought tolerant checks. The proportion of test genotypes that were ranked higher than the best check on the STI ranking was 10.86%.

The mean GY under optimum and stress conditions were 1.72 t ha−<sup>1</sup> and 0.88 t ha<sup>−</sup>1, respectively. This resulted in 51.2% mean yield loss due to drought stress with the highest and lowest percentage yield losses of 64.17% and 29.21%, respectively. Genotype 50 (CML569) was the best yielding check under both drought-stressed and non-stressed conditions with 0.94 t ha−<sup>1</sup> and 2.26 t ha<sup>−</sup>1, respectively. The proportion of test genotypes that yielded higher than the best check was 45.57% and 8.70% under drought stressed and non-stressed conditions, respectively. Mean DA was reduced from 75.03 days under optimum conditions to 69.40 days under water stress conditions. Mean ASI increased from a mean of 1.97 days under optimum conditions to a mean of 8.56 days under drought stress conditions. Mean EPP was reduced by drought stress from 2.24 under optimum to 1.68 under drought stress conditions.

Mean RL increased from 3.04% under optimum conditions to 49.78% under water stress. Stomatal conductance was severely reduced from a mean of 368.94 mmol m−<sup>2</sup> s−<sup>1</sup> under optimum conditions to a mean of 49.78 mmol m−<sup>2</sup> s−<sup>1</sup> under drought stress conditions. Leaf senescence increased due to drought stress from a mean of 11.60% under optimum conditions to 50.18% under water stress conditions. Proline content increased from a mean of 31.83 μg g−<sup>1</sup> under optimum conditions to a mean of 149.23 μg g−<sup>1</sup> under water stress conditions with genotype 31 (CML486) having the highest PC of 230.63 μg g−<sup>1</sup> under water stress conditions.


**Table 1.** Mean squares and significant tests after combined analysis of variance for nine morphophysiological traits of 50 inbred

lines.

ASI—anthesis-silking interval, DA—days to anthesis, EPP—ears per plant, RL—leaf rolling, SEN—leaf senescence, Gs—stomatal conductance, CCI—chlorophyll content PC—prolinecontent,\*\*\*\*\*\*0.001and

 BCC—β-carotene content, *p* < 0.05, *p* < 0.01, *p* < ns, not significant.


#### *Agronomy* **2019** , *9*, 692

S—Water stress conditions and

W—well-watered

 condition.

#### *3.2. Heritability and Variance Components*

Grain yield, DA, EPP, Gs and CC exhibited high heritability and genotypic variance values under non-stressed conditions whilst ASI, RL, SEN and PC had high values under drought conditions (Table 3). Heritability estimates ranged from 0.250 to 0.998 under non-drought conditions and 0.127 to 0.987 under drought conditions. Most traits, except Gs, were characterized by sharp differences in heritability values between the two water regimes. For instance, heritability estimates for GY yield dropped from 0.621 under non-stressed conditions to 0.398 under drought conditions whilst PC had heritability values of 0.263 and 0.777 under non-stressed and drought conditions, respectively. Generally, variance due to genotype by environment interaction was higher under drought conditions than under non-drought conditions.


**Table 3.** Variance components and heritability of measured morphophysiological traits.

#### *3.3. Principal Component Biplot Analysis*

A combined principal components biplot was constructed to show traits that were more outstanding in discriminating among genotypes per each water regime (Figure 1). First and second principal components explained 80.4% of the total variation. Performance of genotypes with respect to the measured nine traits under water-stressed conditions (prefixed with SG) were located on the negative side of the biplot whilst majority of genotypes under optimum conditions (prefixed with WG) were on the positive side of the biplot. Traits ASI, RL, SEN and PC were more discriminating among genotypes under drought conditions whilst GY, Gs, CC, DA and EPP were more discriminating among genotypes under non-stressed conditions. Vector length is relative to the discriminating power of the respective trait. Hence, the trait discriminating power within the non-stressed condition in descending order was EPP, DA, GS, CC and GY. On the other hand, PC, RL, SEN and ASI, was the respective order of discriminating under drought stress conditions. Genotypes positioned at or close to a vector of a traits are more associated with that trait. Also genotypes that are located at the tip end of trait vector excelled in the respective trait. For instance, genotype 39 (TZM1224) excelled in proline content (PC) under drought conditions whilst genotypes 2 (CLHP00306) and 25 (CLHP0113) were further on the GY vector under well-watered conditions. Very few genotypes were associated and excelled in EPP and DA under non-stress conditions.

**Figure 1.** A combined principal component biplot showing genotypes clustering under stress conditions and well-watered conditions. ASI—anthesis-silking interval, CC—chlorophyll content, DA—days to anthesis, EPP—ears per plant, GY—grain yield, Gs—stomatal conductance, RL—leaf rolling, SEN—leaf senescence, PC—proline content, SG—genotypes under water stress conditions, WG -genotypes grown under optimum conditions.

#### *3.4. Phenotypic Correlation Analysis*

Phenotypic correlation coefficients (*r*) among measured morphophysiological and biochemical traits under non-stressed conditions (upper diagonal) and drought conditions (lower diagonal) (Table 4). Under drought stress, GY had significant (*p* < 0.001) and positive correlations that were high with EPP, and low with RL, CC and PC. It also had a significant and negative correlation that was moderate with ASI and Gs, and low with DA, and SEN. Number of ears per plant had significant (*p* < 0.001) and negative correlations which were high with ASI, moderate with DA and Gs, and low with SEN, CC and PC. Stomatal conductance had significant (*p* < 0.001) and moderate positive correlations with RL and SEN. In addition, PC had a significant (*p* < 0.001,) low positive correlation with RL and SEN. It also had a significant (*p* < 0.05), low negative correlation with CC. Under non-stressed conditions, GY had significant (*p* < 0.001) and positive correlation, which was moderate with CC and EPP, low with Gs and DA, and negligible with RL. It also had significant and negative correlation, which was moderate with ASI, low with SEN and negligible with PC. Number of ears per plant had a significant (*p* < 0.001) negative correlation which was moderate with ASI and low with SEN. On the other hand, EPP had a significant (*p* < 0.001) positive correlation, which was moderate with Gs and low with CC. Proline content had negligible correlation with SEN and CC under non-stressed conditions.


PhenotypiccorrelationcoefficientsdescribingassociationoftraitsunderS(lowerdiagonal)andW(upperdiagonal)conditions.

#### *Agronomy* **2019** , *9*, 692

SEN—leaf senescence, \*,\*\*,\*\*\* indicate level of significance of the correlation at *p* < 0.05, *p* < 0.01 and *p* < 0.001, respectively.

Across environments phenotypic correlation coefficients for grain yield only were computed in order to compare greenhouse and field environments (supplementary Table S4). Environments were further subdivided into greenhouse non-stressed, greenhouse-stressed, field-non-stressed and field-stressed. Correlations among all the four environmental subdivisions ranged from moderate to high (*r* = 0.456 to *r* = 0.751). Significant (*p* < 0.001) and high correlation (*r* = 0.751) was observed between greenhouse and field-non-stressed whilst greenhouse and field-stressed had moderate correlations (*r* = 0.643). Significant (*p* < 0.001) moderate correlations were also observed between greenhouse non-stressed and stressed (*r* = 0.553), and field-stressed and non-stressed (*r* = 0.585). Greenhouse non-stressed and field-stressed had a significant (*p* < 0.05) moderate correlation (*r* = 0.502) whilst greenhouse stressed and field-non-stressed had a significant (*p* < 0.05) moderate correlation (*r* = 0.456).

#### *3.5. Relative E*ffi*ciency of Indirect Selection*

Relative efficiency of indirect selection through secondary traits for grain yield ranged from 0.142 for DA to 1.370 for Gs under drought conditions and from 0.012 for PC under to 1.235 for Gs under non-drought conditions (Table 5). Traits that exhibited relative selection efficiency of greater than unity are EPP, RL, Gs, SEN and PC under drought conditions, and EPP, Gs and CC under optimum conditions. All the genetic correlations coefficients used in calculating the relative efficiency of indirect selection are given in supplementary Table S5.

**Table 5.** Genetic correlations and the relative efficiency of indirect selection through secondary traits for grain yield improvement under drought stress and optimum conditions.


\* *p* < 0.05, \*\*\* *p* < 0.001, ns, not significant.

#### **4. Discussion**

Comprehensive screening of germplasm forms the basis of selecting the right materials, which in turn increases the chances of achieving genetic gain in plant breeding. In this study, the ANOVA showed highly significant genotypic differences with respect to all the studied traits. This indicates the feasibility of genetic improvement for drought tolerance and β-carotene content given that genetic variation is a prerequisite for genetic advance [30]. In a related study, [31] also reported a significant genetic variation among inbred lines obtained from CIMMYT and IITA stress-tolerant breeding programs. The significance of the water regime and its interaction with the genotype for all the traits indicate that the imposed drought stress was effective in discriminating among the genotypes. Furthermore, the highly significant genotype by environment interaction effects observed for most of the traits showed that there were performance differences by inbred lines across the environments. This confirms the complexity of drought tolerance breeding as this makes performance ranking difficult [12]. However, the significant high and moderate positive correlations observed between greenhouse and field environments after subdividing them along water regimes suggest that the two screening environments had almost similar discriminating abilities. This validates the ranking of genotypes within water regimes across the four environments.

Trait performance inconsistence was observed between the two water regimes. That is, combined principal component biplot revealed that RL, SEN, PC and ASI were more discriminating under drought conditions while EPP, GY, CC, Gs and DA were more discriminating under non-stressed conditions. Similarly, variables Gs, GY, EPP, DA and CC had higher heritability and genotypic variance estimates under optimum conditions than under drought conditions whilst the opposite pattern was observed for ASI, SEN, PC and RL (Table 2). In this regard, our findings concur with [32], who reported a decline in the genotypic variances for GY and EPP whilst that of ASI increased under drought conditions. This inconsistency justifies the importance of using a combination of multi-traits when screening germplasm for contrasting growing conditions [33]. This also suggests the need of separating the breeding programs to target different growing conditions in which different secondary traits would be utilized for selection. South Africa, like most SSA countries, has mixed growing conditions; therefore, separating the breeding target environments would be ideal.

The low heritability and genotypic variance estimates exhibited by GY under drought stress agree with findings by some of the previous researchers [34,35]. The observed yield loss (51.2%) due to drought stress was lower than the 81% reported by [36] who attributed greater variation under drought stress to kernel number. In this study, greater variation can be attributed to EPP, RL, Gs, SEN and PC under drought conditions, and EPP, Gs and CC under non-stressed conditions. Stomatal conductance was one of the traits that consistently accounted for greater variation under both conditions. This coupled with having a relative efficiency greater than unity suggests that Gs can be effectively used as selection proxy for GY [13]. This concurs with [37], who reported Gs as the major physiological trait that can effectively discriminate among genotypes between drought tolerant and susceptible plant genotypes. However, its low correlation with GY under non-stressed conditions reduces its effective utilization as an indirect selection criterion for production under optimum conditions.

The EPP and ASI are by far the most applied traits in maize drought tolerance studies [21,33,38]. In the current study, EPP was one of the largest contributors to the total genetic variation as demonstrated by high heritability and genotypic variance estimates under both drought stress and optimum conditions. Furthermore, by having higher discriminatory power as indicated by a longer PCA vector and high correlation with GY under both conditions confirms that EPP is an important trait in maize drought tolerance studies as reported by other researchers [9,36]. This, in addition to having a relative efficiency value greater than a unity justifies the use of EPP in indirect selection for GY. Despite having a moderate correlation with GY and moderate to high heritability estimates under both conditions, the effective use of ASI in indirect selection for GY would be limited by having a relative efficiency value of less than a unity. This is contrary to the findings by [32,39], who reported a relative efficiency of indirect selection for ASI, which was greater than a unity. However, it should be noted that although ASI did not demonstrate to be an ideal trait for indirect selection for GY, it could still be utilized in drought tolerant maize breeding by virtue of having moderate correlation with GY and moderate to high heritability estimates.

The observed significant correlations between GY and photosynthesis related traits such as SEN, RL, CC and Gs confirms the importance of photosynthesis for maize yield [40,41]. Leaf senescence and RL had relative efficiency values of greater than unity under drought and therefore can be used in indirect selection for GY. This concurs with [42] who found delayed leaf senescence to be useful in indirectly selecting for GY in maize under drought conditions. Moderate correlation between PC and GY coupled with high heritability estimate and relative efficiency of greater than unity under drought conditions, infers that genotypes that exhibited high PC under drought stress can be selected as drought tolerant. This supports the claim that under drought conditions proline is released to effect plant cell osmotic adjustments which helps to conserve cell turgor [17,43]. In a related study, [44] reported an increase in PC of 47% and 114% after exposing wheat genotypes to reproductive and grain filling drought stresses, respectively. However, lack of high correlation between PC and GY under non-stressed conditions hinders its effective use as a selection proxy for GY for optimum production.

Although this study did not investigate the effects of water regime, environment and the interaction of genotype and environment on BCC, there is enough evidence in literature that provitamin A content cannot be influenced by genotype and environment interaction but can be affected by the environment [45–47]. This makes selection for provitamin A an easy task. Furthermore, in their studies, [47] and [48] reported no significant correlation between GY and BCC, indicating that the two key traits can be improved simultaneously. The maximum BCC of 4.22 μg g−<sup>1</sup> observed in this study was lower than 13.22 μg g−<sup>1</sup> reported by [2]. This difference could be attributed to the natural superiority of temperate maize, which was studied by the later, over the its tropical counterpart [49].

#### **5. Conclusions**

This study showed that BCC, and the morphophysiological and biochemical traits applied in the screening of provitamin-A inbred lines for drought tolerance were effective in discriminating among the evaluated genotypes. There was considerably high genetic variation among the provitamin A genotypes under study that can be utilized when breeding for drought tolerance. The study also demonstrated that EPP, Gs, RL, SEN and PC can be effectively used in indirect selection of GY under drought-stressed conditions whilst EPP, Gs, and CC were ideal traits for GY indirect selection under non-stressed conditions. By having more secondary traits with a relative efficiency greater than unity under drought stress than under non-stressed, the study confirmed that indirect selection would be more useful than direct selection under drought conditions where GY exhibited low heritability and genotypic variance estimates. Consequently, based on BCC >1.5 μg g−1, which was greater than that of ordinary maize and the ranking for STI values, 30 inbred lines were selected to be utilized in developing drought tolerant provitamin A maize varieties.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4395/9/11/692/s1, Table S1: List of maize inbred lines used in the study, Table S2: Monthly weather data during the greenhouse trials (Env & Env2) at UKZN and field trials (Env3 and Env4) at Ukulinga and Makhatini research stations in South Africa, Table S3: Mean performance of all the genotypes evaluated under drought (S) and non-stressed (W) conditions across four environments ranked by the selection index (SI) values, Table S4: Pearson's correlation coefficients (r) describing association of four environments grouped according to water regimes, Table S5: Genetic correlation describing association of traits under S (lower diagonal) and W (upper diagonal) conditions.

**Author Contributions:** Conceptualization, A.K. and J.S.; Data curation, A.K. and C.M.; Formal analysis, A.K., C.M., and S.T.; Funding acquisition, J.S.; Investigation, A.K.; Methodology, A.K., J.S., C.M. and S.T.; Project administration, J.S.; Supervision, J.S. and R.Z.; Validation, A.K., J.S., R.Z., C.M. and S.T.; Writing—original draft, A.K. Writing—review & editing, A.K., J.S., R.Z., C.M. and S.T.

**Funding:** This research was funded by Alliance for a Green Revolution in Africa (AGRA) grant number 2014 PASS 013, The World Academy of Science (TWAS) and the National Research Foundation (NRF) of South Africa grant number SFH150914142525.

**Acknowledgments:** The authors would like to thank Andile Mshengu for administrative duties. CIMMYT and IITA are also acknowledged for providing the germplasm.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
