The Homeodomain-Leucine Zipper Genes Family Regulates the Jinggangmycin Mediated Immune Response of Oryza sativa to Nilaparvata lugens, and Laodelphax striatellus

The homeodomain-leucine zipper (HDZIP) is an important transcription factor family, instrumental not only in growth but in finetuning plant responses to environmental adversaries. Despite the plethora of literature available, the role of HDZIP genes under chewing and sucking insects remains elusive. Herein, we identified 40 OsHDZIP genes from the rice genome database. The evolutionary relationship, gene structure, conserved motifs, and chemical properties highlight the key aspects of OsHDZIP genes in rice. The OsHDZIP family is divided into a further four subfamilies (i.e., HDZIP I, HDZIP II, HDZIP III, and HDZIP IV). Moreover, the protein–protein interaction and Gene Ontology (GO) analysis showed that OsHDZIP genes regulate plant growth and response to various environmental stimuli. Various microRNA (miRNA) families targeted HDZIP III subfamily genes. The microarray data analysis showed that OsHDZIP was expressed in almost all tested tissues. Additionally, the differential expression patterns of the OsHDZIP genes were found under salinity stress and hormonal treatments, whereas under brown planthopper (BPH), striped stem borer (SSB), and rice leaf folder (RLF), only OsHDZIP3, OsHDZIP4, OsHDZIP40, OsHDZIP10, and OsHDZIP20 displayed expression. The qRT-PCR analysis further validated the expression of OsHDZIP20, OsHDZIP40, and OsHDZIP10 under BPH, small brown planthopper (SBPH) infestations, and jinggangmycin (JGM) spraying applications. Our results provide detailed knowledge of the OsHDZIP gene family resistance in rice plants and will facilitate the development of stress-resilient cultivars, particularly against chewing and sucking insect pests.


Introduction
The static nature of plants entails the frequent endurance of various environmental stresses. In response, plants have developed various mechanisms to adjust to the constantly changing environment [1]. These developmental processes are frequently regulated by numerous transcription factors (TFs) [2]. TFs spreads throughout the genome and can bond with certain functional cis-elements, facilitating the plant's response to various environmental stimuli. The homeodomain-leucine zipper (HDZIP) is a transcription factor family that plays a vital role in plant growth, developmental processes, and stress response [3]. The HDZIP is a class of homeobox proteins containing the homeodomain (HD) and leucine zipper (LZ) motifs [3,4]. These two motifs are the signatures of the HDZIP family and have been found in all eukaryotic species. However, their interaction with a single protein is only found in plants and, therefore, HDZIP in Plantae is different than in other organisms [5]. Based on their structure, sequence composition, functional characteristics, and phylogenetic relationship, the HDZIPs are divided into four subfamilies (i.e., HDZIP I, HDZIP II, HDZIP III, and HDZIP IV). Additionally, each subfamily has a unique function and forms a complex interactive network throughout the plant's developmental phases [6][7][8].

Phylogenetic Tree, Motif, and Digital Expression Analysis
To obtain detailed knowledge regarding the evolutionary relationship of the HDZIP gene family in rice with the developed genome species. Here, we investigated the phylogenetic relationships of the HDZIP gene family members of rice with model plants.
Firstly, the HDZIP amino acid sequences were downloaded from their corresponding genome database and were then aligned through ClustalW software (version 2.1) (http: //www.genome.jp/tools/clustalw/) (accessed on 22 February 2022) following the default parameters to examine the evolutionary relationships among the sequences and construct the maximum likelihood phylogenetic tree using MEGA (version 7.0) [35]. Furthermore, the conserved protein motifs of the HDZIP family of O. sativa were predicted using the MEME online server (Version 4.12.0) (http://meme-suite.org/) (accessed 10 March 2022) with the default settings. The details of the top 10 predicted motifs were obtained from the MEME suite. The conserved domains of the HDZIP gene family of O. sativa were predicted using the NCBI-CDD (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (accessed 10 March 2022). The conserved domain and motif distribution were drawn via Microsoft PowerPoint 365 software. Finally, for the visualization of the heatmap, TBtools (Version 1.098765) was used, and the transcriptomic data were retrieved from the rice genome database and GEO dataset platform from the NCBI [27,32].

Cis-Elements and Gene Ontology of the OsHDZIP Genes
To determine the cis-regulatory elements in each OsHDZIP gene, 1.5 kb of an upstream genomic DNA sequence with the starting codon (ATG) was obtained from the rice genome sequence database. Further, we used the plantCARE database (http://bioinformatics. psb.ugent.be/webtools/plantcare/html/) (accessed on 26 February 2022) to identify the cis-elements in the promoter regions of 40 OsHDZIP genes of O. sativa. Furthermore, for the GO analysis, the OsHDZIP protein sequences were downloaded from the iTAK-Plant Transcription Factor and Protein Kinase Identifier and Classifier (http://itak.feilab.net/cgibin/itak/index.cgi) (accessed on 5 March 2022). The obtained OsHDZIP protein sequences were implied in the "CELLO2GO" online server to determine the predicted functions, such as molecular functions, biological processes, and cellular components, and finally, the GO classifications were recovered using Microsoft Excel 365 software.

Interactive Protein Analysis of the OsHDZIP Genes
The online server String (https://string-db.org) (accessed on 8 February 2022) was used for the interactive protein network analysis (accessed on 8 March 2022), using the O. sativa OsHDZIP2 protein as a reference following the default advanced settings [36]. Furthermore, pathway enrichment analysis was carried out by searching for OsHDZIP genes in the rice genome database's online pathway enrichment tool [32].

Prediction of Putative MicroRNAs Targeting OsHDZIP Genes
To predict putative miRNA target sites in the OsHDZIP genes in the rice plants under BPH and SBP infestations and JGM spraying applications (briefly described in Section 2), the sequences of rice miRNAs were downloaded from the rice genome database. Moreover, the OsHDZIP CDS sequences were submitted to the online psRNA Target (Server18) with the default parameters for predicting potential miRNAs in OsHDZIP genes. The interactive Bioengineering 2022, 9, 398 4 of 25 network between the predicted miRNAs and OsHDZIP targeted genes were constructed and visualized using Cytoscape software (version 3.919) by the Institute for Systems Biology (Seattle, Washington 98103, USA) following the same procedure by Rizwan et al. [37].

Insect Rearing and Chemical and Stress Treatment
The insects (i.e., BPH and SBPH) and the rice variety used in the study, Ninjing4, were initially obtained from the China National Rice Research Institute insect repository (Hangzhou, China). The Ninjing4 rice variety, known to be resistanceless to these insect infestations, was used for insect rearing. Initially, the BPH and SBPH colonies were reared on the rice seedlings in cement tanks covered with fine mesh outdoors (i.e., natural conditions) for six months (i.e., April to October) and overwintered in lab-controlled conditions. First, we soaked the seeds for 24 h in a water dip plastic tray with a standard size of one-quarter (60 cm H_100 cm W_200 cm L) in standard conditions of 26 ± 2 • C with the 16 h L:8 h D in relative humidity of 80 ± 10% in the ecological laboratory of Yangzhou University. The germinated seeds were transferred to cement tanks covered with fine mesh in an outdoor natural environment and were grown until the six-leaf seedling stage. Secondly, the seedlings were then transferred into plastic pots (dimensions = R 1 4 16 cm). Finally, the stress treatments proceeded at the tillering stage (40 ± 2, 40 ± 4, and 40 ± 8 days).
The JGM technical grade of 61.7% used in this study was obtained from the Qianjiang Biochemistry Co., Ltd. (Haining, Zhejiang, China). Following the protocol of in a previous report by Ahmad et al. (2022), the two hundred parts per million (PPM) solution was Tween 20 obtained from the Sinopsin Group Chemical Reagent Company (Shanghai, China). The fungicide was then sprayed on the rice seedlings, following the procedure in a previous study [27,38].

Expression Profiling of HDZIP Genes in Oryza sativa
Forty day-old (40 ± 2) rice plants were exposed to BPH and SBPH stress, and samples were collected at 2, 4, and 8 days after infestation. Similarly, JGM was sprayed on the rice plants, and samples were taken 2, 4, and 8 days after treatment. The samples were then stored at −80 • C degrees until further experiments. After that, the total RNA was extracted from the samples using kits (Vazyme, Nanjing, China). First, the DNA was removed using DNase I, the concentration and purity were measured with a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Rockford, IL, USA), and the integrity was checked using 1.5% agarose gel electrophoresis. Finally, the resulting cDNA was used as a template for qPCR (quantitative real-time PCR) analysis using SYBR Green real-time PCR master mix (Vazyme, Nanjing, China). The qPCR assays were performed in triplicate using a real-time PCR system (Bio-Rad, Hercules, CA, USA) following the manufacturer's protocol [27].
Furthermore, the 2 µL aliquots of cDNA were amplified by qPCR in 20 µL reaction volumes using the SYBR Premix Ex Taq TM II (TaKaRa, Dalian, China). The cDNAs were amplified at 95 • C for 2 min, followed by 35 cycles of 10 s at 95 • C, then for 30 s, and at 72 • C for 30 s, with a final extension step of 72 • C for 10 min in a CFX96 real-time PCR system (Bio-Rad Co., Ltd., CA, USA). The mRNA amounts of all genes were separately quantified with the stable expression of the constitutive reference gene, actin. The specific primers are listed in (Table S1). After amplification, the target gene cycle threshold (Ct) values were normalized to the reference gene by the 2 −∆∆CT method [39]. The data's mean values of three biologically independent replicates were used for the final graphs [27].

Statistical Analysis
The data presented in this paper were analyzed using the SPSS software (version 25.0, SPSS Inc., Chicago, IL, USA) for statistical analysis (ANOVA), statistical significance, and a 95% confidence interval (p ≤ 0.05). The data were analyzed and are expressed as the mean ± standard deviation (SD) of three biologically independent replicates in all measured parameters, and finally, GraphPad Prism (Version 8.0.2) (GraphPad Software, Inc., LA Jolla, CA, USA) was used for graphical representation [27].

Identification and Sequence Analysis of OsHDZIP Genes in Oryza sativa
We retrieved 40 rice OsHDZIP transcription factors from the rice genome database. All genes had a nomenclature of OsHDZ1 to OsHDZ40 (Table 1). Among the 40 OsHDZIP genes family members, 35 OsHDZIP genes resided in the nucleus, whereas OsHDZIP3, OsHDZIP7, OsHDZIP13, and OsHDZIP40 were in the plasma membrane, while 2 genes, OsHDZIP28 and OsHDZIP37, resided in the cytoplasm. Other features of Oryza sativa OsHDZIP protein were identified, such as locus ID, subfamilies, chromosomal coordinates, molecular weight, chemical properties, and isoelectric point (PI), and tabulated.

The OsHDZIP Genes Conservative Domain Analysis
The HDZIP gene family consisted of two functional domains: homeodomain (HD) and leucine zipper (LZ). Based on their sequence conservation and functional properties, the HDZIP gene family was divided into four subfamilies (i.e., HDZIP I, HDZIP II, HDZIP III, and HDZIP IV). The HDZIP I subfamily contains the HD and LZ domains among the four subfamilies; HDZIP II contains a similar HD and LZ domain with an additional CPSCE motif; HDZIP III and HDZIP IV contain HD and LZ with an additional START and SAD domain. Only HDZIP III possesses the highly conserved MEKHLA domain ( Figure 1). The domain distribution and structure of the OsHDZIP genes in Oryza sativa are shown in Table S3.

The OsHDZIP Genes Conservative Domain Analysis
The HDZIP gene family consisted of two functional domains: homeodomain (HD) and leucine zipper (LZ). Based on their sequence conservation and functional properties, the HDZIP gene family was divided into four subfamilies (i.e., HDZIP I, HDZIP II, HDZIP III, and HDZIP IV). The HDZIP I subfamily contains the HD and LZ domains among the four subfamilies; HDZIP II contains a similar HD and LZ domain with an additional CPSCE motif; HDZIP III and HDZIP IV contain HD and LZ with an additional START and SAD domain. Only HDZIP III possesses the highly conserved MEKHLA domain (Figure 1). The domain distribution and structure of the OsHDZIP genes in Oryza sativa are shown in Table S3.

Evolutionary Relationship of Oryza sativa HDZIP Genes
After alignment, a maximum likelihood phylogenetic tree was constructed to gain insightful knowledge regarding the evolutionary relationship of the homeodomain-leucine zipper in Arabidopsis thaliana [30], Cucumis sativus [8], and O. sativa [32] (Figure 2). The protein sequence of the rice (40 OsHDZIP), A. thaliana (47 AtHDZIP) [40], and cucumber (40 CsHDZIP) [8] were downloaded from their respective databases and through the maximum likelihood method, the phylogenetic tree was constructed. Following the same procedure as our previously published article, first, all sequences were aligned using ClustalX software with default parameters, then the phylogenetic tree was constructed using MEGA6 software, and the final tree was programmed using Interactive Tree Of Life (iTOL) (version 5) (accessed on 20 March 2022) [41]. Additionally, the OsHDZIP proteins were clustered into four subfamilies based on their phylogenetic relationships: HDZIPI, HDZIPII, HDZIPIII, and HDZIPIV, and the number of OsHDZIP proteins was measured in each subfamily. Subfamily I accounted for 13 proteins, followed by subgroup II with 12 proteins, subfamily IV with 11 proteins, and subfamily III with the least number proteins at 4.
(iTOL) (version 5) (accessed on 20 March 2022) [41]. Additionally, the OsHDZIP proteins were clustered into four subfamilies based on their phylogenetic relationships: HDZIPI, HDZIPII, HDZIPIII, and HDZIPIV, and the number of OsHDZIP proteins was measured in each subfamily. Subfamily I accounted for 13 proteins, followed by subgroup II with 12 proteins, subfamily IV with 11 proteins, and subfamily III with the least number proteins at 4.

Figure 2.
Phylogenetic analysis of HDZIP: the phylogenetic tree was generated using the aminoacid sequences of selected HDZIPs via the maximum likelihood tree method. All Oryza sativa HDZIPs, Arabidopsis thaliana, and Cucumis sativus, with their counterparts, were classified into four subfamilies, and the final tree was displayed using the Interactive Tree Of Life (iTOL) (version 5). Phylogenetic analysis of HDZIP: the phylogenetic tree was generated using the amino-acid sequences of selected HDZIPs via the maximum likelihood tree method. All Oryza sativa HDZIPs, Arabidopsis thaliana, and Cucumis sativus, with their counterparts, were classified into four subfamilies, and the final tree was displayed using the Interactive Tree Of Life (iTOL) (version 5).

An Interactive Network of OsHDZIP Protein
The protein interaction analysis revealed various other proteins interacting with the OsHDZIP orthologous gene OsHDZIP2 ( Figure 3). The OsHDZIP2 protein from the homeodomain-leucine zipper gene family plays a crucial role in plant growth and stress response. For instance, the YUCCA pathway is the most important and well-characterized pathway that plants deploy to produce auxin, which is the essential hormone in plant development and stress response [42]. In addition, our reference protein, OsHDZIP2, was found to be highly interactive with rice LAZY1 proteins and had an essential role in auxin biosynthesis, the predicted functional partners of OsHDZIP proteins (Table S4).
OsHDZIP orthologous gene OsHDZIP2 ( Figure 3). The OsHDZIP2 protein from the homeodomain-leucine zipper gene family plays a crucial role in plant growth and stress response. For instance, the YUCCA pathway is the most important and well-characterized pathway that plants deploy to produce auxin, which is the essential hormone in plant development and stress response [42]. In addition, our reference protein, OsHDZIP2, was found to be highly interactive with rice LAZY1 proteins and had an essential role in auxin biosynthesis, the predicted functional partners of OsHDZIP proteins (Table S4).

Prediction of the Potential MicroRNAs Targeting OsHDZIP Genes
MicroRNAs are a class of small noncoding regulatory RNAs that control gene expression by directing target mRNA cleavage or translational repression [43]. In recent decades, several investigations have reported that the miRNAs regulate numerous stresses, plant development, and signal transduction. Therefore, to better understand the regulatory mechanism of miRNAs involved in the regulation of OsHDZIP genes, 56 putative miRNAs targeting four OsHDZIP genes were identified, as shown in the network illustration ( Figure 4).

Prediction of the Potential MicroRNAs Targeting OsHDZIP Genes
MicroRNAs are a class of small noncoding regulatory RNAs that control gene expression by directing target mRNA cleavage or translational repression [43]. In recent decades, several investigations have reported that the miRNAs regulate numerous stresses, plant development, and signal transduction. Therefore, to better understand the regulatory mechanism of miRNAs involved in the regulation of OsHDZIP genes, 56 putative miR-NAs targeting four OsHDZIP genes were identified, as shown in the network illustration ( Figure 4).

Figure 4.
Targeted miRNA sites of the OsHDZIP subfamily III genes in rice represents the functional network assembly of the HDZIP genes in Oryza sativa. The subfamily III (i.e., OsHDZIP9, OsHDZIP13, OsHDZIP37, and OsHDZIP40) were mapped to the co-expression database. This analysis revealed 56 unique miRNAs that exhibited various physical/functional interactions, and a network was then assembled based on these interactions.

Gene Ontology (GO) Analysis
The Gene Ontology (GO) enrichment pathway analysis showed various key functions of the OsHDZIP genes in Oryza sativa. Three functional predictions were analyzed into biological, molecular, and cellular processes [27]. According to the predicted biological processes, OsHDZIP genes play a crucial role in growth-related activities via hormonal and metabolic modulation ( Figure 6). Additionally, the response to external stimuli and the cellular analysis proved that 36 of the OsHDZIP genes were located in the nucleus, 4 genes resided in the plasma membrane, whereas OsHDZIP28 and OsHDZIP37 resided in the cytoplasm and may be involved in many cellular-based activities. At the same time, many molecular predictions regarding OsHDZIP genes indicated that they are involved in DNA-binding activities.

Gene Ontology (GO) Analysis
The Gene Ontology (GO) enrichment pathway analysis showed various key functions of the OsHDZIP genes in Oryza sativa. Three functional predictions were analyzed into biological, molecular, and cellular processes [27]. According to the predicted biological processes, OsHDZIP genes play a crucial role in growth-related activities via hormonal and metabolic modulation ( Figure 6). Additionally, the response to external stimuli and the cellular analysis proved that 36 of the OsHDZIP genes were located in the nucleus, 4 genes resided in the plasma membrane, whereas OsHDZIP28 and OsHDZIP37 resided in the cytoplasm and may be involved in many cellular-based activities. At the same time, many molecular predictions regarding OsHDZIP genes indicated that they are involved in DNA-binding activities.

Identified cis-Regulatory Elements in OsHDZIP Genes
The OsHDZIP in silico analysis revealed that the upstream region of HDZIP genes possesses various stress, hormonal, and growth receptive cis-regulatory elements. Here, we identified twenty-two cis-regulatory elements. Among these 22 cis-elements, 13 were related to stress and growth changes, and 9 were responsive to hormonal changes (Table  S2). Further, the hormonal responsive cis-regulatory elements, such as AuxRR-Core (auxin-responsive), the GCTCZ-motif (MeJA responsive), and the TCA-salicylic acid element, were responsive. Moreover, ABRE (abscisic acid responsive) was found in most of the OsHDZIP genes; meanwhile, the GARE-motif and P-box (a gibberellin-responsive ciselement) were identified in several genes upstream regions. Following a drought, the anaerobic induction-responsive cis-elements, MYB and ARE, and the light-responsive ciselements, such as the ZTCT-motif, G-box, and ACE, were also found in the majority of OsHDZIP genes. In addition, various stress-and growth-responsive cis-elements, such as CAT-box (involved in meristem expression), TC-rich repeats (stress-and defense-responsive cis-regulatory elements), MBS1 (regulates flavonoid biosynthesis gene expression), and LTR (low-temperature responsive), were observed in the promoter region of OsHDZIP genes. These clustering cis-regulatory elements in the promoter region of OsHDZIP genes imply their role in regulating gene expression during different growth stages and under environmental stimuli; the same results were reported in [8].

Gene Structure, and Motif Patterns of OsHDZIP Genes
We identified the exons-introns distribution using the CDS and genomic sequence of OsHDZIP genes. The OsHDZIP family genes contain multiple exons and have varied intron lengths (Figure 7). Among these, the OsHDZIP subfamily IV and subgroup III were observed with the highest number of intron and exon distribution in which all members (i.e., OsHDZIP9, OsHDZIP13, OsHDZIP37, and OsHDZIP40) had 20 exons and 17 introns; followed by subfamily IV in which the OsHDZIP3 and OsHDZIP17 had 13 exons and 10 introns; OsHDZIP16, OsHDZIP20, OsHDZIP32, and OsHDZIP39 had 11 exons and 8 introns in total. However, OsHDZIP2 had 10 exons and 9 introns, and with OsHDZIP24, 10 exons and 7 introns were observed. Furthermore, the third-highest number of exons and introns was recorded in subfamily II in which the OsHDZIP1, OsHDZIP15, OsHDZIP25, OsHDZIP27, and OsHDZIP30 had six exons and three introns, followed by OsHDZIP34 and OsHDZIP38 with five exons and two introns; meanwhile, OsHDZIP18, OsHDZIP19, and OsHDZIP30 had four exons and two introns; OsHDZIP21 had three exons and two introns; only OsHDZIP4 had the least number with two exons and a single intron. Additionally, the OsHDZIP subfamily I was counted as having the least number of introns and

Identified cis-Regulatory Elements in OsHDZIP Genes
The OsHDZIP in silico analysis revealed that the upstream region of HDZIP genes possesses various stress, hormonal, and growth receptive cis-regulatory elements. Here, we identified twenty-two cis-regulatory elements. Among these 22 cis-elements, 13 were related to stress and growth changes, and 9 were responsive to hormonal changes (Table  S2). Further, the hormonal responsive cis-regulatory elements, such as AuxRR-Core (auxinresponsive), the GCTCZ-motif (MeJA responsive), and the TCA-salicylic acid element, were responsive. Moreover, ABRE (abscisic acid responsive) was found in most of the OsHDZIP genes; meanwhile, the GARE-motif and P-box (a gibberellin-responsive ciselement) were identified in several gene's upstream regions. Following a drought, the anaerobic induction-responsive cis-elements, MYB and ARE, and the light-responsive ciselements, such as the ZTCT-motif, G-box, and ACE, were also found in the majority of OsHDZIP genes. In addition, various stress-and growth-responsive cis-elements, such as CAT-box (involved in meristem expression), TC-rich repeats (stress-and defense-responsive cis-regulatory elements), MBS1 (regulates flavonoid biosynthesis gene expression), and LTR (low-temperature responsive), were observed in the promoter region of OsHDZIP genes. These clustering cis-regulatory elements in the promoter region of OsHDZIP genes imply their role in regulating gene expression during different growth stages and under environmental stimuli; the same results were reported in [8].

Gene Structure, and Motif Patterns of OsHDZIP Genes
We identified the exons-introns distribution using the CDS and genomic sequence of OsHDZIP genes. The OsHDZIP family genes contain multiple exons and have varied intron lengths (Figure 7). Among these, the OsHDZIP subfamily IV and subgroup III were observed with the highest number of intron and exon distribution in which all members (i.e., OsHDZIP9, OsHDZIP13, OsHDZIP37, and OsHDZIP40) had 20 exons and 17 introns; followed by subfamily IV in which the OsHDZIP3 and OsHDZIP17 had 13 exons and 10 introns; OsHDZIP16, OsHDZIP20, OsHDZIP32, and OsHDZIP39 had 11 exons and 8 introns in total. However, OsHDZIP2 had 10 exons and 9 introns, and with OsHDZIP24, 10 exons and 7 introns were observed. Furthermore, the third-highest number of exons and introns was recorded in subfamily II, in which the OsHDZIP1, OsHDZIP15, OsHDZIP25, OsHDZIP27, and OsHDZIP30 had six exons and three introns, followed by OsHDZIP34 and OsHDZIP38 with five exons and two introns; meanwhile, OsHDZIP18, OsHDZIP19, and OsHDZIP30 had four exons and two introns; OsHDZIP21 had three exons and two introns; only OsHDZIP4 had the least number with two exons and a single intron. Additionally, the OsHDZIP subfamily I was counted as having the least number of introns and exons distributed in which OsHDZIP11, OsHDZIP25, OsHDZIP26, OsHDZIP28, OsHDZIP29, and OsHDZIP31 accounted for five exons and two introns; followed by OsHDZIP6, OsHDZIP10, OsHDZIP14, OsHDZIP33, OsHDZIP35, and OsHDZIP36 with four exons and a single intron; a single gene OsHDZIP8 was observed with six exons and three introns.  A total of 40 conserved motifs were discovered using the MEME online server (Version 5.4.1) (accessed on 10 April 2022) [44], and they were found to be appropriate for explaining the HDZIPs gene's structure (Figure 8)  A total of 40 conserved motifs were discovered using the MEME online server (Version 5.4.1) (accessed on 10 April 2022) [44], and they were found to be appropriate for explaining the HDZIPs gene's structure (Figure 8). Among the 40 OsHDZIP genes, the OsHDZIP7, OsHDZIP8, OsHDZIP9, OsHDZIP16, OsHDZIP17, OsHDZIP22, OsHDZIP23, and OsHDZIP32 were counted as having the highest number with 10 motifs; followed by OsHDZIP2 with 8 motifs and OsHDZIP3 with 7 motifs. Furthermore, OsHDZIP9, OsHDZIP13, OsHDZIP24, OsHDZIP34, OsHDZIP37, and OsHDZIP40 had five motifs. However, the OsHDZIP1, OsHDZIP4, OsHDZIP5, OsHDZIP12, OsHDZIP15, OsHDZIP18, OsHDZIP19, OsHDZIP21, OsHDZIP27, OsHDZIP30, and OsHDZIP38 in addition to all subfamily I members (i.e., OsHDZIP6, OsHDZIP8, OsHDZIP10, OsHDZIP11, OsHDZIP14, OsHDZIP25, OsHDZIP26, OsHDZIP28, OsHDZIP29, OsHDZIP31, OsHDZIP33, OsHDZIP35, and OsHDZIP36) counted has having three motifs in total. Bioengineering 2022, 9, x FOR PEER REVIEW 13 of 25

Microarray Expression Analysis of HDZIP Genes in Developmental Stages
We examined the different developmental stages and tissue-specific expressions to study the biological roles of OsHDZIP genes in plant growth and development based on a set of microarray data obtained from the RiceXPro expression database (version 3.0) (accessed on 13 April 2022) [45]. The microarray expression data analysis of the rice OsHDZIP genes family is presented as a heatmap, with blue to red colors reflecting the expression pattern. In twelve tissues (i.e., leaf blade, leaf sheath, roots, stem, inflorescences, anther, pistil, lama, plea, ovary, embryo, and endosperm), the OsHDZIP gene family members showed moderate to high expressions, respectively (Figure 9). Among the twelve tissues, OsHDZIP27 showed dominant expression in leaf blade, leaf sheath, and roots, followed by OsHDZIP18, which showed high transcription in roots. Furthermore, OsHDZIP1, OsHDZIP3, OsHDZIP5, OsHDZIP13, OsHDZIP33, OsHDZIP34, OsHDZIP35, OsHDZIP39, and OsHDZIP40 showed moderate expression in leaf blade, root, inflorescences, and endosperms. All developmental stages were observed with no or extremely low transcripts, particularly in embryos and endosperm. On the contrary, the leaf blade, root, inflorescence, and anther had high transcription levels. Additionally, the developmental stages, including pistil, lama, and plea, revealed the response of many OsHDZIP gene expressions. This expression analysis revealed the essential role of the OsHDZIP gene family in developmental stages.

Microarray Expression Analysis of HDZIP Genes in Developmental Stages
We examined the different developmental stages and tissue-specific expressions to study the biological roles of OsHDZIP genes in plant growth and development based on a set of microarray data obtained from the RiceXPro expression database (version 3.0) (accessed on 13 April 2022) [45]. The microarray expression data analysis of the rice OsHDZIP genes family is presented as a heatmap, with blue to red colors reflecting the expression pattern. In twelve tissues (i.e., leaf blade, leaf sheath, roots, stem, inflorescences, anther, pistil, lama, plea, ovary, embryo, and endosperm), the OsHDZIP gene family members showed moderate to high expressions, respectively (Figure 9). Among the twelve tissues, OsHDZIP27 showed dominant expression in leaf blade, leaf sheath, and roots, followed by OsHDZIP18, which showed high transcription in roots. Furthermore, OsHDZIP1, Os-HDZIP3, OsHDZIP5, OsHDZIP13, OsHDZIP33, OsHDZIP34, OsHDZIP35, OsHDZIP39, and OsHDZIP40 showed moderate expression in leaf blade, root, inflorescences, and endosperms. All developmental stages were observed with no or extremely low transcripts, particularly in embryos and endosperm. On the contrary, the leaf blade, root, inflorescence, and anther had high transcription levels. Additionally, the developmental stages, including pistil, lama, and plea, revealed the response of many OsHDZIP gene expressions. This expression analysis revealed the essential role of the OsHDZIP gene family in developmental stages.

Expression Analysis of OsHDZIP Genes under Salinity
Salinity is an important stress that hinders plant growth and yield. The injurious effects of salinity can be noted over the whole plant. To obtain insightful knowledge regarding HDZIP genes' responsiveness to high salinity, the transcriptomic expression data were obtained from the publicly available in the NCBI (GSE102152) (accessed on 17 April 2022) [46]. In the heatmap, the dark orange color on the scale bar represents high expression, whereas the light orange color is moderate, and the blue color genes have no expressions ( Figure 10A). Among

Expression Analysis of OsHDZIP Genes under BPH, SSB, and SSB_BPH
The transcriptomic expression data were obtained from the NCBI (GSE167872) (accessed on 11 May 2022) [47]. The analysis provided insightful predictions of OsHDZIP genes involved in rice plant defense against brown planthopper (BPH), rice striped stem borer (SSB), Chilo suppressalis, and combined SSB_BPH stresses ( Figure 10B). However, in response to BPH, SSB, and the combined stress of SSB_BPH, only two genes, OsHDZIP4 and OsHDZIP10, had dominant expressions. The expression analysis revealed the role of the OsHDZIP gene family in plants' defense against pest infestations.

Expression Analysis of OsHDZIP Genes under Cnaphalocrocis medinalis
The transcriptomic expression data were taken from the NCBI (GSE159259) (accessed on 15 May 2022) [48]. The analysis provided insightful predictions of the OsHDZIP gene's role in rice plant defense against rice leaf folder (RLF) Cnaphalocrocis medinalis Guenée (Lepidoptera: Crambidae) ( Figure 10C). Further, the response of the HDZIP gene family was moderate; among the 40 OsHDZIP genes, only OsHDZIP4 and OsHDZIP10 were highly expressed at all time points (i.e., 0, 6, 12, and 24 h), and the rest of the genes did not show transcription. This expression analysis suggests that the HDZIP gene family plays a role in the rice defense system against biotic stress.

Expression Analysis of OsHDZIP Genes under Brassinosteroids and Jasmonic Acid
The expression data of the OsHDZIP genes family under Jasmonic acid (JA) and brassinosteroids (BRs) were obtained from the RiceXPro expression database (version 3.0) (accessed on 25 May 2022) [45]. The microarray expression data analysis of the rice HDZIP genes family is presented as a heatmap, with blue to red colors reflecting the expression

Expression Analysis of OsHDZIP Genes under BPH, SSB, and SSB_BPH
The transcriptomic expression data were obtained from the NCBI (GSE167872) (accessed on 11 May 2022) [47]. The analysis provided insightful predictions of OsHDZIP genes involved in rice plant defense against brown planthopper (BPH), rice striped stem borer (SSB), Chilo suppressalis, and combined SSB_BPH stresses ( Figure 10B). However, in response to BPH, SSB, and the combined stress of SSB_BPH, only two genes, OsHDZIP4 and OsHDZIP10, had dominant expressions. The expression analysis revealed the role of the OsHDZIP gene family in plants' defense against pest infestations.

Expression Analysis of OsHDZIP Genes under Cnaphalocrocis medinalis
The transcriptomic expression data were taken from the NCBI (GSE159259) (accessed on 15 May 2022) [48]. The analysis provided insightful predictions of the OsHDZIP gene's role in rice plant defense against rice leaf folder (RLF) Cnaphalocrocis medinalis Guenée (Lepidoptera: Crambidae) ( Figure 10C). Further, the response of the HDZIP gene family was moderate; among the 40 OsHDZIP genes, only OsHDZIP4 and OsHDZIP10 were highly expressed at all time points (i.e., 0, 6, 12, and 24 h), and the rest of the genes did not show transcription. This expression analysis suggests that the HDZIP gene family plays a role in the rice defense system against biotic stress.

Differential Expression of OsHDZIP Genes in Response to Nilaparvata lugens, Laodelphax striatellus Infestations, and JGM Spraying
To further investigate the response of OsHDZIP genes under BPH and SBPH infestations and botanical fungicide JGM applications, herein, we performed quantitative real-time polymerase chain reaction (qRT-PCR) to analyze the expression patterns at three time points over 2 day (2D), 4D, and 8D long treatments during BPH and SBPH infestations and JGM spraying (Figures 13-15). The results revealed that the response of ten candidate genes were expressed under all stress conditions. The BPH is a severe rice pest and causes a huge loss in annual rice production. The expression analysis revealed the important aspect of the O. sativa HDZIP transcription factor through moderate and high expressions. Among the eight candidate genes (two genes from each subfamily) at three different time points (Figure 13), moderate to high expression patterns were found, particularly for OsHDZIP20, which had the highest expression of six-fold over 2 days of infestation followed by 4 days with four-fold expression and two-fold after 8 days of BPH infestation. However, OsHDZIP03, OsHDZIP28, and OsHDZIP40 were found with low transcription; in addition, OsHDZIP04, OsHDZIP10, OsHDZIP15, and OsHDZIP37 were found with moderate expressions.

Differential Expression of OsHDZIP Genes in Response to Nilaparvata lugens, Laodelphax striatellus Infestations, and JGM Spraying
To further investigate the response of OsHDZIP genes under BPH and SBPH infestations and botanical fungicide JGM applications, herein, we performed quantitative real-time polymerase chain reaction (qRT-PCR) to analyze the expression patterns at three-time points over 2 days (2D), 4D, and 8D long treatments during BPH and SBPH infestations and JGM spraying (Figures 13-15). The results revealed that the response of ten candidate genes were expressed under all stress conditions. The BPH is a severe rice pest and causes a huge loss in annual rice production. The expression analysis revealed the important aspect of the O. sativa HDZIP transcription factor through moderate and high expressions. Among the eight candidate genes (two genes from each subfamily) at three different time points (Figure 13), moderate to high expression patterns were found, particularly for OsHDZIP20, which had the highest expression of six-fold over 2 days of infestation, followed by 4 days with four-fold expression and twofold after 8 days of BPH infestation. However, OsHDZIP03, OsHDZIP28, and OsHDZIP40 were found with low transcription; in addition, OsHDZIP04, OsHDZIP10, OsHDZIP15, and OsHDZIP37 were found with moderate expressions.
The SBPH is the second most important pest of rice after BPH, which causes a drastic loss to rice plants and their production. To validate the OsHDZIP transcription factor response through qRT-PCR under SBPH infestation, we used eight candidate genes at three time points (Figure 14). Among them, OsHDZIP03 displayed moderate expressions at 2 and 4 days of treatment and had the highest expression at 8 days of infestation up to 10-fold. This was followed by OsHDZIP40, with a six-fold expression at eight days. Furthermore, OsHDZIP04, OsHDZIP10, OsHDZIP15, OsHDZIP28, and OsHDZIP37 were found to have moderate expressions; however, OsHDZIP20 was found to have the lowest expression pattern.
JGM is a synthetic antibiotic applied to treat rice sheath blight disease, and it is also reported to enhance BPH fecundity; herein, we performed qRT-PCR analysis to validate the OsHDZIP gene family's response under JGM treatment ( Figure 15). The obtained results revealed that OsHDZIP genes, such as OsHDZIP10, OsHDZIP20, and OsHDZIP40, present upregulated expressions, whereas OsHDZIP04, OsHDZIP15, and OsHDZIP37 had moderate expressions and a single gene, and OsHDZIP03, OsHDZIP04, and OsHDZIP10 had negligible expression patterns. These results suggest that the OsHDZIP participated in rice immunity regulations against JGM spraying applications. The SBPH is the second most important pest of rice after BPH, which causes a drastic loss to rice plants and their production. To validate the OsHDZIP transcription factor response through qRT-PCR under SBPH infestation, we used eight candidate genes at three time points (Figure 14). Among them, OsHDZIP03 displayed moderate expressions at 2 and 4 days of treatment and had the highest expression at 8 days of infestation up to 10fold. This was followed by OsHDZIP40, with a six-fold expression at eight days. Furthermore, OsHDZIP04, OsHDZIP10, OsHDZIP15, OsHDZIP28, and OsHDZIP37 were found to have moderate expressions; however, OsHDZIP20 was found to have the lowest expression pattern. The SBPH is the second most important pest of rice after BPH, which causes a drastic loss to rice plants and their production. To validate the OsHDZIP transcription factor response through qRT-PCR under SBPH infestation, we used eight candidate genes at three time points (Figure 14). Among them, OsHDZIP03 displayed moderate expressions at 2 and 4 days of treatment and had the highest expression at 8 days of infestation up to 10fold. This was followed by OsHDZIP40, with a six-fold expression at eight days. Furthermore, OsHDZIP04, OsHDZIP10, OsHDZIP15, OsHDZIP28, and OsHDZIP37 were found to have moderate expressions; however, OsHDZIP20 was found to have the lowest expression pattern. reported to enhance BPH fecundity; herein, we performed qRT-PCR analysis to validate the OsHDZIP gene family's response under JGM treatment ( Figure 15). The obtained results revealed that OsHDZIP genes, such as OsHDZIP10, OsHDZIP20, and OsHDZIP40, present upregulated expressions, whereas OsHDZIP04, OsHDZIP15, and OsHDZIP37 had moderate expressions and a single gene, and OsHDZIP03, OsHDZIP04, and OsHDZIP10 had negligible expression patterns. These results suggest that the OsHDZIP participated in rice immunity regulations against JGM spraying applications. Figure 15. Represents the differential expression analysis of OsHDZIP genes (i.e., OsHDZIP3, OsHDZIP4, OsHDZIP10, OsHDZIP15, OsHDZIP20, OsHDZIP28, OsHDZIP37, and OsHDZIP40) under JGM spraying treatment in Oryza sativa.

Discussion
Plants counter various biotic and abiotic stresses during their life cycle, impairing their biochemical and physiological processes. The plant develops its mechanism to tackle these adverse conditions, such as activation of stress-related TFs and intensified metabolic activities. HDZIP TFs are distributed widely across the plant kingdom and have been recognized for their role in the growth and developmental activities and response to environmental stimuli [7,8,49].

OsHDZIP Genes Are Widely Distributed in the Rice Genome
Phylogenetic trees represent an established method for determining evolutionary changes and functional relationships [50]. The HDZIP protein family has been identified in many species, from mosses to higher plants, such as Ceratopteris richardii [51] and Physcomitrella patens [52], angiosperms, and gymnosperms [53]. In our study, we conducted a genome-wide survey to determine the phylogenetic relationships and investigate their potential role via qRT-PCR expression analysis; there was a total of 127 HDZIP proteins from rice (40 OsHDZIPs), Cucumis sativus (40 CsHZIPs) [8], and Arabidopsis thaliana (47 AtHDZIPs) [40]. These proteins were further divided into four subfamilies ( Figure  2). The obtained results aligned with previously reported studies on Arabidopsis thaliana and Cucumis sativus [8,40]. The HDZIP III subfamily was more highly conserved than the other subfamilies. On the other hand, HDZIP I, HDZIP II, and HDZIP IV varied in the numbers of different species [54][55][56]. The conserved number of motifs in OsHDZIP genes was analyzed, and the results stated that the HDZIP III and HDZIP IV subfamilies exhibited the highest number of motifs of the other two subfamilies (Figure 8

Discussion
Plants counter various biotic and abiotic stresses during their life cycle, impairing their biochemical and physiological processes. The plant develops its mechanism to tackle these adverse conditions, such as activation of stress-related TFs and intensified metabolic activities. HDZIP TFs are distributed widely across the plant kingdom and have been recognized for their role in the growth and developmental activities and response to environmental stimuli [7,8,49].

OsHDZIP Genes Are Widely Distributed in the Rice Genome
Phylogenetic trees represent an established method for determining evolutionary changes and functional relationships [50]. The HDZIP protein family has been identified in many species, from mosses to higher plants, such as Ceratopteris richardii [51] and Physcomitrella patens [52], angiosperms, and gymnosperms [53]. In our study, we conducted a genome-wide survey to determine the phylogenetic relationships and investigate their potential role via qRT-PCR expression analysis; there was a total of 127 HDZIP proteins from rice (40 OsHDZIPs), Cucumis sativus (40 CsHZIPs) [8], and Arabidopsis thaliana (47 AtHDZIPs) [40]. These proteins were further divided into four subfamilies (Figure 2). The obtained results aligned with previously reported studies on Arabidopsis thaliana and Cucumis sativus [8,40]. The HDZIP III subfamily was more highly conserved than the other subfamilies. On the other hand, HDZIP I, HDZIP II, and HDZIP IV varied in the numbers of different species [54][55][56]. The conserved number of motifs in OsHDZIP genes was analyzed, and the results stated that the HDZIP III and HDZIP IV subfamilies exhibited the highest number of motifs of the other two subfamilies (Figure 8). A motif's position in each subfamily was highly conserved, supporting their evolutionary classification into different subfamilies [40,55,57].

OsHDZIP Genes Have Tissue Specificity and Play an Integral Role in the Development of Oryza sativa
The HDZIP I TFs were documented for their involvement in plant developmental processes such as root growth and stem elongation, leaf morphology, flowering induction, and pollen hydration [8]. The OsHDZIP gene family members, including OsHDZIP18 and OsHDZIP27, displayed dominant expression in all tested organs of O. sativa in the developmental stages. Similar results were reported for the OsHDZIP18, OsHDZIP27 homologs CsHDZ02, and CsHDZ33, respectively. These results indicate the potential importance of OsHDZIP18, OsHDZIP27, and other subfamily members in regulating the growth and developmental activities of O. sativa [8].
Members of the HDZIP III subfamily were shown to be important for appropriate morphogenesis, embryo development, giving support during the commencement of lateral organ development, control of the shoot apical meristem, and water and nutrient transport throughout the plant's overall body [6,8,59,60]. Similarly, all of the genes in the HDZIP III subfamily showed transcription in practically all of the investigated tissues (Figure 9), suggesting that they may be important in sustaining rice's normal growth activities.
Members of the HDZIP IV subfamily have been implicated in developing root hairs, trichome production, anthocyanin biosynthesis, and flowering regulation in plants. [49,61]. OsHDZIP11 and OsHDZIP37 were expressed in all the tissues (Figure 9), suggesting that subfamily HDZIP IV plays a potential role in the developmental activity regulations of O. sativa.

OsHDZIP Genes Regulate Plant Response to Chewing Insects and other Abiotic Stresses
The homeodomain-leucine zipper transcription factors play a potential role in a plant's growth and development and are highly responsive to the pronounced effects of stresses including drought [62], cold [63], salinity [64], heat stress [65], heavy metal [66], flooding stress [67], and nutrient stress (iron deficiency) [68].
The BPH and SBPH are serious rice pests inflicting damage on a massive scale across Asia [19,21]. Their nymphs and adults cause direct damage by feeding on phloem sap from the tillering to the milking stages of rice and, in the process, transmitting viral pathogens, such as rice ragged stunt virus (RRSV) and rice grassy stunt virus (RGSV) [69,70]. Several studies have identified the potential role of the HDZIP gene family. Until now, no comprehensive study has been conducted to unfold the response of OsHDZIP genes to BPH and SBPH infestations. Herein, we performed qRT-PCR on eight candidate genes (two from each subfamily) to validate their expressions. The results revealed the potential role of OsHDZIP genes. Furthermore, among these eight candidate genes, OsHDZIP20 was observed with dominant expressions, and the rest of the genes (i.e., OsHDZIP04, Os-HDZIP10, OsHDZIP15, and OsHDZIP37) had moderate expressions, whereas OsHDZIP03, OsHDZIP28, and OsHDZIP40 had low transcription (Figures 13 and 14). Collectively, the response of these genes uncovered the crucial role of OsHDZIPs in rice plant immunity in response to insect pest infestations. However, a functional study is required to investigate the underlying mechanisms inducing the HDZIP gene family response against insect pests.
The antibiotic JGM was developed in recent decades in China, and it is usually applied two to three times in rice fields to treat rice sheath blight disease (Rhizoctonia solani) and fungal infections [25]. Moreover, JGM is reported to be an enhanced controlling agent of the sheath blight disease; however, JGM also has consequences because of its potential role in inducing BPH fecundity [25]. For instance, JGM was applied to the rice plants at the rate of 200 parts per million (ppm), and the results revealed that JGM increased the rice's resistance to rice sheath blight disease by disrupting the fungal cell wall and reduced sporulation [25]. In addition, JGM was also reported to enhance BPH fecundity [25]. In the current study, the OsHDZIP gene showed a differential expression pattern where the candidate genes from each subfamily, including OsHDZIP40, OsHDZIP28, and OsHDZIP10, were dominantly expressed, whereas the OsHDZIP3, OsHDZIP4, OsHDZIP15, and OsHDZIP37 were moderately expressed ( Figure 15). We speculate that these OsHDZIP genes are essential for inducing plant immunity to fungal pathogens in rice plants. Further studies are required to elucidate the underlying mechanisms of OsHDZIP genes boosting rice immune responses under various insect pest infestations.
The osa-miR166 family targeted four OsHDZIP genes, and miR166 has been reported to be involved in drought stress in maize [72], cowpea [73], soybean (Glycine max) seed development [74], peanut disease-resistance [75], and plant growth, development, and stress response in apple [77]. miR444 has been reported to be involved in the cadmium stress regulation of rice [78]. Previously, Anca et al. (2012) reported that the expression profiles of three miRNAs (i.e., osa-miR414, osa-miR408, and osa-miR164) targeting the OsABP, OsDBH and OsDSHCT genes in rice enhanced the miRNAs response against salinity stress [79]. Our results found that four OsHDZIP genes (i.e., OsHDZIP9, OsHDZIP13, OsHDZIP37, and OsHDZIP40) were targeted by these miRNAs providing another aspect of these miRNAs' potential role in regulating rice response under salinity ( Figure 4) and their target sites ( Figure 5). Collectively, these findings suggest that miRNAs may play vital roles in numerous growth and developmental processes and stress regulation processes by altering the transcriptional level of HDZIP genes in O. sativa. Among these 28 miRNAs families targeting the HDZIP gene family, several miRNAs predicted expression levels and functions have been identified, whereas various miRNAs targeting OsHDZIP genes that regulate biotic/abiotic stress responses and other important agronomic traits remain to be clarified.

Expression Analysis of HDZIP Gene Family under Hormonal Applications
JA is synthesized from linolenic acid through the action of several enzymes in plant chloroplast membranes, and the current evidence indicates that it induces resistance against necrotrophic pathogens and chewing herbivores [80,81]. In comparison, the microarray data showed that the expression of OsHDZP3, OsHDZP9, and OsHDZP28 was enhanced, suggesting their possible involvement in JA-mediated rice resistance against BPH ( Figure 11).
BRs are a group of polyhydroxylated steroidal phytohormones that are mandatory for plant development, growth, and productivity [82]. In addition to their significant involvement in growth-related activities, BRs are a key stress hormone [27]. Herein, the upregulated expressions of OsHDZIP15 and OsHDZIP20 suggest the essential role of these genes and their possible participation in the immunity regulation of the rice plant; however, functional studies are required to unfold the underlying mechanism of these crucial hormones (Figure 12).

Conclusions
Forty OsHDZIP family TFs were identified in the rice genome database and classified into four subfamilies based on their domain and structural properties. Furthermore, heatmap analysis revealed the expression of OsHDZIP TFs in distinct O. sativa tissues.
The expression of OsHDZIP genes under salinity and hormone stress indicated that they might play a role in modulating O. sativa resistance. Furthermore, differential expression patterns were found under JA and BR treatments, indicating that OsHDZIP genes may be critical in hormonal-mediated rice immunity against BPH and SBPH. On the other hand, stress control is a complicated system, and our results suggested the potential application of HDZIP biomarkers in developing stress-resilient rice lines. These results provide the missing role of HDZIP genes in regulating rice immunity against SBPH and particularly BPH.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/bioengineering9080398/s1, Table S1: Primer sequences for RT-qPCR, Table S2: Types and number of cis-acting regulatory elements analysis involved in growth, development, and stress and hormonal response; Table S3: Gene structure and domain distribution; Table S4: Predicted functional partners of OsHDZIP proteins; Table S5: A detailed list of miRNAs.
Author Contributions: S.A. and L.G. designed the research; S.A. and Y.C. conducted the experimental work; A.Z.S., H.Z., H.W. and C.X. contributed to the preparation of biological materials; S.A. performed bioinformatics analysis and wrote the main manuscript. All authors have read and agreed to the published version of the manuscript.