An investigation into the beneficial effects and molecular mechanisms of humic acid on foxtail millet under drought conditions

The aim of this study was to determine the effects and underlying molecular mechanisms of humic acid (HA) on foxtail millet (Setaria italica Beauv.) under drought conditions. The rainless climate of the Shanxi Province (37°42'N, 112°58'E) in China provides a natural simulation of drought conditions. Two foxtail millet cultivars, Jingu21 and Zhangza10, were cultivated in Shanxi for two consecutive years (2017–2018) based on a split-plot design. Plant growth, grain quality, and mineral elements were analyzed in foxtail millet treated with HA (50, 100, 200, 300, and 400 mg L-1) and those treated with clear water. Transcriptome sequencing followed by bioinformatics analysis was performed on plants in the normal control (CK), drought treatment (D), and drought + HA treatment (DHA) groups. Results were verified using real-time quantitative PCR (RT-qPCR). HA at a concentration of 100–200 mg L-1 caused a significant increase in the yield of foxtail millet and had a positive effect on dry weight and root-shoot ratio. HA also significantly increased P, Fe, Cu, Zn, and Mg content in grains. Moreover, a total of 1098 and 409 differentially expressed genes (DEGs) were identified in group D vs. CK and D vs. DHA, respectively. A protein-protein interaction network and two modules were constructed based on DEGs (such as SETIT_016654mg) between groups D and DHA. These DEGs were mainly enriched in the metabolic pathway. In conclusion, HA (100 mg L-1) was found to promote the growth of foxtail millet under drought conditions. Furthermore, SETIT_016654mg may play a role in the effect of HA on foxtail millet via control of the metabolic pathway. This study lays the foundation for research into the molecular mechanisms that underlie the alleviating effects of HA on foxtail millet under drought conditions.


Introduction
Foxtail millet (Setaria italica Beauv.) is a grain crop that grows in arid and semi-arid regions. It is an important crop in many areas of Africa and Asia due to its ability to grow in harsh environments [1], and is the dominant food crop in many provinces of China, including Shanxi

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0234029 June 2, 2020 1 / 14 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [1]. Shanxi Province has a typical continental arid climate with annual rainfall ranging between 400-650 mm, which is below standard irrigation conditions [2]. It is therefore of great practical importance to investigate the drought-resistant growth mechanisms of foxtail millet and identify ways to increase its yield in arid and semi-arid conditions. Humic Acid (HA) is a natural organic polymer compound and a major component of humus, which plays an important role in crop quality, yield, and stress resistance [3,4]. According to recent reports, HA affects plant quality in two major ways: helping to resist stress by controlling the amount of reactive oxygen species present [5], and promoting growth, photosynthesis, nitrogen assimilation, and amino acid metabolism [4,6]. Parađiković et al. found that using a mixture of biostimulants containing HA increased the yield of yellow pepper crops. Peptides and amino acids contained in this mixed preparation promoted high temperature resistance in the peppers, and stimulated root growth and development, while vitamins and HA supported fruit growth [7]. A study showed that high concentrations of HA had a significant positive effect on the morphological characteristics of cucumber, including plant height, leaf number, and fresh weight and yield. Additionally, the percentage of total chemical components [nitrogenium (N), phosphorus (P), potassium kalium (K), calcium (Ca), and magnesium (Mg)] in the leaves of cucumber plants increased as HA concentration increased [8]. Maji et al. showed that HA facilitates plant growth by improving the microbial community structure of soil and increasing mycorrhizal colonization in the roots of Pisum sativum [9], while another study found that foliar-applied HA improved dry weight and mineral nutrient uptake [such as manganese (Mn), copper (Cu), zinc (Zn), and calcium (Ca)] of maize [10]. HA has also been shown to stimulate nitrogen assimilation and amino acid metabolism in maize at both the physiological and molecular level [6]. Olaetxea et al. demonstrated that root hydraulic conductivity and aquaporin-related gene expression are important for plant shoot outgrowth induced by HA [11], and Kuşvuran et al. described the effects of different HA treatments on yield and performance in common millet. Together, these results indicate that HA treatment can significantly improve plant yield and quality [12]. However, despite all of these promising reports, the effects of HA on foxtail millet, and the underlying molecular mechanisms, remain unknown In the present study, we utilized the natural drought conditions in Shanxi Province to investigate the effects of HA on crop yield and quality of foxtail millet. The optimal HA concentration was determined by measuring growth and nutritional indicators in foxtail millet samples after different periods of HA treatment. Moreover, high-throughput sequencing analysis was used to assess changes in gene expression related to photosynthetic assimilation and nutrient indicators. This allowed the molecular mechanisms that underlie the effects of HA on dry matter and nutrient accumulation in foxtail millet to be further elucidated.

Test material
Test materials were selected from ordinary high-quality foxtail millet Jingu21 (Shanxi Academy of Agricultural Sciences Economic Crops Research Institute, China) and hybrid highyield foxtail millet Zhangza10 (Zhangjiakou City Academy of Agricultural Sciences, Hebei Province, China). rainfall, thus naturally simulating a drought environment. Crops from the previous season were not foxtail millet, and continuous cropping was therefore avoided.
A split-plot design was adopted for this study. Briefly, clear water was used as a blank control for the main plot, while HA at concentrations of 50, 100, 200, 300, and 400 mg L -1 (T1, T2,  T3, T4, and T5, respectively) was used for the secondary plot. Both water and HA were sprayed onto the leaves of the foxtail millet at the jointing and filling stages (800 L ha -1 ). The experiment was repeated three times in all 36 districts (10 m 2 district -1 ).

Plant quality evaluation
The foxtail millet was harvested on October 13, 2017, and October 1, 2018. The ear, stem, leaf, and root from three representative mature plants were baked at 105˚C for 30 min, followed by 80˚C for 24-48 h, then the dry weight was measured. Plant height, stem diameter, ear length, ear diameter, ear yardage, and 1000-grain weight were measured on three representative mature plants. Grain ears were manually threshed, shelled, comminuted (0.5 mm sieve), and dried. Protein content was determined using an MPA Fourier transform near-infrared spectrometer (Bruke, Germany). Total P and K content was determined via the molybdenum antimony colorimetric method and flame photometry. Presence of mineral elements, such as iron (Fe), Mn, Cu, Zn, Ca, and Magnesium (Mg), was determined using the diethylenetriaminepentaacetic acid (DTPA) extraction inductively-coupled plasma spectroscopy (ICP) method.

HA treatment evaluation based on the membership function method
The yield and quality of foxtail millet Jingu21 and Zhangza10 after treatment with different HA concentrations were comprehensively analyzed using the membership function method in fuzzy mathematics [13]. The membership function formula is: U ðX i Þ ¼ ðXi À X min Þ=ðX max À X min Þ ðindicator trait positively correlates with nutrient uptake and yieldÞ U ðX i Þ ¼ 1 À ðXi À X min Þ=ðXmax À X min Þ ðindicator trait negatively correlates with nutrient uptake and yieldÞ where U (X i ) is the membership function value; X i is the measured value of an index at each processing level; X max and X min are the maximum and minimum values within an indicator at all processing levels, respectively.

Grouping for molecular analysis
Leaves from 3-5 leaf potted foxtail millet seedlings were sequenced to ensure accuracy of the test. Sequencing samples (Jingu21) were divided into three groups; normal control group (CK), drought treatment group (D), and drought + HA treatment group (DHA). Plants in the CK group were cultured under normal conditions. Culture conditions for plants in groups D and DHA were in accordance with field trials. The optimal HA concentration identified during field trials was used as the default concentration for HA treatment (100 mg L -1 ). Seedlings in groups D and DHA were harvested and analyzed after five days of drought.

Transcriptome sequencing
Total RNA was extracted from samples from all groups using the TRIzol-based method (RNAiso Plus, TaKaRa, 9109). RNA quantification and purity were determined using diethyl phosphorocyanidate (DEPC) H 2 O as the blank control. Transcriptome sequencing was then performed based on Illumina high throughput sequencing. The design and detailed operation for sequencing were validated by Beijing Novogene Technology Co., Ltd (Project Number: P101SC18122767-01). Sequencing data were retained for subsequent analysis.

Quality control and preprocessing
Quality control and preprocessing were performed on the original sequencing data. A clean read was obtained by filtering joint and low-quality sequences. Quality control of the clean reads was performed using fastqc (version: 0.11.5, http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/) [14]. TopHat software (version: 2.1.0, http://ccb.jhu.edu/software/tophat/) [15] was used to locate clean reads on the reference genome Setaria_italica_v2.0 (version: Ensembl Plants) [16] for foxtail millet. Finally, the featureCounts tool (version: 1.6.0, http:// subread.sourceforge.net/) [17] was used to annotate the samples with a Gens genome annotation file (Ensembl Plants), and read information for each gene alignment was obtained.

DEGs analysis
According to the RNA read count data, the read count was pretreated using the TMM (trimmed mean of M values) normalization method within the edgeR package in R software [18,19]. DEGs between D vs. CK and D vs. DHA were then revealed using the quasi-likelihood (QL) F-test in the edgeR package. Results were visualized via a heat map using pheatmap software (https://cran.r-project.org/web/packages/pheatmap). DEGs in both the D vs. CK and D vs. DHA groups were then visualized using a Venn diagram generated using VENNY software (http://bioinfogp.cnb.csic.es/tools/venny/index.html) [20].

PPI network construction and module analysis
Protein interaction information was obtained according to the search tool for the retrieval of interacting genes (STING) database (version: 11.0, http://www.string-db.org/) [21], and PPI pairs among DEGs between groups D and DHA were predicted with median confidence (score) = 0.4. Next, a PPI network was constructed using Cytoscape software (version: 3.7.0, http://www.cytoscape.org/) [22], and molecular complex detection (MCODE, Version1.5.1 http://apps.cytoscape.org/apps/MCODE) [23], a plug-in of Cytoscape software, was used to screen significantly enriched modules from the PPI network with a module score � 4.

Real-time quantitative PCR
RT-qPCR was used to verify DEGs obtained in this study. Expression of SETIT_009509mg, SETIT_021707mg, SETIT_016840mg, SETIT_015030mg, SETIT_004913mg, and SETIT_ 016654mg was assessed by qPCR. Briefly, total RNA was extracted from samples from each group (CK1, CK2, CK3, D1, D2, D3, DHA1, DHA2, and DHA3) using TRIzol reagent (RNAiso Plus, TaKaRa, 9109) and quantified. Reverse transcription was performed using 5×primeScript RT Master Mix (Perfect Real Time, TAKARA, RR036A). Actin was used as a reference. Primers are listed in S1 Table. Reaction conditions were as follows: 50˚C for 3 min, 95˚C for 3 min, 40 cycles at 95˚C for 10 s, and 60˚C for 30 s. A fluorescence signal was recorded at the end of each cycle, and an amplification curve was generated. Relative expression of candidate genes was calculated using the 2 -ΔΔCT method [26].

Statistical analysis
Statistical analysis was carried out using SPSS 16.0 (Inc., Chicago, IL, USA) and GraphPad Prism 5 software (GraphPad Software, San Diego, CA). Results from field trials were processed and plotted using Microsoft Excel 2010, and data are presented as mean ± standard error (SE). RT-qPCR results are presented as mean ± standard deviation (SD). p < 0.05 and p < 0.01 were the screening criteria for significant and extremely significant differences, respectively.

Effect of HA on growth and yield characteristics of foxtail millet
As the HA concentration increased, the stem diameter of Zhangza10 initially increased then decreased, and this was true for samples collected in both 2017 and 2018. The effects of HA on plant height and stem diameter of Jingu21, and plant height of Zhangza10, were not significant (all p > 0.05) (Fig 1). In 2017, the underground dry weight of Jingu21 and Zhangza10 initially

PLOS ONE
increased then decreased with increasing HA concentrations (S2 Table). Jingu21 reached its maximum root-to-shoot ratio with T3 in 2017 and 2018, while Zhangza10 reached its maximum root-to-shoot ratio with T2 in both years. These results show that HA at a concentration of 100-200 mg L -1 significantly improved the dry weight and root-shoot ratio of foxtail millet. Moreover, there was an initial increase in ear yardage, ear weight, and yield in both Jingu21 and Zhangza10 as HA concentration increased, but they then decreased (Table 1). Jingu21 had the largest ear diameter and yield with T3 in 2017 and with T2 in 2018, while Zhangza10 had the largest ear diameter and yield with T3 in 2017. In 2018, the maximum ear diameter of Zhangza10 was recorded with T2, while the maximum yield occurred with T3, and this treatment also resulted in the maximum ear length. The yield of Jingu21 increased by 16.96% (twoyear average) with T3 and T2, compared to CK. Meanwhile, Zhangza10 yield increased by 14.48% (two-year average) with T3, compared to CK. These results show that HA at a concentration of 100-200 mg L -1 causes a significant increase in foxtail millet yield, and has a positive effect on dry weight and root-shoot ratio. The correlation coefficients between growth indicators and yield traits are listed in Table 2. There was a strong negative correlation between yield, plant height and aboveground dry weight (p < 0.01), and a positive correlation between yield and root-shoot ratio (p < 0.01). Thousand kernel weight showed a strong positive correlation with stem diameter and underground dry weight (p < 0.01).

Seed quality and membership function evaluation
As HA concentration increased, the protein content of Jingu21 and Zhangza10 grains initially increased and then decreased, in both 2017 and 2018 (Fig 2). Results of the membership function and comprehensive evaluation are shown in Table 3, and reveal an initial increase and subsequent decrease in nutrient uptake and yield in foxtail millet treated with HA. The order of action concentration for Jingu21 was T3 > T2 > T1 > T4 > T5, and for Zhangza10 was T2 > T3 > T1 > T4 > T5. The average membership function with T2 was 0.65 for Jingu21 and 0.80 for Zhangza10, and with T3 was 0.81 for Jingu21 and 0.69 for Zhangza10. The comprehensive analysis showed that the order of HA concentrations that improved foxtail millet yield and quality was T3 > T2 > T1 > T4 > T5.

DEG analysis based on mRNA sequencing data
Sequencing data identified 453 genes that were up-regulated and 645 that were down-regulated between groups D and CK. Meanwhile, a total of 272 up-regulated and 137 down-regulated genes were identified between groups D and DHA. A heat map for the union of the two groups is shown in Fig 3A, and indicates the relative expression values of DEGs among samples in different groups. Furthermore, a Venn plot analysis was performed on all DEGs ( Fig  3B), and revealed 79 common DEGs in both D vs. CK and D vs. DHA. Meanwhile, 1019 and 330 DEGs were identified for D vs. CK and D vs. DHA, respectively. Three hundred and thirty DEGs between groups D and DHA were selected for further investigation into the molecular mechanisms that underlie the effect of HA on dry matter and nutrient accumulation in foxtail millet.

PPI network and module investigation
A PPI network was constructed based on DEGs between groups D and DHA (Fig 4A) with a combined score of 0.4. Moreover, with score � 0.4, two modules, A (score = 4.8, six nodes, twelve interactions) and B (score = 4, four nodes, six interactions), were found to form the current PPI network (Fig 4B). Characteristic genes with top ten degrees were within the PPI network (such as SETIT_026527mg), module A (such as SETIT_001087MG), and module B (such as SETIT_025844mg).

Gene expression analysis
Expression of human SETIT_009509mg, SETIT_021707mg, SETIT_016840mg, SETIT_ 015030mg, SETIT_004913mg, and SETIT_016654mg genes was analyzed by RT-qPCR. Results show that SETIT_016654mg was significantly up-regulated in group D compared with the CK group (p < 0.05). Moreover, SETIT_021707mg, SETIT_016840mg, and SETIT_ 015030mg were significantly up-regulated, while SETIT_004913mg and SETIT_016654mg were significantly down-regulated in the DHA group compared with group D (all p < 0.05) PLOS ONE (Fig 6). Finally, SETIT-009509mg was significantly up-regulated in the DHA group compared to groups CK and D (all p < 0.05).

Discussion
Organic HA contains macronutrients and micronutrients, growth-promoting substances, vitamins, and beneficial microorganisms, all of which are important for maximizing plant yield. Previous research showed that HA treatment resulted in an increase in grain yield in wheat PLOS ONE [27]. For the two varieties of millet in the present study, there was no difference in P content between CK and the T4 treatment, whereas other concentrations of HA resulted in an increase in P content. This increase in phosphorus after HA treatment could be due its forming a complex with iron (Fe) [28]. Some trace elements, including Fe, have low solubility in soil at higher pH values, which can result in Fe deficiency in plants. The addition of HA can reduce the soil pH [29], which in turn would not only promote the growth of crops, but could also increase nutrient absorption from the soil [27]. Humic acid can either activate or inhibit enzyme activity by promoting the absorption of underground soil nutrients by the plant, thus altering cell

PLOS ONE
membrane permeability, which in turn leads to protein synthesis and stimulates growth, ultimately increasing the yield.
In the current study, HA was found to alleviate drought stress in foxtail millet. More specifically, treatment with 100-200 mg L -1 HA had a positive effect on underground dry weight and root-shoot ratio., and also increased the yield. There was a strong positive correlation between root-shoot ratio and yield (p < 0.01). An increase in underground dry weight is an indication of enhanced root growth, which could be a good way to alleviate moisture loss in the soil [30]. It has been shown that addition of HA to soil increases proliferation and branching of root hairs [31], which could in part explain the enhanced nutrient absorption and increased yield observed in the present study. HA treatment had little influence on ear length, ear diameter or ear yardage of millet, which could be because the HA was applied during the jointing and filling stages, when these ear characteristics are already defined [32].

PLOS ONE
Drought stress is one of the major abiotic stresses known to affect crop production worldwide. In order to understand the mechanism by which plants cope with water deficiency, it is necessary to study natural drought-tolerant plants and identify the molecular mechanisms that underlie their drought stress tolerance. The biological function of HA is commonly studied in terms of regulation at either the gene or molecular level [33]. In the present study, 330 DEGs were identified between groups D and DHA, and were found to be enriched mainly in metabolic pathways (SETIT_009509mg, SETIT_017919mg, and SETIT_021617mg), secondary metabolite biosynthesis (SETIT_014105mg, SETIT_017811mg, and SETIT_022131mg), and starch and sucrose metabolism (SETIT_009543mg, SETIT_034299mg, and SETIT_009543mg) [34,35], indicating that these pathways are dominant when the plant is under drought stress. The expression of six genes, which are involved in important metabolic pathways and with high degrees of interaction in a PPI network, were verified by RTqPCR. Of these, SETIT_016654mg was significantly up-regulated in group D compared to groups CK and DHA (p < 0.0001). SETIT_016654mg encodes arginine decarboxylase 2 (ADC2), which has been shown to be induced by drought stress in Arabidopsis thaliana [36]. Under high permeability conditions, such as after treatment with mannitol, ADC activity was found to be increased in the leaves and roots of wheat. Furthermore, putrescine is known to play an important role in salt tolerance in plants [37], and levels of putrescine and spermine were increased in the leaves and roots of plants treated with mannitol,.

Conclusions
In conclusion, HA at a concentration of 100 mg L -1 can increase the yield and protein and mineral content of foxtail millet grain under drought conditions. Moreover, SETIT_016654mg may play an important role in the effect of HA on foxtail millet by regulating metabolic pathways. This study identified the treatment dose of HA that should be used in millet under drought conditions, and lays the foundation for research into the molecular mechanisms that underlie the alleviating effects of HA on foxtail millet under drought conditions.