text:
Genome-Wide Characterization and Phylogenetic Analysis of GSK Genes in Maize and Elucidation of Their General Role in Interaction with BZR1

Glycogen synthase kinase-3 (GSK-3) is a nonreceptor serine/threonine protein kinase that is involved in diverse processes, including cell development, photomorphogenesis, biotic and abiotic stress responses, and hormone signaling. In contrast with the deeply researched GSK family in Arabidopsis and rice, maize GSKs’ common bioinformatic features and protein functions are poorly understood. In this study, we identified 11 GSK genes in the maize (Zea mays L.) genome via homologous alignment, which we named Zeama;GSKs (ZmGSKs). The results of ZmGSK protein sequences, conserved motifs, and gene structures showed high similarities with each other. The phylogenetic analyses showed that a total of 11 genes from maize were divided into four clades. Furthermore, semi-quantitative RT-PCR analysis of the GSKs genes showed that ZmGSK1, ZmGSK2, ZmGSK4, ZmGSK5, ZmGSK8, ZmGSK9, ZmGSK10, and ZmGSK11 were expressed in all tissues; ZmGSK3, ZmGSK6, and ZmGSK7 were expressed in a specific organization. In addition, GSK expression profiles under hormone treatments demonstrated that the ZmGSK genes were induced under BR conditions, except for ZmGSK2 and ZmGSK5. ZmGSK genes were regulated under ABA conditions, except for ZmGSK1 and ZmGSK8. Finally, using the yeast two-hybrid and BiFC assay, we determined that clads II (ZmGSK1, ZmGSK4, ZmGSK7, ZmGSK8, and ZmGSK11) could interact with ZmBZR1. The results suggest that clade II of ZmGSKs is important for BR signaling and that ZmGSK1 may play a dominant role in BR signaling as the counterpart to BIN2. This study provides a foundation for the further study of GSK3 functions and could be helpful in devising strategies for improving maize.

Introduction
Glycogen synthase kinase-3 (GSK-3) proteins constitute an important group of kinases. GSK-3 regulates fundamental biological pathways through inhibition and lies downstream of multiple cell signaling factors in eukaryotes [1][2][3]. GSK-3 kinases are defined by a conserved serine/threonine protein kinase catalytic domain (PF0069.24) [4]. GSK-3 was first reported in mammals as a kinase [5], and two isoforms have been reported in mammals: GSK-3α and GSK-3β. Both are involved in multiple physiological processes, including the regulation of the cell cycle, the induction of cell differentiation, and glycogen metabolism [6,7]. In contrast, GSK-3 kinase genes have been shown to have more diverse isoforms and functions in plants [1,8]. However, the isoform composition and functions of ZmGSK-3 proteins in maize are unclear. Notably, whether they are involved in hormone signaling is unknown.
BIN2 is a member of the GSK-3 kinase family that has been identified in many plants [9,10]. In Arabidopsis, BIN2 participates in the regulation of stress response by interacting with and phosphorylating RD26 [11]. TaGSK3, the wheat ortholog of Arabidopsis BIN2, induces round grains [12]. In addition, the knockout and overexpression of OsGSK1, an ortholog of Arabidopsis BIN2, showed that OsGAK1 plays important roles in the stress-response signaling pathway and reproductive organ development [13]. In Arabidopsis, a GSK-3 kinase mutant named dwf12 exhibits significant dwarfism and leaf curling [14]. GSK-3 knockouts exhibit morphological similarities to Brassinosteroid (BR)signaling mutants and display insensitivity to brassionazole (BRZ; an inhibitor of BR biosynthesis) [15]. Similar to the BR-signaling mutant, other GSK-3 mutant or overexpression models have indicated a negative role of GSK-3 in the BR signaling pathway, which led to its being recognized as a phytohormone.
Numerous studies have revealed that plant GSK-3 kinases generally function through interaction with transcription factors, especially BZR1 [16]. For instance, BIN2 mainly regulates downstream BR-related transcription networks by controlling the phosphorylation status of BZR1 to mediate the BR signaling pathway in Arabidopsis [17]. In the absence of BRs, BZR1/BES1 function can be inhibited due to phosphorylation by BIN2 [18,19]. In the presence of BRs, BZR1/BES1 dephosphorylates and accumulates in the nucleus to regulate target-gene expression because BIN2 kinase activity is inhibited [20]. BZR1 transcription factors play critical roles in BR signaling [21]. BZR1 regulates CBF1 and CBF2 expression to modulate Arabidopsis' freezing tolerance [22]. BZR1 can interact with other transcription factors to participate in BR-induced gene regulation. For example, AtIWS1 can interact with BES1 both in vitro and in vivo to promote plant-steroid-induced gene expression [23]. Another example is AtMYB30, an MYB transcription factor that cooperatively interacts with BES1 to promote BR-target-gene expression and is targeted by BES1 to amplify BR signaling [24]. In addition, BES1 can recruit other transcription factors to coordinate BR responses [25].
In contrast, little is known about whether the GSK-3 family could interact with BZR1 to participate in the BR signaling pathway in maize. Maize is the main cereal grain crop that contributes to solving the food security problem, and it is also known to be a social security concern for farmers [26][27][28][29]. Thus, understanding the GSK-3 family is significant. In this report, we aimed to investigate the role of ZmGSK in maize using a comprehensive analysis of phylogeny, gene structure, conserved motifs, cis elements, chromosomal location, and protein-protein interaction. Furthermore, we examined the expression pattern of ZmGSK genes in various developmental tissues, investigating their function using different hormones, and probing their interaction with BZR1.

Genome-Wide Identification of ZmGSKs in Maize Genome
After removing redundant sequences, 11 GSK genes were identified in the maize genome. The parameters used to describe the ZmGSK proteins are listed in Table 1. Putative ZmGSK proteins varied from 403 (ZmGSK7) to 1076 (ZmGSK3) amino acids (AAs) in length, with molecular weights (M W s) ranging from 45.17 to 120.83 KDa. Additionally, the theoretical pI of most ZmGSK proteins was higher than 7, except for those of ZmGSK6 (6.79) and ZmGSK3 (6.33).

Comparative and Phylogenetic Analyses of ZmGSKs
In order to precisely reveal the evolutionary relationships of the ZmGSK proteins, we performed phylogenetic analyses on 30 GSK proteins, including 11 maize, 10 Arabidopsis, and 9 rice GSKs, using the neighbor-joining method ( Figure 1). All proteins were categorized into four groups (I, II, III, and IV). Meanwhile, the number and protein sequences of ZmGSKs were more similar to those of OsGSKs than to those of AtASKs, indicating that GSKs from maize have a closer relationship with GSKs from rice than with those from Arabidopsis. This shows that the GSK gene is not highly conserved in the evolution of species and diversified during its evolution.

Comparative and Phylogenetic Analyses of ZmGSKs
In order to precisely reveal the evolutionary relationships of the ZmGSK proteins, we performed phylogenetic analyses on 30 GSK proteins, including 11 maize, 10 Arabidopsis, and 9 rice GSKs, using the neighbor-joining method ( Figure 1). All proteins were categorized into four groups (I, II, III, and IV). Meanwhile, the number and protein sequences of ZmGSKs were more similar to those of OsGSKs than to those of AtASKs, indicating that GSKs from maize have a closer relationship with GSKs from rice than with those from Arabidopsis. This shows that the GSK gene is not highly conserved in the evolution of species and diversified during its evolution.

The Analyses of Chromosomal Locations, Gene Structures, and Conserved Motifs
The ZmGSK genes were unevenly distributed among the 10 chromosomes of the maize genome ( Figure 2). The number of genes on each chromosome ranged from 0 (chromosomes 2 and 7) to 2 (chromosomes 3, 5, and 8), and the other five ZmGSK genes were evenly distributed in chromosomes 1, 4, 6, 9, and 10 (one gene on each chromosome).
The ZmGSK genes were unevenly distributed among the 10 chromosomes of the maize genome ( Figure 2). The number of genes on each chromosome ranged from 0 (chromosomes 2 and 7) to 2 (chromosomes 3, 5, and 8), and the other five ZmGSK genes were evenly distributed in chromosomes 1, 4, 6, 9, and 10 (one gene on each chromosome). The sequence alignment of the ZmGSK proteins is shown in Figure S1. The ~280amino acid is near the carboxyl terminus. The phylogenetic tree branches of maize GSKs were categorized into four groups (Figure 3a). MEME was used to analyze conserved motifs in the 11 ZmGSK proteins (Figure 3b, Table S1). In total, we identified 15 putative conserved motifs, with 6-50 residues and E-values < 1.2×10 −22 . Motifs 1-8 were found in almost all members of the ZmGSK family. All GSKs contained a highly conserved S_TKc domain ( Figure S2). In addition, ZmGSK3 contained two KH domains and a few PPR domains. To gain insights into the structural features of the predicted ZmGSK genes, we determined their intron/exon distribution ( Figure 3c). The structural analysis indicated that the coding regions of all maize GSK genes were interrupted by 10-12 introns. ZmGSK3 was an exception to this, with 19 introns in the gene. The UTR regions of ZmGSK3 were undetermined, while other maize GSK genes had UTR-region information in the genomic annotation files. Overall, highly similar gene structures were observed for the four groups of GSK family genes, especially for groups I, II, and III, according to the number and length of exons. The sequence alignment of the ZmGSK proteins is shown in Figure S1. The~280-amino acid is near the carboxyl terminus. The phylogenetic tree branches of maize GSKs were categorized into four groups ( Figure 3a). MEME was used to analyze conserved motifs in the 11 ZmGSK proteins ( Figure 3b, Table S1). In total, we identified 15 putative conserved motifs, with 6-50 residues and E-values < 1.2 × 10 −22 . Motifs 1-8 were found in almost all members of the ZmGSK family. All GSKs contained a highly conserved S_TKc domain ( Figure S2). In addition, ZmGSK3 contained two KH domains and a few PPR domains. To gain insights into the structural features of the predicted ZmGSK genes, we determined their intron/exon distribution ( Figure 3c). The structural analysis indicated that the coding regions of all maize GSK genes were interrupted by 10-12 introns. ZmGSK3 was an exception to this, with 19 introns in the gene. The UTR regions of ZmGSK3 were undetermined, while other maize GSK genes had UTR-region information in the genomic annotation files. Overall, highly similar gene structures were observed for the four groups of GSK family genes, especially for groups I, II, and III, according to the number and length of exons.

Cis-Acting Elements of ZmGSK Genes in Maize
To further understand how ZmGSK genes are regulated, we detected their upstream promoter sequences (2000 bp upstream of the coding region) and subjected them to cisacting element prediction using PlantCARE. Twenty cis-acting elements, including ab-

Expression Patterns of ZmGSK Genes in Different Tissues and Organs
To further define their distinct biological functions in development, the expression patterns of all ZmGSKs were analyzed using public RNA-seq data [30] (Figure 5a) and semiquantitative RT-PCR (Figure 5b). Our results showed that ZmGSK9 was expressed in different tissues and that the expression was much higher in mature leaves, especially after pollination, indicating that its function is mainly related to source regulation. ZmGSK5 was highly expressed in anthers, showing that this gene is involved in anther development. ZmGSK4 and ZmGSK11 exhibited ubiquitous expression in all tissues, but ZmGSK4 was expressed at higher levels in the endosperm, as shown by RT-PCR. The transcript level of ZmGSK8 was much higher in the cob and silk tissues than in other tissues. ZmGSK1 and ZmGSK3 expression were similar in roots, stems, and leaves. ZmGSK1 was expressed at a higher level in silk and decreased with seed development, while ZmGSK3 expression was increased with seed and endosperm development, suggesting that ZmGSK1 is involved in fertilization and that ZmGSK3 focuses on endosperm filling and maturation. However, ZmGSK3 was detected at a low level in all tissues using RT-PCR. The transcript levels of ZmGSK7 were much higher in silk and developing embryos.

Expression Patterns of ZmGSK Genes in Different Tissues and Organs
To further define their distinct biological functions in development, the expression patterns of all ZmGSKs were analyzed using public RNA-seq data [30] (Figure 5a) and semiquantitative RT-PCR (Figure 5b). Our results showed that ZmGSK9 was expressed in different tissues and that the expression was much higher in mature leaves, especially after pollination, indicating that its function is mainly related to source regulation. ZmGSK5 was highly expressed in anthers, showing that this gene is involved in anther development. ZmGSK4 and ZmGSK11 exhibited ubiquitous expression in all tissues, but ZmGSK4 was expressed at higher levels in the endosperm, as shown by RT-PCR. The transcript level of ZmGSK8 was much higher in the cob and silk tissues than in other tissues. ZmGSK1 and ZmGSK3 expression were similar in roots, stems, and leaves. ZmGSK1 was expressed at a higher level in silk and decreased with seed development, while ZmGSK3 expression was increased with seed and endosperm development, suggesting that ZmGSK1 is involved in fertilization and that ZmGSK3 focuses on endosperm filling and maturation. However, ZmGSK3 was detected at a low level in all tissues using RT-PCR. The transcript levels of ZmGSK7 were much higher in silk and developing embryos. ZmGSk6 and ZmGSK10 exhibited similar expression patterns in the root, stem, leaf, cob, seed, and endosperm, but ZmGSK10 had higher expression levels in the embryo. ZmGSK2 was detected at a high level in the endosperm, while it was missing from the MaizeGDB database.

Protein-Protein Interaction of ZmGSKs and BZR1
Previously, BZR1 was identified to interact with BIN2 and could be phosphorylated by BIN in Arabidopsis [31,32]. As BZR1 is a transcription factor and has to activate transcription to function, we investigated whether the fused BD-GSK3 protein interacted with AD-BZR1 using a yeast two-hybrid experiment. The results showed that ZmGSK1, ZmGSK3, ZmGSK4, ZmGSK7, ZmGSK8, and ZmGSK11 could interact with BZR1 ( Figure  6a). The other ZmGSK proteins did not interact with BZR1. The ZmGSKs-ZmBZR1 interaction was further verified with bimolecular fluorescence complementation (BiFC) assays. ZmGSKs were fused to a C-terminal yellow fluorescent protein (YFP) fragment (ZmG-SKs-cYFP), and the ZmBZR1 protein was fused to an N-terminal YFP fragment (ZmBZR1-cYFP). When fused ZmGSKs-cYFP was co-expressed with ZmBZR1-nYFP in onion, a YFP fluorescence signal was detected in transformed-cell nuclei (Figure 6b). No fluorescence was observed in negative-control experiments (Figure 6b). These results suggest that ZmGSKs kinase physically interacts with the ZmBZR1 transcription factor in plant-cell nuclei to participate in the BR pathway. To investigate the response of these genes to BR and ABA starvation, we harvested Mo17 endosperms 9 or 10 DAP; grew them on normal, BR, and ABA 1/2 MS media for 48 h; and subjected them to qRT-PCR analysis. Our data showed that eight genes (ZmGSK3, ZMGSK4, ZmGSK6, ZmGSK7, ZmGSK8, ZmGSK10, and ZmGSK11) were downregulated after BR treatment for 48 h compared with after treatment with CK ( Figure 5c). Meanwhile, ZmGSK1 and ZmGSK9 were upregulated by more than two-fold after BR treatment. Five genes (ZmGSK4, ZmGSK5, ZmGSK6, ZmGSK9, and ZmGSK10) were upregulated, and four genes (ZmGSK2, ZmGSK3, ZmGSK7, and ZmGSK11) were downregulated after ABA treatment, compared with after CK treatment.

Protein-Protein Interaction of ZmGSKs and BZR1
Previously, BZR1 was identified to interact with BIN2 and could be phosphorylated by BIN in Arabidopsis [31,32]. As BZR1 is a transcription factor and has to activate transcription to function, we investigated whether the fused BD-GSK3 protein interacted with AD-BZR1 using a yeast two-hybrid experiment. The results showed that ZmGSK1, ZmGSK3, ZmGSK4, ZmGSK7, ZmGSK8, and ZmGSK11 could interact with BZR1 (Figure 6a). The other ZmGSK proteins did not interact with BZR1. The ZmGSKs-ZmBZR1 interaction was further verified with bimolecular fluorescence complementation (BiFC) assays. ZmGSKs were fused to a C-terminal yellow fluorescent protein (YFP) fragment (ZmGSKs-cYFP), and the ZmBZR1 protein was fused to an N-terminal YFP fragment (ZmBZR1-cYFP). When fused ZmGSKs-cYFP was co-expressed with ZmBZR1-nYFP in onion, a YFP fluorescence signal was detected in transformed-cell nuclei (Figure 6b). No fluorescence was observed in negative-control experiments (Figure 6b). These results suggest that ZmGSKs kinase physically interacts with the ZmBZR1 transcription factor in plant-cell nuclei to participate in the BR pathway.

Discussion
The plant GSK3 protein family plays an important role in growth and development [33], and participates in biological-and abiotic-stress responses [34]. BIN2, a member of the GSK3 family, has been reported as a negative regulator of BR signaling in different species [10,[35][36][37]. Our findings could accelerate the further functional characterization of ZmGSK genes and provide a framework to clarify the roles of these genes in biology and hormone response.

Classification and Phylogenetic Analysis of the Maize GSK Family
A number of GSK genes have been identified in Arabidopsis [4,38,39] and rice [13]. In contrast with mammals, GSK proteins are much more diverse in plants; for example, there are ten GSK isoforms in Arabidopsis [9] and nine in rice [40]. Here, we found eleven GSK genes in maize, which might reflect the larger size of the maize genome.
Our phylogenetic comparison of the GSK proteins included maize, rice, and Ara-

Discussion
The plant GSK3 protein family plays an important role in growth and development [33], and participates in biological-and abiotic-stress responses [34]. BIN2, a member of the GSK3 family, has been reported as a negative regulator of BR signaling in different species [10,[35][36][37]. Our findings could accelerate the further functional characterization of ZmGSK genes and provide a framework to clarify the roles of these genes in biology and hormone response.

Classification and Phylogenetic Analysis of the Maize GSK Family
A number of GSK genes have been identified in Arabidopsis [4,38,39] and rice [13]. In contrast with mammals, GSK proteins are much more diverse in plants; for example, there are ten GSK isoforms in Arabidopsis [9] and nine in rice [40]. Here, we found eleven GSK genes in maize, which might reflect the larger size of the maize genome.
Our phylogenetic comparison of the GSK proteins included maize, rice, and Arabidopsis, and all 30 GSK genes were divided into four subgroups (Figure 1). The diverse subgroups of gene motifs were different, but there was a highly conserved motif common to all subgroups (Figure 3b). For example, subgroup I of ZmGSKs contained a conserved domain (motif 12) in the N terminus and motif 11 only existed in subgroup III. Meanwhile, the relationship of maize GSK genes in the two trees was comparable. For example, the branches of ZmGSK1, ZmGSK4, ZmGSK7, ZmGSK8, and ZmGSK11 in the maize NJ tree were also included in subgroup II in the total NJ tree. A genetic analysis of T-DNA insertional mutants of group II in Arabidopsis revealed that their function is redundant in BR signaling but that BIN2 plays a dominant role [15,41]. In addition, both OsGSK1 and OsGSK2 are orthologs of ASK21 (BIN2), whereas only OsGSK2 appeared to play a predominant role in BR signaling [35,42]. Thus, OsGSK2 is considered to be the rice counterpart of Arabidopsis BIN2 because it has the highest sequence similarity [42]. A previous study found that ZmASK1, named ZmGSK3 in our study, belongs to group IV in maize. However, BIN2 belongs to group II in rice and Arabidopsis. Thus, the maize counterpart of BIN2 remains to be identified. Here, we found that both ZmGSK1 and ZmGSK8 have a high sequence similarity with BIN2 of rice and Arabidopsis. Therefore, one of these two genes is likely maize BIN2. Our logical classification may be helpful in further studies of GSK protein function.

Diverse Expression Patterns and Hormone Responses Generate Differential Functions
Detailed analyses of gene expression patterns and hormone responses can provide considerable insight into the critical regulatory points in which they are implicated [43][44][45]. Based on the maizeGDB database and our semi-quantitative RT-PCR results, the expression patterns of the candidate GSK genes were differentiated. Most of the GSK genes in maize were expressed in various tissues during plant development, but some genes showed higher transcript levels in specific organs (Figure 5a). For example, ZmGSK2, ZmGSK4, and ZmGSK7 were mainly expressed in the endosperm, whereas ZmGSK1, ZmGSK8, and ZmGSK10 were consistently weakly expressed in all analyzed organs. The observed expression patterns were consistent with a previous analysis of ZmASK1 [46]. Diverse expression patterns may imply differential functions. For example, AtSK11 and AtSK12 are mainly expressed during flower development and are required for the patterning of gynecium [38]. In an OsGSK1 knockout mutant, OsGSK1 expression was the highest in young panicles, and spikelets were larger [13]. In our study, the relatively high expression of maize GSK genes in reproductive organs was more obvious, indicating that these GSK genes are involved in the development of reproductive organs. Meanwhile, we found that some subgroups of the GSK gene showed similar expression patterns, such as ZmGSK1 and ZmGSK8, indicating that they perform similar functions. However, the expression of ZmGSK1 was higher than that of ZmGSK8. These data suggest that ZmGSK1 is more important than ZmGSK8 and may be the counterpart of BIN2 in maize. In addition, group I members showed diverse expression patterns, similar to Arabidopsis [4]. ZmGSK2 was mainly expressed in the endosperm, and ZmGSK5 was expressed in anthers, while ZmGSK9 was detected at high levels in leaves.
In our study, we focused on identifying the genes involved in seed development. Thus, we analyzed the expression of eleven genes under BR and ABA conditions. The quantitative results indicated that all the GSK genes were induced in seeds by BR (except for ZmGSK5) and that nine GSK genes were induced under ABA conditions (Figure 5b). Meanwhile, for most of the GSK genes, the expression levels were the opposite under BR and ABA conditions, and a high BR expression was accompanied by low expression under ABA, and vice versa. Although these results indicated that GSK genes have an antagonistic effect on BR and ABA, which regulate seed development, the effect of the BR and ABA signaling pathways is apparently complicated [47]. These results are consistent with previous research on seed germination [48,49]. In addition to studying expression, other studies are needed to clarify the biological functions of the isolated genes with respect to their diverse expression patterns and hormonal responses.
Previous studies on the function of GSK proteins mainly focused on interactions with BZR1 [20]. In our study, we found that members of groups II and IV in maize interacted with BZR1 (Figure 6a). A study on Arabidopsis also found that members of group II ASKs could bind to BZR1 as a negative regulator in BR signaling [50]. Meanwhile, members of group III ASKs in Arabidopsis interact with BEH2 [51]. Thus, it is suggested that members of groups II and IV mediate the hormonal pathway through interaction with BZR1 and that groups I and III may interact with other transcription factors to participate in hormonal responses. Previous studies showed that ZmBZR1could bind to the KRP6 and GRACE promoter to regulate seed size [52]. Here, we found that ZmGSK1, ZmGSK4, ZmGSK7, and ZmGSK8, which could interact with ZmBZR1, were expressed in the endosperm. These results show that ZmGSK1, ZmGSK4, ZmGSK7, and ZmGSK8 play important roles in seed development.

Plant Materials
Seeds from maize inbred line Mo17 were obtained from Maize Research Institute of Sichuan Agricultural University (latitude, 30.705207; altitude, 511 m). To analyze tissue-specific gene expression patterns in Mo17, plants were grown under standard cropmanagement conditions in Chengdu, China. Roots, stems, and leaves were collected at the shooting stage. Additionally, seeds, endosperms, and embryos were harvested 15 days after pollination (DAP). Anthers and silks were gathered during plant pollination. All collected samples were immediately frozen in liquid nitrogen and stored at −80 • C until use.

Identification of Maize GSK Genes
Putative ZmGSK genes were identified via homologous alignment. The ZmGSK family genes from Oryza sativa and Arabidopsis thaliana were downloaded from the Tair database (https://www.arabidopsis.org (accessed on 1 July 2019) and the rice annotation project database (http://rapdb.dna.affrc.go.jp/download/irgsp1.html (accessed on 1 July 2019), respectively. The versions of the maize genome sequence obtained from NCBI (https://www.ncbi.nlm.nih.gov/ (accessed on 1 July 2019) were used to identify ZmGSK proteins and their corresponding nucleotide sequences. All the A. thaliana ASK genes and encoded protein sequences, which were used as queries, were blasted against the NCBI database to search for GSK homologs in maize, with the E-value cutoff set as 1.00 × 10 −10 and the coverage ratio as 50%. The Hidden Markov Model (HMM) profile of Pkinase (PF0069.24) was obtained from the Pfam protein families database (http://pfam.xfam.org/ (accessed on 28 July 2019) and used to verify maize GSKs. The candidate genes were further checked using pfamscan and the Pfam A database to ensure the presence of the Pkinase domain (PF0069.24).

Analysis of Chromosomal Localization
The maize GSK genes were located on the corresponding chromosomes using MapInspect software (https://mapinspect1.software.informer.com/ (accessed on 28 July 2019) according to the start positions indicated in the maize database.

Multiple-Sequence Alignment and Phylogenetic Analysis
The full-length amino acid sequences of GSK proteins from maize, rice, and Arabidopsis were aligned using MUSCLE (http://www.ebi.ac.uk/Tools/msa/muscle/ (accessed on 31 July 2019) and saved in the ClustalW format. Then, the unrooted phylogenetic tree was constructed using the neighbor-joining method in MEGA5.10 software with the bootstrap test replicated 1000 times.

Analysis of Promoter Cis-Regulatory Elements
To investigate the putative cis-acting elements of candidate ZmGSK genes, their promoter sequences (2000 bp upstream of the initiation codon "ATG") were obtained with BLAST searches of the maize genome data using whole gene IDs. Potential cis-acting regulatory elements of the extracted sequences were subsequently examined using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/ (accessed on 13 August 2019) [53].

Analysis of Biophysical Properties of the GSK Genes
The maize GSK protein sequences were analyzed using the ExPASy-ProtParam tool (http://web.expasy.org/protparam/ (accessed on 17 August 2019) to calculate the number of amino acids, molecular weight, and theoretical pI. Meanwhile, the subcellular localization of these genes was predicted using ProComp 9.0 (http://www.softberry. com/berry.phtml?group=programs&subgroup=proloc&topic=protcomppl (accessed on 17 August 2019).

Analysis of Expression Patterns of the GSK Genes
Expression data of the eleven maize GSK genes in different maize tissues were obtained from the MaizeGDB database. The expression patterns are presented as a heat map, which reflects log 2 (RPKM + 1), with red, white, and blue indicating high, medium, and low expression, respectively. Total RNA was extracted using TRIzol reagent (Invitrogen, Waltham, MA, USA). The isolated RNA (1.5 µg) was treated with DNase I, and cDNA synthesis was conducted using reverse-transcription PCR (TaKaRa, Dalian, China). The transcription levels of ZmGSK genes in various tissues were measured with semiquantitative RT-PCR, and ZmGSK gene expression levels following BR (10 nM) and ABA (100 µM) treatment of maize seeds (9-10 DAP) were measured using qRT-PCR. The primers used for semiquantitative RT-PCR and qRT-PCR are shown in Table S2.

Yeast Two-Hybrid Assay
Yeast two-hybrid assays were performed as previously described. Full-length coding sequences of ZmGSKs were cloned into the pGBKT7 vector (BD) and transformed into yeast strain AH109 to test for auto-activation. Yeast on SD/-Trp and SD/-Trp-Ade-His agar plates was grown at 28 • C for 3 days. For the protein-protein interaction assay, ZmBZR1 was ligated to the pGADT7 vector (AD). Plasmid pGADT7-ZmBZR1 with pGBKT7-ZmGSKs were co-transformed into AH109 cells. Yeast strain AH109 was planted on SD/-Trp-Leu and SD/-Trp-Leu-Ade-His for 3 days. The primers used for cloned genes are shown in Tables S3 and S4.

Biomolecular Fluorescence Complementation (BiFC) Assays
The cDNA sequences of the ZmGSK1,4,7,8,11 fragments were PCR-amplified and cloned into pSAT6-cEYFP, respectively. Full-length ZmBZR1-encoding sequences were inserted into pSAT6-nEYFP to generate an N-terminal fusion with nEYFP. For the BiFC protein-protein interaction analysis, the above constructs were introduced into onion with a helium biolistic gun transformation system. The bombarded samples were then cultured for from 24 to 48 h at 28 • C. Then, the samples were observed using BX61 fluorescent microscopy. The primers used for cloned genes are shown in Table S5.

Conclusions
In summary, we identified and cloned 11 GSK genes in maize, which were distributed into four groups. We found that GSK genes are highly conserved among Arabidopsis, rice, and maize. Additionally, groups II and IV play important roles in the regulation of seed development by interacting with BZR1. Multiple lines of evidence suggest that ZmGSK1 plays a dominant role in BR signaling as the counterpart of BIN2. Overall, our results indicate that ZmGSK proteins are generally involved in the interaction with the ZmBZR1-mediated hormone pathway. Further, we show that some GSK gene members are involved in the BR, ABA, and BR-ABA hormone pathways.


s2fieldsofstudy:
Agricultural And Food Sciences
Biology

provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
In vitro antimicrobial appraisal of the potentials of Morinda lucida against some selected bacteria

The bioactive compounds from plant extract represent a potent alternative to combating the increasing antimicrobial resistance (AMR) to synthetic drugs. This study aimed to investigate the phytochemical, and antibacterial activities of Morinda lucida . Leaves of M. lucida was extracted successively using n - hexane, acetone, methanol and distilled water. The crude plant extracts were evaluated for their antibacterial activity against pathogens found in food and other sources such as Staphylococcus aureus , Bacillus subtilis , Pseudomonas aeruginosa , Salmonella enterica serovar Typhi, Escherichia coli , and Klebsiella spp. using the agar well diffusion technique. Minimum inhibitory concentration (MIC) was carried out using the microbroth dilution technique. The phytochemical analysis showed the presence of alkaloids and sterols in all extracts. Flavonoid was only present in the n - Hexane extract. The result showed that n - Hexane extract had the widest zone of inhibition against P. aeruginosa - 13.67 mm (12.5 mg/mL MIC), B. subtilis - 15 mm (12.5 mg/mL MIC), while acetone extract was most active against S. enterica ser. Typhi - 24.67 mm (6.25 mg/mL MIC), S. aureus - 24.83 mm (3.13 mg/mL MIC), E. coli -13.33 mm (1.56 mg/mL MIC) and Klebsiella spp. - 18.33 mm (12.5 mg/mL MIC). These findings indicated that n - Hexane and acetone extract of M. lucida exhibited significant antimicrobial activity. Therefore, further research should be carried out to isolate the bioactive compounds with potentials in developing new drug leads against drug - resistant microbial strains. This current study revealed good antibacterial activity of the different plant extracts against selected gram - positive and gram - negative bacteria. The presence of phytochemicals such as alkaloids, saponin, cardiac glycoside, sterols, phenols, terpenoids, quinones, flavonoids and tannins in these plant extracts could be responsible for the antimicrobial activities observed. antimalarial, anticancer antiviral,

Introduction
Antimicrobial resistance (AMR) is a situation where microorganisms evolve to a state of non-susceptibility upon exposure to established antimicrobial agents resulting in persistent infection by resistant microbes (Hay et al., 2018). The impact of AMR on economies and health care systems of a nation is significant and limits the productivity of patients and handlers as a result of the extended hospital admission which may require a sophisticated treatment approach that may be more expensive (Dadgostar, 2019). The continued occurrence of AMR and multidrug resistance (MDR) in pathogenic microbes are of great concern as there are no effective antimicrobials to combat such pathogens (Tanwar et al., 2014;Aslam et al., 2018;Hay et al., 2018;Talebi Bezmin Abadi et al., 2019). Some of these pathogens of concern are transmitted through contaminated food. Thus, the need to find new antimicrobial agents that are effective against pathogenic microbes. Before the 16 th century, plants were used for prophylactic and therapeutic purposes and around the 16 th century the advent of iatrochemistry changed the approach to the treatment of patients. Recent trends have seen a reduction in the efficacy of synthetic drugs and their increasing detrimental side effect. This has resulted in sourcing for treatment with the use of natural products (Petrovska, 2012). Medicinal plants have been engaged in treating ailments in Africa for centuries. In line with World Health Organization's (WHO) report, nearly 80% of people living in developing countries adopt the administration of herbal medicine in primary healthcare practice (World Health Organization, 2013).
South-west Nigeria is host to Morinda lucida. The plant is rich in nutrients and is available all year round. It has been reported to be rich in Vitamin A, E and K. It has been reported to have antibiotic, anti-parasitic and antiviral components. Infusions or plasters and decoctions of parts of the plant are used to treat different fevers, stomach aches, diabetes, cerebral congestion, gonorrhoea, dysentery, ulcers, hypertension, and leprosy (Agyare et al., 2013;Adeleye et al., 2018). The plant is locally known in the southwest region of Nigeria as "Oruwo" (Osuntokun, 2015). This current study explored continuous successive extraction for stratification of extract based on polarity, to establish the antimicrobial properties of the different phytochemical components of M. lucida.

Materials and methods
The leaves of M. lucida were collected from the Omu-Aran forest in Kwara State, Nigeria. Identification of the plants was done at the Department of Plant Biology, the University of Ilorin and a voucher specimen was archived with the reference number UILH/001/913/2021 (Morinda lucida).

Extraction
A total of 600 g of M. lucida leaves were washed and dried at 16-18°C. The leaves were pulverized to a fine powder and stored in an air-tight container. The pulverized material was extracted using the modified method by Rubab et al. (2021). The dried pulverized plant leaves were extracted successively in n-hexane, acetone, methanol and water, by soaking in 6 L of each solvent. Maceration was done in each solvent for three days while shaking at 100 rpm to increase the extraction rate at ambient temperature (27±2.5°C). Extraction solvent was removed from each extract with the aid of a rotary evaporator under reduced pressure and the concentrated extracts were weighed and stored in a glass container (Rubab et al., 2021).

Phytochemical screening 2.2.1 Qualitative phytochemical screening
Crude solvent-free extracts of M. lucida leaves were analyzed for the presence of phytochemicals -alkaloids, saponins, cardiac glycoside, sterols, phenols, terpenoids, quinones, flavonoids and tannins using the established method described by .

Determination of total alkaloid
The gravimetric technique previously described in the literature was modified for use in this work (Onochie et al., 2020). A total of 5 g of each extract was dissolved in 50 mL of 10% ethanolic acetic acid solution. The solution was stirred for about 10 mins and allowed to stand for a period of 4 hrs before filtering. The filtrate collected from the solution was allowed to evaporate to 25% of its initial volume after which the alkaloid present in the solution was precipitated with concentrated NH 4 OH (drop-wise addition). The precipitate was removed from the solution with a filter paper and washed with 1% NH 4 OH solution. The resulting precipitate was dried at 60°C in a hot air oven for 30 mins to evaporate all liquid from the precipitate and reweighed to get the alkaloid content of each extract.

Determination of total phenolic content
The total phenolic content of the leaf extracts was quantified using the Folin Ciocalteu technique as reported by Singleton et al. (1999). The standard was gallic acid. A gallic acid calibration curve was made using varying concentrations of 0.2, 0.4, 0.6, 0.8 and 1.0 mg/mL. Then, 50 μL of crude extract was measured into 250 μL Folin-Ciocalteu reagent and mixed with 3 mL of double distilled water and allowed to stand for 5 mins on the bench. About 750 μL of 20% Na 2 CO 3 was added to the standing mixture and mixed for 2-4 mins after which it was incubated at room temperature for 30 mins. A UVvisible spectrophotometer was used to measure the absorbance at 760 nm in triplicate. The total phenolic content was evaluated from the calibration curve.

Determination of total tannin content
The total tannin content of the leaf extracts was derived using a slightly modified Folin-Ciocalteu assay technique. A Gallic acid calibration curve was made using varying concentrations of 0.2, 0.4, 0.6, 0.8 and 1.0 mg/mL. Then, 5 mL of already prepared leaf extract, was added to 5 mL of distilled water, Na 2 CO 3 (35%) and Folin-Ciocalteu reagent (500 μL) was mixed for 15 mins and allowed to stand for 30 mins at ambient temperature. A UV-visible spectrophotometer was used to measure the absorbance for standard solutions against blank at 725 nm and the total tannin content was extrapolated from the calibration curve and expressed as mg/g (Singleton et al., 1999).

Determination of total saponin content
A total of 5 g of the leaf extract was dissolved in 100 mL of 20% acetic acid (in ethanol) and left to stand in a water bath at 50°C for 24 hrs. The resulting solution was filtered and concentrated in a water bath to a quarter of the initial volume. Concentrated NH 4 OH was added drop -wise to the solution until it was completely precipitated. The solution containing the precipitate was allowed to stand for the precipitate to settle before the precipitate was filtered and weighed. The saponin content was eISSN: 2550-2166 © 2022 The Authors. Published by Rynnye Lyan Resources FULL PAPER weighed and calculated as mg/g of the sample .

Determination of total terpenoid content
A total of 2 g of the leaf extract were weighed into a beaker containing 50 mL of 95% ethanol for 24 hrs. The leaf extract was filtered and the collected filtrate was reextracted using petroleum ether (60-80 o C) and concentrated to dryness. The resulting mass from drying was taken as total terpenoids and estimated in mg/g (Bashua and Oladunmoye, 2017).

Antimicrobial screening 2.4.1 Antimicrobial activity of extracts
The antimicrobial potential of all crude extracts (200 mg/mL) was tested using the agar well diffusion technique. A modified method by Zofuo et al. (2019) was used. A 24-hr old culture on nutrient broth was used as inoculum and the closed mesh was streaked on solidified Mueller-Hinton agar plates for every organism in replicates. A 9 mm cork borer was used to make four wells in the Mueller-Hinton agar. Each well was filled with 100 µL of each leaf extract. Ciprofloxacin was used as the positive control and DMSO as the negative control. The prepared plates were allowed to stand for the diffusion of the extract into the medium before incubating at 37°C for 18 -24 hr. The plates were observed for zones of inhibition and the antimicrobial activity was recorded as the mean of the diameters of the zones of inhibition (Zofou et al., 2019).

INT colourimetric assay for MIC determinations
The MIC of leaf extracts of M. lucida on S. aureus (ATCC 6538), B. subtilis (ATCC 6633), P. aeruginosa (ATCC 9027), S. enterica serovar Typhi (ATCC 20971), E. coli (ATCC 25922) Klebsiella spp. were done using INT colourimetric micro broth dilution as described by Ohikhena et al. (2017) with some modification. The leaf extracts were dissolved in DMSO to a concentration of 200 mg/mL. Mueller Hinton Broth (100 μL) was dispensed into all 96 wells of the microtitre plates. A two -fold dilution of plant extract (100 μL) was made across the well rows, till the tenth dilution was achieved. A total of 20 μL broth containing the test organism was introduced into each well. The plates were covered and incubated at 37°C for 16-18 hrs. The assay was done in replicates. DMSO served as the negative control. 10 μL of 0.2 mg/mL p-iodonitrotetrazolium chloride (INT) was added to all wells after 18 h and further incubated for another 30 mins. The MIC was recorded as the least concentration in wells that could not reduce the yellow INT solution to pink (Ohikhena, Wintola, and Afolayan, 2017).

Qualitative phytochemical screening
The extraction procedure yielded four different extracts one for each extraction liquid used (n-hexane, acetone, methanol and water). The agents responsible for antimicrobial activities in plants are the phytochemicals. The qualitative phytochemical analysis of the leaf extract from the four solvents was carried out. The analysis revealed the presence of Alkaloids and sterols in all extracts and water extract has the highest number of phytochemicals including alkaloids, saponins, cardiac glycoside, sterol, phenols terpenoids, quinones and tannins lacking only flavonoids. N-hexane extract tested positive for only a few phytochemicals, namely, alkaloids, sterols and flavonoids which could be a result of the polarity of the solvent (Table 1).

Quantitative phytochemical screening
The availability and quantity of the phytochemical constituent of plants play a major role in their activity when used in local concoctions. Some selected phytochemicals present in the extracts were quantified (Table 2) and it was observed that the acetone extract had the highest quantity of phenols (0.06±0.34 mg/g) and tannins (0.06±0.32 mg/g) per gram of extract. The methanol extract had the highest quantity of Alkaloids (0.06±0.004 mg/g) and terpenoids (0.05±0.01 mg/g), while the water extract had the highest quantity of saponins (0.12±0.1 mg/g) ( Table 2).  (Table  3).

INT colourimetric assay for MIC determinations
The MIC of the leaf extract was examined against six different isolates and the MIC varied against the different isolates. The n-hexane leaf extract had the lowest MIC (6.25 mg/mL) against S. enterica ser. Typhi (Table 4) while the acetone leaf extract had the lowest MIC (12.5 mg/mL) against P. aeruginosa (Table 6). In the case of S. aureus (Table 5), Acetone extract had the lowest MIC (3.125 mg/mL) next to which was n-hexane extract (6.25 mg/mL). The MIC of all four extracts in broth against Klebsiella spp. was 12.5 mg/mL showing similar potentials across all extracts (Table 7). Acetone extract had a MIC of 1.56 mg/mL against E. coli while other extracts had MICs of 12.5 mg/mL ( Table 8). Three of the leaf extracts (in n-hexane, acetone and water) shared the lowest MIC (12.5 mg/mL) against B. subtilis (Table 9).

Discussion
Time and evolution have allowed bacterial pathogens to attain a state of resistance to previously potent antibacterial agents. Some of these organisms access the human host through food, water and other sources. These trends call for the discovery of new, effective and safe antibacterial agents. Plants extracts have been reported to be a suitable source for new antibacterial agents due to their use in traditional medicine for centuries. Plants are rich in phytochemicals that could have a prophylactic or therapeutic effect on health. They have been studied widely to ascertain potency and to decipher how they confer these health benefits.  This current study revealed good antibacterial activity of the different plant extracts against selected gram-positive and gram-negative bacteria. The presence of phytochemicals such as alkaloids, saponin, cardiac glycoside, sterols, phenols, terpenoids, quinones, flavonoids and tannins in these plant extracts could be responsible for the antimicrobial activities observed. Reports have shown that alkaloids exhibit antibacterial, antiproliferation, insecticidal, antiviral and antimetastatic effects on diverse cancer cells both in vitro and in vivo (Qiu et al., 2014;Hussain et al., 2018;Chikowe et al., 2020). The glycoside and saponin constituent of the plant leaf extracts indicates the plant's potential to confer protection against pathogens and exhibit antimicrobial, antiulcer, anti-inflammatory, and antidiarrheal activities (Prakash et al., 2020). Terpenoids have been implicated to possess antimicrobial, anticarcinogenic, hepaticidal, antimalarial, anticancer activities and anti-ulcer activities. Flavonoids (e.g. rutin and quercetin) have antiviral, anti-inflammatory and antihistaminic effects. The presence of this phytochemical may be responsible for the health benefits derived from M. lucida when used in traditional medicine to treat diverse ailments as claimed by the African Traditional Medicine practitioners.
Quantitative analysis of the phytochemicals present in the leaf extract of this plant revealed the acetone extract had the highest phenolic content (total phenols and tannins) per mg of GAE/g of leaf extract, which could be responsible for the superior antimicrobial effect when compared to extract from other solvents. Plantderived phenolics including tannins have been implicated for their diverse biological activities which include their inhibitory and bactericidal activities against selected bacteria of medical importance (Mandal et al., 2017;Bouarab-Chibane et al., 2019;Maisetta et al., 2019).
The extraction solvent plays an important role in identifying the active fractions or compounds in plant material. The extraction of the plant leaf was done in solvent successions, starting from the least polar solvent to the most polar. Though the plant leaf is traditionally extracted using water, this study reported a better activity from acetone extracts in-vitro. This is in tandem with reports that most antimicrobial compounds found in plant materials are non-polar or have intermediate polarity, thus making acetone extract most potent against   bacteria. This may be a result of the polar and nonpolar characteristics of acetone as an extraction solvent (Eloff, 2019).
The agar well antimicrobial susceptibility testing revealed that the acetone extract had a higher zone of inhibition against S. enterica ser. Typhi, S. aureus, E. coli and Klebsiella spp. (24.67 mm, 24.83 mm, 13.33 mm and 18.33 mm respectively). The activity against S. aureus was almost similar to ciprofloxacin (26.33). The result is in line with the report from several works on the antimicrobial potentials of M. lucida (Fakoya et al., 2014;Temitope et al., 2016).
The acetone extracts had the lowest MIC of all extracts against Escherichia coli and Staphylococcus aureus (1.56 mg/mL and 3.1 mg/mL respectively) and were only equalled by n-hexane extract against S. enterica ser. Typhi and P. aeruginosa, the activities against E. coli are in tandem with earlier reports (Osuntokun, 2015;Olawuwo et al., 2019). This is a confirmation of Eloff's report on the antimicrobial activities of plant compounds characterized by intermediate polarity or non-polar (Eloff, 2019).

Conclusion
In conclusion, the findings in this work have been able to support the ethnobotanical claims on the use of M. lucida by rural dwellers in Nigeria. The activity of the extracts against all selected isolates which cuts across Gram-positive and negative bacteria of medical importance with its best activity against E. coli (acetone extracts) showed the plant can be considered as a source of a new antibiotic in the search for alternative antimicrobial agents in the continued rise of drugresistant bacteria. The method of extraction engaged in this study can be compared with others for the best approach to the extraction of active ingredients in the plant leaf. It is recommended that further studies be done to elucidate the active compound in the plant and isolate it for further work. This is the first time the extraction technique used in this work has been used for this plant and the result has shown that it can be considered a viable technique to study the activity of plant extracts against microbes. This work provides a phytochemical explanation for the herbal use of M. lucida in African Traditional Medicine. The effective inhibition of pathogenic bacteria by acetone extract at low MIC suggests that further detailed study could lead to the discovery of some potent antibiotic agents.

Conflicts of interest
The authors declare no conflict of interest.


text:
Dynamic Measurement of Water Use Level Based on SBM-DEA Model and Its Matching Characteristics with Economic and Social Development: A Case Study of the Yellow River Basin, China

: Enhancing the level of water use and alleviating the constraints of water shortage on economic and social development are powerful supports to realize the harmonious balance of water and economic society. In this study, the data envelopment analysis (DEA) window analysis method is applied to the study of water use level, and the SBM-DEA model (slack based measure, SBM) is combined to explore the spatial and temporal evolution characteristics of composite water use index (CWUI) in nine provinces from 2012 to 2018. The Malmquist index model is used to decompose the intrinsic causes of total factor productivity (TFP) changes, and the spatial matching degree calculation method is applied to study the matching degree between CWUI and economic and social development levels (E-SDL). The results showed that: (1) the overall trend of CWUI in the nine provinces from 2012 to 2018 was increasing, with signiﬁcant spatial variability in water use levels; (2) the improvement of TFP of water in the nine provinces was mainly driven by technological change (TC), and the main factor limiting the improvement of TFP of water was technical efﬁciency change (EC); (3) E-SDL of the nine provinces showed an increasing trend, with the spatial distribution characteristics of decreasing E-SDL of the downstream, midstream and upstream provinces in sequence; (4) the degree of matching between CWUI and E-SDL shows strong regional differences, with different types of matching.

Introduction
UNESCO (United Nations Educational, Scientific and Cultural Organization) defines the water resource as a source of water of a certain quality and quantity that is available and potentially available to meet the utilization needs of a site in the long term [1]. Water resources are the basic support for economic and social development, as well as a necessary guarantee for the construction of ecological civilization [2]. Since the 21st century, a green, coordinated, and sustainable water resource utilization model has gradually become a hot spot of concern for countries around the world [3]. The United Nations Sustainable Development Summit in 2015 released the 2030 Agenda for Sustainable Development, which identifies 17 Sustainable Development Goals (SDGs) covering economic, social, resource, and ecological dimensions [4]. It is internationally recognized that an important part of achieving SDGs is to conduct a comprehensive and in-depth exploration of water DEA model to assess the water use level of the three major industries in 30 provinces in mainland China. In studies on DEA models, traditional CCR and BCC models [33], super-efficient DEA models [34], SBM-DEA models [27,30,35], Malmquist models [24,36], etc., are involved. For example, Bai et al. [36], constructed a non-radial method based on the theory of Malmquist to dynamically measure changes in water use level of the Bohai Bay urban cluster in China. In recent years, with the gradual development of the green development concept. It is worth noting that DEA models considering non-expected output have started to appear more often in resource and energy use level studies. For example, Yang et al. [30] chose a non-expected output SBM model to measure the industrial water use level in mainland China. It can be found that the DEA method was widely applied in the field of water resources. The above high-quality research results provided a solid basis for the smooth development of subsequent studies.
Economic and social development and the level of resource and energy utilization are inextricably linked, especially under the dual effect of intensifying climate change and the impact of human activities. Whether the level of resource and energy utilization can match the level of economic and social development has become an important indicator to measure the sustainable development capability of the region. With the further implementation of sustainable development policies, the relationship between the level of resource and energy use and sustainable development has gradually become a hot issue of global concern. In this context, many scholars have adopted various methods to explore the relationship between the two. The subjects of this kind of study are mainly related to energy consumption [37][38][39], pollution emissions [40], land resource use [41,42], etc. For example, Sarkodie et al. [38] examined the effects of energy consumption on GHG emissions using U-test estimation methods and non-additive fixed effects panel quantile regression. On World Water Day, March 22, 2021, UNESCO released the World Water Development Report 2021 [43], which focuses on "the value of water" to explain the important role of water resources in economic and social development in different dimensions. There are not many studies on the relationship between water resources utilization and economic and social development, but some scholars conducted relevant studies. For example, Sun [44] constructed a system dynamics model to study sustainable water resources utilization considering economic development, and simulated the water resources supply and demand situation from 2005 to 2020 and the changes of future supply and demand gap. The above research results have contributed value in promoting sustainable economic and social development.
After combing through the representative literature in recent years, it can be found that the important results revealed by the above-mentioned studies have greatly promoted the development of the field of water use level and the field of water resources-economic and social relations, and expanded the scope and depth of research in this field. However, there are some shortcomings. For example, most studies on water use levels still use the DEA cross-sectional model to deal with long time series panel data. This defect leads to the comparison of different decision units only within a specific year, and the evolution characteristics of water use level in time series cannot be analyzed, which reduces the reference value of the results to some extent. Meanwhile, some studies used traditional radial DEA models to measure water use levels, which failed to fully consider slack variables [45]. Finally, studies on the matching relationship between water use level and economic and social development in the Yellow River Basin are still relatively lacking.
However, the combination of SBM-DEA model and Window-DEA model can effectively solve the above problems. This study attempts to apply the DEA window analysis method to the study of water use level, combining the SBM-DEA model to explore the spatial and temporal evolution characteristics of water use level in nine provinces in the Yellow River basin. It can realize the dynamic analysis of CWUI under the premise of considering slack variables. On this basis, the Malmquist index model is used to analyze the deep-seated reasons, leading to the variation of total factor productivity of water, and the spatial matching degree calculation method based on series distance is introduced for exploring the matching relationship between E-SDL and CWUI. This study aims to achieve the following objectives: (1) apply the SBM-DEA model combined with Window-DEA model to the study of water use level; (2) conduct dynamic evaluation of water use level in nine provinces with consideration of relaxation improvement; (3) clarify the deepseated reasons for the changes in total factor productivity of water in the nine provinces; (4) analyze the matching characteristics between the level of water utilization and the level of economic and social development in the nine provinces under different time scales.
The rest of this study is structured as follows: Part 2 introduces the four main research methods, Part 3 presents an overview of the study area and the selection of relevant indicators, and Part 4 analyzes the spatial and temporal evolution characteristics of CWUI, the changes in TFP of water, and the matching characteristics of CWUI and E-SDL. Part 5 is the conclusion and the outlook for future research.

SBM-DEA Model
The DEA model can be used to evaluate the relative efficiency of different DMUs (decision-making units) through means of a specific mathematical programming model. The basic principle is to determine the production frontier surface with the help of a linear model and to determine the relative efficiency value of each decision unit by comparing the deviation of each decision unit relative to the production frontier surface. DEA models are widely used in efficiency assessment studies because they do not require assumptions about the functional relationships between variables and avoid too much subjectivity. However, the traditional radial DEA model is also partially flawed in that it does not consider slack variables in the efficiency measure of inefficient DMUs. In traditional radial DEA models (e.g., CCR-DEA, BCC-DEA, etc.), the assessment of the degree of DMU inefficiency only includes the proportional change of outputs and inputs. For the ineffective DMU, the gap between its current state and the state of the effective DMU on the production frontier should contain both equal proportional improvement and slack improvement, and the slack improvement part is not reflected in the traditional radial DEA model. Especially when the number of input-output indicators is relatively large, there will be more invalid DMUs, and the traditional radial DEA model loses the function of measuring the slack variables of invalid DMUs, and in such cases, the measurement results are not accurate to a certain extent. The slack variable-based efficiency measurement model (SBM) proposed by Tone [46] incorporated slack variables of input-output indicators into the calculation of decision unit efficiency, which effectively improved the problem of lack of consideration of slack variables in traditional DEA models.
Assuming there are n DMUs, each DMU contains m input X and q output Y. The a-th input of the k-th DMU is expressed as x ak (a = 1, 2, ..., m; k = 1, 2, . . . , n), the b-th output is expressed as y bk (b = 1, 2, . . . , q; k = 1, 2, . . . , n), then the input-oriented SBM-DEA model with constant returns to scale can be expressed as [46]: In Equation (1), X k = (x 1k , x 2k , . . . , x mk ), Y k = y 1k , y 2k , . . . , y qk ; s − is the slack variable, i.e., the slack improvement value of the input index; θ and λ are the relative efficiency value and weight of the decision unit respectively, if θ ≥ 1, it stands that the DMU is in production frontier and belongs to the effective DMU.

Window-DEA Model
In most cases, efficiency values need to be dynamically assessed for multiple regions in different years. However, the SBM-DEA model cannot measure the whole panel data directly because the production frontier is different in different years, and the model can only decompose the panel data into cross-sectional data for static measurement separately. Therefore, the efficiency results measured by SBM-DEA model alone are not comparable between years. However, the Window-DEA model, also known as the DEA window analysis method, is a good solution to these problems [47].
The basic principle of DEA window analysis method is that the same decision unit in different years is regarded as multiple independent DMUs to participate in the calculation. Then, a number of reference sets are constructed based on the moving average method to realize the dynamic evaluation of efficiency values under the conditions of multiple decision units and long time series. Through the window analysis, it can meet the need of dynamic comparative analysis of each decision unit in time series. Moreover, the same cross-sectional data are repeatedly involved in the calculation, which can more fully explore the data value and reflect the real level. Thus, the DEA window analysis method is fully applied in the dynamic assessment of efficiency or performance of long time series in many fields. The difference between the Window-DEA model and the SBM-DEA model is only the change of reference set thus the DEA window analysis method can be well combined with the SBM-DEA model.

Malmquist Index Model
In order to analyze the evolution of CWUI in more depth and to explore the deepseated causes of efficiency changes, this study uses the Malmquist index model to explore the changes of TFP of water. Malmquist index is a non-parametric indicator used to dynamically measure the change in TFP of each DMU [48]. Moreover, Färe et al. [49] firstly used the Malmquist index for efficiency research in conjunction with DEA theory and proposed a decomposition method for total factor productivity change factors. t + 1 period compared to t period productivity change can be decomposed into integrated technical efficiency change (EC) and technical change (TC). TC refers to technological progress, which can also be understood as the forward movement of the production frontier. EC and TC are the two main causes of changes in TFP. Moreover, EC can be decomposed into deeper pure technical efficiency change (PEC) and scale efficiency change (SEC). The decomposition method can be expressed as [49]: In Equation (2), MI is the Malmquist index, MI t,t+1 > 1 means that the TFP of water in period t + 1 has increased compared to period t; EC t,t+1 > 1 means that the DMU is closer to the production frontier in period t + 1 compared to in period t, representing an increase in overall technical efficiency; and TC t,t+1 > 1 means that the production frontier in period t + 1 has moved forward compared to period t, representing technological progress.
In this study, there is a strong link between the SBM-DEA model and the Malmquist index model. The SBM-DEA model measures integrated technical efficiency, and the Malmquist index model measures the rate of change in total factor productivity. The rate of change in TFP is not TFP itself. The Malmquist index model is a nonparametric model that cannot calculate TFP itself but can calculate the rate of change in TFP. The rate of change in TFP is also known as the Malmquist index (MI), and it can be further decomposed into the rate of change of integrated technical efficiency EC and technical progress TC. In efficiency studies involving long time series, the production frontier is constantly moving forward and can be understood as technical progress TC. The analysis of technical efficiency using only the SBM-DEA model is incomplete. Therefore, this study uses the Malmquist index model as an effective complement to the SBM-DEA model, and the two models together constitute the quantitative assessment model of water use level in this study.

Spatial Matching Degree Calculation Method Based on Series Distance
The matching degree calculation method based on the distance of the series proposed by Zuo et al. [50] can quantitatively describe the matching degree between different variables. The basic principle of this method is to create a new data series by using the proportion of a value in the series to the total value of the series and then characterize the match according to the ratio of the difference between the values of different variables and the maximum distance in the new data series. In this paper, the methodology is used to quantitatively assess the matching relationship between CWUI and economic and social development level (E-SDL) in nine provinces of the Yellow River Basin. Assuming that the study object contains a total of N subunits, a quantitative study of the match between two variables, A 1 and A 2 is carried out using a spatial matching degree calculation method based on the series distance, which is calculated as follows: In Equation (3), , (r = 1, 2, · · · , N); MD r is the degree of the match between the two variables A 1 and A 2 in the rth subunits.
It is important to note that Equation (3) only applies to the quantification of the degree of the match when the two variables are positively correlated (i.e., the larger the value of A 1 and the larger the value of A 2 , the better the match between the two variables). The formula for calculating the degree of the match when the two variables are negatively correlated can be given similarly (i.e., the larger the value of A 1 and the smaller the value of A 2 , the better the match between the two variables): In Equation (4), the meaning of each symbol is the same as before.

Study Area
The Yellow River is the longest river in north China. The Yellow River basin is an extremely important economic and cultural corridor and natural ecological barrier in China. It contains several basic energy bases and key ecological function areas, and it is also one of the basic water supply sources for northern China, playing an important role in the construction of China's ecological civilization and economic and social development. From the source to the mouth, the Yellow River passes through Qinghai, Sichuan, Gansu, Ningxia, Inner Mongolia, Shaanxi, Shanxi, Henan, and Shandong in turn, called the "nine provinces of the Yellow River basin," hereinafter referred to as the "nine provinces." Most of the nine provinces are arid or semi-arid areas, with a large geographical span and significant differences in altitude. There are large differences in resource and environmental endowments and economic development levels in different provinces and regions, and the difference in per capita GDP between the source and the mouth of the Yellow River is more than 10 times [51]. The elevation distribution of nine provinces of the Yellow River basin is shown in Figure 1. The current situation of economic development of the nine provinces is lower than the overall level of the whole China. According to the "2019 Monitoring Report on China's Cities Completely Built a Well-off Society" issued by the Chinese government, only one of the top 20 prefecture-level cities in the Overall Well-off Index (which can represent the level of economic and social development to some extent) is located in the nine provinces. The nine provinces of the Yellow River Basin account for nearly one-third of the population and 34.8% of the arable land in China. However, the annual runoff of this river accounts for only 2% of China's total, with serious water supply and demand conflicts constraining regional development [52]. Especially in the past 20 years, under the general trend of highly rapid economic and social development in China, the water demand of the nine provinces has also increased dramatically. The Yellow River Basin is facing a severe situation of contradiction between water supply and demand. Shrinking natural incoming water and low levels of water utilization are also the causes of water shortages in the basin. The comprehensive improvement of CWUI has significance to the sustainable development of the nine provinces, especially to these provinces in the upper reach.

Indicator Description and Data Sources
When applying the DEA model, the reasonableness of the input-output indicators will directly affect the accuracy of the final efficiency measurement results. Thus, it is very important to construct a comprehensive and effective set of CWUI input-output indicators. The basic idea of efficient utilization of water is to satisfy both minimizations of resource consumption and maximization of production value. Guided by this idea, this study selects four representative indicators from four dimensions of resources, capital, labor, and ecological environment as input variables. The most direct input to reflect the level of water use is water consumption, and the total regional water consumption is used to represent water consumption [27,30,31,34]. The level of capital investment is also an important aspect in quantifying the level of resource utilization, using fixed-asset investment to represent social capital investment [27,[29][30][31]34]. This study discusses the overall level of regional water utilization. Considering that all industries are inseparable from water utilization, the number of employees is taken as a variable reflecting human resource input [27,[29][30][31]. The water use level without considering environmental pollution cannot reflect the real level of regional water utilization [30,34]. In order to evaluate the CWUI of nine provinces in the Yellow River basin more scientifically and comprehensively, this study adds ammonia nitrogen emission in wastewater as an input indicator reflecting environmental pollution. The most intuitive manifestation of a high-level area of water utilization is the ability to obtain higher economic benefits under certain conditions of resource inputs. Therefore, the gross regional product was selected as the output variable reflecting the value of water resources production [27,[29][30][31]34]. The final input-output indicator selection results are listed in Table 1. Table 1. Water use level input-output indicator system.

Type
Indicators Necessity

Input indicators
Total water consumption Reflect natural resource input [27,30,31,34] Investment in fixed assets a Reflect social capital input [27,[29][30][31]34] Number of employees Reflect human capital input [27,[29][30][31] Ammonia nitrogen emissions in wastewater Reflect environmental load-bearing input [30,34] Output indicators GDP b Reflect economic benefit output [27,[29][30][31]34] Notes: a Investment in fixed assets data is processed through the fixed asset investment price index to eliminate price effects. b GDP data are processed through the GDP index to eliminate price effects.
High-quality economic and social development involves economic, social, and ecological aspects. Thus, the level of regional economic and social development needs to be characterized by a combination of indicators. In order to investigate the matching relationship between CWUI and E-SDL in the nine provinces, this study refers to relevant research results [14,50,51], as well as the "Statistical Monitoring Program for Building a Well-off Society" and the "Statistical Monitoring Indicator System for Building a Well-off Society in All Respects" issued by the Chinese government, and follows the principles of representativeness and dynamism to select a total of 12 representative indicators to characterize the relative level of economic and social development in terms of economic development, social harmony, and ecological friendliness. The indicator system is shown in Table 2. Notes: a The "+" in the indicator characteristic represents the "positive indicator," and the larger the value, the stronger the positive effect on E-SDL, and the "-" represents the "negative indicator," and the larger the value, the stronger the negative effect on E-SDL. b The weight of each indicator is combined with the corresponding indicator weights in the "Statistical Monitoring Program for Building a Well-off Society" to be converted to determine.
In order to ensure the reliability of data and the integrity of the time series, a total of seven research years from 2012 to 2018 were finally selected in this study to carry out the case study of nine provinces. The raw data used in this study were obtained by collating data from the China Statistical Yearbook, the China Science and Technology Statistical Yearbook, the Water Resources Bulletin of nine provinces, and the Statistical Yearbook of nine provinces.

Spatial-Temporal Evolution Characteristics of Water Use Level
In this study, the values of CWUI of nine provinces during the period 2012-2018 needed to be dynamically evaluated. When applying the DEA window analysis method, there was no technical progress within a single window because all decision units within each window were involved in the calculation [53]. To increase the credibility of the results, DEA window analysis should be applied with as narrow a window width as possible [54]. Therefore, the optimal window width for the Window-DEA model was determined to be 3 years. The total length of the time series in this study was 7 years, thus for each decision unit, five windows needed to be established, namely Window1   As shown in Figure 2, there is an overall decreasing trend in the values for different windows within a given year during the study period. Based on the data distribution characteristics of CWUI, CWUI was classified into five levels of (0-0.5), (0.5-0.7), (0.7-0.8), (0.8-0.9), and (0.9-1.0), representing poor, medium, good, outstanding, and excellent water use levels, respectively. Taking Qinghai province as an example, in 2015, the value under Window2 was 0.506, CWUI was at a medium level. However, the value under Window3 dropped to 0.425, and under Window4 its value dropped to 0.394, under both windows, CWUI was at a poor level. Compared to Window3 and Window4, the CWUI under Window2 was at a relatively good level, but under Window4, which had a better production frontier, there was a significant decrease in its CWUI value. The CWUI values measured by the data envelopment analysis model were not absolute but were derived by comparing the decision-making unit with the production frontier. The CWUI values of the same decision-making unit in the same year can vary depending on the production frontier, thus, the CWUI value has relativity. On the one hand, it shows that the production frontier surface of the nine provinces is moving forward, and the level of water use in the advanced provinces is improving rapidly; on the other hand, it also confirms the relative nature of CWUI values. The results of the window analysis of the SBM-DEA model in Table 3 were compiled to obtain the final results of the CWUI values of the nine provinces from 2012 to 2018, as shown in Figure 2f. The multi-year average values and rankings of CWUI for the nine provinces are listed in Table 3. From an overall perspective: the overall CWUI values of the nine provinces were above 0.55 from 2012 to 2018, showing a trend of first decreasing and then increasing. The inflection point occurred in 2015, indicating that the overall CWUI level of the nine provinces deviated to the greatest extent from the production frontier surface in 2015, after which the degree of deviation gradually decreased. Subsequently, the CWUI value began to gradually rebound, rising significantly to 0.704 in the three years of 2016-2018, reaching the highest value in the study period.
From the perspective of each province: among the nine provinces, Shandong, Inner Mongolia, and Shanxi had a relatively high level of water utilization, with a multi-year average CWUI above 0.7, at a good level and above. Shandong was the only province that water use was at an excellent level and at a significant advantage in water use level compared to other provinces. Followed by Henan and Shaanxi provinces, with multi-year average CWUI above 0.6 and small fluctuations in CWUI, both provinces were at the medium level. The multi-year average CWUI of Sichuan, Gansu, Qinghai, and Ningxia were all below 0.5, which deviated from the production frontier surface to a large extent, and the water use level was at a poor level among the nine provinces. During the period 2012-2017, only Shandong province reached the production front surface among the nine provinces, and its CWUI was 1 in 2013, 2014, and 2015, indicating it was the only effective DMU among the nine provinces for three consecutive years. In 2018, the CWUI of Inner Mongolia and Shanxi also exceeded 1 and were jointly located on the production frontier surface with Shandong, while the remaining six provinces failed to reach the production frontier surface during 2012-2018, and the CWUI values of some provinces showed a decreasing trend in multiple years, which to some extent indicates while that the spatial variability of CWUI in the nine provinces of the Yellow River Basin was obvious, the synergy of water use level improvement was poor, and the gap between high and low CWUI provinces is gradually expanding, which is not conducive to the balanced and coordinated development of the nine provinces.
From the perspective of different regions: according to the dividing points of the upper, middle, and lower reaches of the Yellow River, the provinces of Inner Mongolia and above were divided into upstream areas, Shaanxi and Shanxi were divided into midstream areas. In order to facilitate the analysis of regional differences in water use levels in the Yellow River basin, Henan and Shandong were divided into downstream areas in this paper. Finally, the temporal evolution trends of CWUI in different regions are shown in Figure 3. The CWUI evolution trends of upstream and downstream provinces were basically the same, while the CWUI evolution trends of midstream provinces had more obvious fluctuations, and their water utilization levels were significantly improved in multiple years. The spatial distribution characteristics of CWUI in upstream, middle, and downstream provinces were basically consistent with the distribution of their economic and social development levels. As can be seen from Figure 3, the CWUI of middle and downstream provinces was higher than that of upstream provinces in all years. Moreover, in 2014, the difference between the CWUI of upstream and downstream provinces reached 0.32, indicating that the industrial structure of upstream provinces needs to be optimized, there is a large redundancy of resource inputs. There is still a gap between the level of intensive use of water in the upstream provinces and the middle and downstream provinces.

TFP Analysis of Water Based on Malmquist Index Model
Based on the panel data of input-output variables of nine provinces from 2012 to 2018, the Malmquist index model was applied to measure the productivity index of nine provinces, and the total factor productivity index and its decomposition result of different years and provinces were obtained, as shown in Table 4 and Figure 4. The operation was implemented by DEAP2.1 software (Université Laval, Quebec City, Canada) [56]. From the results in Table 4, the multi-year average value of MI was 1.078, and MI greater than 1 indicates that the overall TFP of water of the nine provinces showed an increasing trend during the period of 2012-2018. Among them, the overall total factor productivity of water in the nine provinces declined year by year before 2016, but the declining trend gradually leveled off. Moreover, the overall total factor productivity of water has increased year by year since 2016. In order to explore the intrinsic causes of TFP changes, the results of total factor productivity decomposition were further analyzed. Moreover, the multi-year average values of all indices were less than 1, except for technical progress TC, which indicates that technical progress TC was the dominant factor driving TFP improvement. EC was greater than 1 only in 2016 and 2017, indicating that changes in integrated technical efficiency only played a driving role in individual years. However, EC, in general, constrains the improvement of total factor productivity of water thus the industrial structure and water resources management mode need to be optimized to improve the scale efficiency and pure technical efficiency of water utilization in the nine provinces. As can be seen from Figure 4, the MI of Inner Mongolia, Shanxi, Shaanxi, and Shandong provinces are all higher than the nine-province average of 1.078, indicating that the water use level of these four provinces is improving faster. Among them, Inner Mongolia and Shandong provinces, because they are at the production frontier themselves, the improvement of TFP of water is entirely driven by the technological progress TC. In comparison, Shanxi and Shaanxi are driven by both technological progress TC and integrated technical efficiency change EC. The remaining five provinces have been slow to improve their TFP of water, mainly due to the constraints of the integrated technical efficiency change EC, and the distance between their water utilization levels and the production frontier surface has gradually increased. Further decomposition of EC can be found that Qinghai, Gansu, and Ningxia are mainly affected by the scale efficiency change SEC, and the scale efficiency needs to be further optimized, while Sichuan and Henan are mainly constrained by the pure technical efficiency change PEC, and need to focus on strengthening the water resources management level and optimizing the water use structure.
Comparing the Malmquist index and its decomposition results from the perspective of different regions. The comparison of upstream, midstream, and downstream provinces is shown in Figure 5. The overall productivity level of the midstream provinces is higher than that of the upstream and downstream provinces. It is not difficult to find that PEC and SEC of the upstream provinces are both at a lower level, especially the poor performance of SEC, which in turn leads to a much lower integrated technical efficiency change EC than that of the midstream provinces. Thus, the upstream provinces need to reasonably allocate water resources, optimize industrial layout, improve industrial concentration, and narrow the gap between their SEC and that of the midstream and downstream provinces. Technological progress TC is the main driving force of TFP of water in the three regions, and the difference of TC in the three regions is relatively small, with the upstream provinces even slightly higher than the downstream provinces. The production frontier of each of the three regions has moved forward significantly, indicating that the nine provinces have all invested more in scientific research in the field of efficient water utilization in recent years, effectively promoting the transformation of advanced technology into productivity [57].

Analysis of the Relative Level of Economic-Social Development
Based on the quantitative index system of the relative level of economic and social development in the nine provinces in Table 2, the data of each index in the nine provinces during the period of 2012-2018 were normalized, and then all the indexes were weighted and summed up by combining the weights to obtain the quantitative characterization value of E-SDL. Arcgis10.2 software was used to map the E-SDL calculation results of the nine provinces for multiple years into a spatial distribution map, which was used to characterize the spatial and temporal distribution characteristics of the economic and social development levels of the nine provinces, as shown in Figure 6. On the time scale, the average value of E-SDL in the nine provinces showed an increasing yearly trend during the seven years, from 0.352 in 2012 to 0.657 in 2018, indicating that with the gradual implementation of the concept of sustainable and green development, the overall level of economic and social development in the nine provinces is continuously improving. Gansu, Qinghai, Sichuan, and Henan, have seen the fastest improvement in economic and social development, with E-SDL more than doubling in seven years. From the three criteria, it was found that the economic development and social harmony of the above four provinces rapidly improved in the past seven years under the premise of maintaining the steady growth of the level of ecological friendliness, thus realizing the obvious improvement of E-SDL. On the spatial scale, there were large spatial differences in the economic and social development levels of the nine provinces. In terms of the multi-year average value of E-SDL, the economic and social development level of Shandong was relatively high (E-SDL > 0.7), the E-SDLs of Inner Mongolia, Shaanxi, Henan, Sichuan, Shanxi, and Ningxia were at an intermediate level (0.4 < E-SDL < 0.6), and the economic and social development levels of Gansu and Qinghai were low (E-SDL < 0.4). The E-SDL of the nine provinces as a whole shows the spatial distribution characteristics of downstream, midstream, and upstream provinces in descending order, which is basically consistent with the regional step structure of economic and social development in the Yellow River basin [14].

Analysis of Matching Degree between CWUI and E-SDL
The spatial matching degree of the multi-year average CWUI and E-SDL of nine provinces was calculated using the spatial matching degree calculation method based on the series distance, and the results are shown in Figure 7. To further analyze the temporal change characteristics of the matching degree, the spatial matching degree of CWUI and E-SDL of the nine provinces and the upstream, midstream, and downstream provinces of the Yellow River in 2012, 2015, and 2018 were measured, and the results are listed in Table 5.  From the results in Figure 7, the CWUI of Gansu, Sichuan, Henan, Qinghai, and Shandong matches well with the E-SDL (MD > 0.9), but the matching categories in different provinces are not exactly the same. Shandong and Henan provinces have a high level of economic and social development, and their water use level is relatively good among the nine provinces, which is a "high-high" match. Gansu and Qinghai provinces have a relatively low level of economic and social development, and their CWUI is also at the back of the nine provinces, which is a "low-low" match. The degree of matching CWUI with the E-SDL in Inner Mongolia and Shaanxi provinces is medium among the nine provinces (0.8 < MD < 0.9). The match between CWUI and E-SDL in Ningxia and Shanxi provinces is poor (MD < 0.8), and the reason for this is that Shanxi Province has a high CWUI among the nine provinces, while the level of economic and social development is relatively low, resulting in the two not matching in value. Ningxia, however, is at the bottom of the nine provinces in terms of water use level, but its economic and social development is at the middle level among the nine provinces, and it needs to focus on improving water resources utilization in the future economic and social development process. In another study on the coordination between water use and urbanization [58], the results show that there are significant differences in the synergistic effects between water use level and urbanization in different provinces in China. Moreover, the differences can be divided into three different types of synergy, with only individual provinces showing a significant synergistic effect of urbanization level on water use. The results are similar to those of this paper, and to a certain extent, they verify the spatial variability of CWUI and E-SDL among nine provinces in the Yellow River basin.
From the results in Table 5, the matching degree of CWUI and E-SDL in Gansu Province has been increasing within the three study years of 2012, 2015, and 2018, and its MD has increased from 0.885 in 2012 to 0.990 in 2018. The matching degree of water use level and economic and social development level in Gansu Province in 2018 is the highest among the nine provinces, but it belongs to the "low-low" matching type, which still needs to improve CWUI and E-SDL comprehensively. The match between CWUI and E-SDL in Ningxia and Inner Mongolia provinces has been decreasing, with MD in Ningxia decreasing from 0.761 in 2012 to 0.729 in 2018 and MD in Inner Mongolia decreasing from 0.889 in 2012 to 0.602 in 2018. It indicates that compared with other provinces, the synergy between the level of water utilization and the improvement of economic and social development level in these two provinces is poor, and should focus on the coordinated and harmonious development of water resources utilization and economy and society. The main problem in Ningxia is the low level of water utilization, and the level of water utilization in Inner Mongolia is at a high level, but its economic and social development level is not outstanding, and it needs to pay attention to the all-round balanced development of economy, society, resources, and ecology.
The MD of the six provinces of Qinghai, Sichuan, Henan, Shaanxi, Shandong, and Shanxi showed the characteristics of first rising and then falling within three years. CWUI and E-SDL reached the best matching status in 2015. Among the six provinces mentioned above, the characteristics of MD changes can be explained in three cases. The first case is that the improvement of water use level lags behind the improvement of economic and social development level, such as Sichuan, Shaanxi, and Qinghai provinces. Sichuan's CWUI was ranked 6th in all three years 2012, 2015, and 2018, but its E-SDL improved significantly from 7th in 2012 to 4th in 2018. Shaanxi's CWUI ranked 5th in 2012 and 2018 and 4th in 2015, but its economic and social development level improved from 3rd in 2012 to 2nd in 2018. Qinghai's CWUI dropped from the 7th in 2012 to the 9th in 2018, and its E-SDL remained at the 8th for three years. The second case is that the E-SDL is rising faster, but the level of water utilization is rising slower, such as Shanxi Province. Its CWUI rose from 4th in 2012 to 2nd in 2018, but its E-SDL dropped from 5th in 2012 to 7th in 2018. The third case is that the level of water utilization and the level of economic and social development has maintained simultaneous improvement to some extent, such as Henan and Shandong. Henan's CWUI decreased from 3rd in 2012 to 4th in 2018, and its E-SDL decreased from 4th in 2012 to 5th in 2018. CWUI and E-SDL of Shandong Province have maintained 1st place in 2012, 2015, and 2018, and although there are fluctuations in MD, the fluctuations are not obvious, and the water use level and economic and social development level have maintained a relatively coordinated development. Analyzed from the perspective of different regions in the Yellow River Basin, the overall situation is that the MD of downstream provinces is higher than that of upstream provinces than that of midstream provinces. Except for the downstream provinces, the MDs of the remaining two regions showed a decreasing trend year by year during 2012, 2015, and 2018, and their economic and social development and water use level failed to achieve a harmonious and balanced improvement, and the upstream and midstream provinces of the Yellow River Basin still need to focus on the simultaneous improvement of water use level under the premise of ensuring economic and social development.

Conclusions
Based on the SBM-DEA model and combined with the Window-DEA model, this study measured the water use level of nine provinces in the Yellow River Basin from 2012 to 2018 and decomposed the TFP changes of water in the nine provinces using the Malmquist index model. A quantitative index system of economic and social development levels was constructed, the spatial and temporal variation characteristics of the relative levels of economic and social development in the nine provinces were analyzed, and finally the matching relationship between water use level and economic and social development levels in the nine provinces of the Yellow River Basin was explored. The results were concluded as follows: (1) The nine provinces showed an overall increasing trend in water use level during 2012-2018, with the midstream and downstream provinces having higher overall water use levels than the upstream provinces. The water use level of Shandong Province was at the production frontier of the nine provinces for a long time, and the gap between the water use levels of different provinces is gradually expanding. The upstream and midstream provinces of the Yellow River Basin should optimize the water resources utilization pattern, combine the advanced experience of water resources utilization in Shandong Province, and explore the effective ways to improve the level of regional water resources utilization according to local conditions. Efforts should be made to narrow the gap between water use level and that of downstream provinces. (2) Technological progress TC is the dominant factor driving the TFP of water in the nine provinces, and integrated technical efficiency change EC overall constrains the TFP of water. Insufficient SEC and PEC are the main reasons for the productivity level gap between upstream and downstream provinces and midstream provinces. The relevant departments of water resources management in the Yellow River Basin should strengthen the implementation of the most stringent water resources management system. While ensuring technological progress, relevant departments should strengthen the level of water resources management, promote the adjustment of its industrial structure in areas lagging behind in water utilization, and force the improvement of comprehensive technical efficiency of regional water resources utilization. (3) There are significant regional differences in the match between the level of economic and social development and the level of water utilization in the nine provinces, with Gansu, Sichuan, Henan, Qinghai, and Shandong having a better match, but with different match types. Among them, Gansu and Qinghai belong to the "low-low" matching type, which needs to improve both the level of water utilization and the level of economic and social development. The matching degree of CWUI and E-SDL in Ningxia and Shanxi is poor, and there is no dynamic synergy between the improvement of economic and social development level and the improvement of water use level. The harmonious relationship between resource utilization and economic and social development is an important foundation for regional sustainability. Moreover, the nine provinces in the Yellow River Basin should explore new paths for high-quality development based on a harmonious balance between the two.
The findings of the study clarify the characteristics of temporal changes and spatial distribution patterns of water use levels in the nine provinces of the Yellow River Basin in recent years. The key factors leading to the variation of total factor productivity of water resources in the nine provinces are also identified, and the results of the study can provide some basis for the strictest water resources management in different regions of the Yellow River Basin. In addition, this paper analyzes the adaptation characteristics between water resources utilization and economic and social development in the nine provinces, which can provide some reference for the future high-quality development layout of the Yellow River basin to a certain extent. It is also a useful exploration to carry out more complex research on the relationship between water resources utilization and economic and social development in the future.
Despite that, limitations also exist in the present study, which could be further improved. First, the SBM-DEA model cannot allow further comparison of the CWUI values of multiple effective DMUs. For example, in 2018, we learned that the CWUI of Inner Mongolia, Shanxi, and Shandong are all one, but we cannot make further comparisons of the three levels. A possible solution to this problem is to explore the applicability of more comprehensive DEA models (e.g., super-efficient SBM-DEA models) in water use level studies. Second, the water use level involves many aspects such as resources, ecology, economy, and society, etc. Although the input-output index system of water use level constructed in this study covers as many representative indicators as possible, the extent to which it can represent the real level of water utilization still needs to be further explored. Third, this study explores the spatial matching characteristics between CWUI and E-SDL using the method of calculating spatial matching based on the series distance. However, this method cannot explore the degree of coordination of different variables on the time series. In future research, a combination of multiple coordination relationship exploration methods (e.g., Tapio decoupling model, coupling coordination degree model) can be used to comprehensively explore the adaptation relationship between CWUI and E-SDL in both time and space dimensions.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: [http://www.stats.gov.cn/tjsj/ndsj/].