Identification and expression analysis of GARP superfamily genes in response to nitrogen and phosphorus stress in Spirodela polyrhiza

GARP transcription factors perform critical roles in plant development and response to environmental stimulus, especially in the phosphorus (P) and nitrogen (N) sensing and uptake. Spirodela polyrhiza (giant duckweed) is widely used for phytoremediation and biomass production due to its rapid growth and efficient N and P removal capacities. However, there has not yet been a comprehensive analysis of the GRAP gene family in S. polyrhiza. We conducted a comprehensive study of GRAP superfamily genes in S. polyrhiza. First, we investigated 35 SpGARP genes which have been classified into three groups based on their gene structures, conserved motifs, and phylogenetic relationship. Then, we identified the duplication events, performed the synteny analysis, and calculated the Ka/Ks ratio in these SpGARP genes. The regulatory and co-expression networks of SpGARPs were further constructed using cis-acting element analysis and weighted correlation network analysis (WGCNA). Finally, the expression pattern of SpGARP genes were analyzed using RNA-seq data and qRT-PCR, and several NIGT1 transcription factors were found to be involved in both N and P starvation responses. The study provides insight into the evolution and function of GARP superfamily in S. polyrhiza, and lays the foundation for the further functional verification of SpGARP genes.

Background GARP (G: GOLDEN2 in Zea mays, AR: ARR-B in Arabidopsis thaliana, P: Psr1 in Chlamydomonas reinhardtii) is a plant specific transcription factor (TF) superfamily with diverse functions involved in nutrients (N and P) sensing, chloroplast development, circadian clock oscillation, and hormonal signaling [1,2]. GARP superfamily TFs are made of type-B authentic response regulators (ARR-B) and GOLDEN2-like (G2-like, GLK) TFs, and first be named by Riechmann in A. thaliana [1,2]. Both of these two TFs contain GARP motif (also named B-motif ), which is a domain of about 60 amino acids that forms a three α-helices 3D structure containing a helix-turn-helix (HTH) motif [2,3]. B-motif is the multifunctional domain in GARP proteins that responsible for both nuclear localization and DNA binding [2,3].

Open Access
*Correspondence: houhw@ihb.ac.cn 1 The State Key Laboratory of Freshwater Ecology and Biotechnology, The Key Laboratory of Aquatic Biodiversity and Conservation of Chinese Academy of Sciences, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan 430072, China Full list of author information is available at the end of the article GOLDEN2 gene was first reported to function in chloroplast biogenesis in Z. mays [4]. Then, the orthologues of GOLDEN2 were identified in A. thaliana, Physcomitrella patens, Capsicum annuum, Oryza sativa, and Solanum lycopersicum [4][5][6][7][8]. AtGLK1/2 in A. thaliana are the orthologous genes of ZmGOLDEN2, expressing in photosynthetic tissue and exhibit functional redundancy, however, the double mutants reduce the accumulation of photosynthetic gene products and thylakoids in chloroplasts [7]. ZmGOLDEN2 and their orthologues in O. sativa and A. thaliana are regulated by light [6,7]. Subsequently, GOLDEN2 and its orthologues were classified into the GARP family with ARR-B and Psr1-like genes [1]. A gene can be classified in the GARP family if the derived protein contain B-motif [2]. ARR-B proteins contain a B-motif in their C-terminal, act as positive regulators in the two-component cytokinin signaling pathway and play pivotal roles in plant development, including vascular development, light sensitivity, chlorophyll production, hypocotyl elongation, and cell division in the root and shoot [9]. GLK TFs play important roles in biotic and abiotic stresses, such as pathogen infection, salt stress, drought, and nutrient deficiency. AtGLK1 regulates genes that involved in disease resistance and effect on different pathogen [10], nine AtGLK genes were responded to salt stress in A. thaliana [11]. Both of Phosphorus Starvation Response 1 (Psr1) in C. reinhardtii and Phosphate Starvation Response 1 (PHR1) in A. thaliana play the central regulatory roles in inorganic phosphate (Pi) sensing system [12,13].
The duckweed family (Lemnaceae) is a group of fastgrowing free-floating aquatic plants, and distributed in various fresh water environments throughout the world besides the most extreme habitats [14]. Based on the morphological characteristics and molecular taxonomy, 36 duckweed species were recognized belonging to five genera: Spirodela (2), Landoltia (1), Lemna (12), Wolffiella (10), and Wolffia (11) [15,16]. Duckweed is an attractive model in plant research for the convenient cultivation system, clear genetic background, and robust transformation methods [17]. They are some of the fastest growing flowering plants (doubling time < 30 h under the optimal growth conditions) in the world, with high productivity of dry mass (80-100 tons per hectare per year) which is more than five times that of maize [18][19][20]. Duckweed also exhibits efficient N and P removal capacities from wastewater with about 1.3 g•m −2 •d −1 and 0.18 g•m −2 •d −1 respectively [21]. Thus, duckweed is also considered as ideal plant in phytoremediation to recover nutrients (N and P) from eutrophic water. S. polyrhiza occupies the ancestral phylogenetic position among duckweeds, possesses the largest individual and the smallest genome size in the Lemnaceae [22]. The prominent performance of S. polyrhiza in nutrient removal from wastewater has been observed in previous studies [23,24]. The available genomic data and robust transformation method of S. polyrhiza provide theoretical and technical supports for the research of molecular mechanism and germplasm improvement in S. polyrhiza [25][26][27][28][29][30]. Thus, S. polyrhiza is one of the best prospective species in Lemnaceae for phytoremediation and biomass production.
N and P are two major essential elements required for plant growth. Nitrate (NO 3 − ) and ammonium (NH 4 + ) are the main forms of N source in soils, and inorganic phosphate (Pi) is the main form of P source for plant uptake in the environment [31,32]. The distributions of N and P in the environment are affected by many factors which create a variable spatiotemporal landscape at the local and global scale [33]. Thus, it is crucial for plants to adapt themself to various nutrient environments (N and P content), however, the underlying molecular mechanism remains to be elucidated in-depth. The GARPs TFs are involved in the responses to nutrients and include probable nutrient sensors of plants. CrPsr1 is the first reported GARP TF involving in nutritional responses, it is critical for the acclimation of C. reinhardtii to P starvation [12]. AtPHR1 was the central regulator in the downstream of Pi starvation signaling pathway [13,34]. Then, Nitrate-Inducible Garp-Type Transcriptional Repressor 1 (NIGT1)/Hypersensitivity to Low Phosphate-Elicited Primary Root Shortening 1 (HRS1)/HRS1 Homolog (HHO) were found to be the most robustly and quickly NO 3 − regulated genes [35][36][37]. Furthermore, NIGT1 TFs are involved in both N and P sensing, uptake, and assimilation through the interactions with PHR1, Nodule Inception (NIN)-like protein (NLP), SYG1-Pho81-XPR1 (SPX) domain proteins, phosphate transporter 1 (PHT1), nitrate transporter 1 (NTR1), and NRT2 [38,39].
The GARP superfamily or GLK subfamily has been identified and characterized in A. thaliana, O. sativa, Z. mays, Nicotiana tabacum, and Gossypium hirsutum [1,[40][41][42][43]. Bhutia et al. (2020) found that some GARP members (OsGLK10, OsGLK15, OsGLK22, and OsGLK30) respond to P starvation in rice [41]. However, the GARP gene family has not been thoroughly examined in S. polyrhiza, to the best of our knowledge. In current study, we perform a comprehensive analysis of GARP superfamily in the giant duckweed S. polyrhiza, including chromosomal locations, evolutionary perspectives, structural arrangement, and their functional role through gene expression analysis. The results of this study offer a robust platform for further functional studies of the candidate SpGARP genes, so as to understand the molecular mechanism of N/P sensing, acquisition and balance, and enable their efficient use in the germplasm improvement to enhance the phytoremediation potential of S. polyrhiza.

Genome-wide identification of GARP superfamily genes in S. polyrhiza
The GARP proteins in Arabidopsis and rice were used as queries to perform the BLASTP search in S. polyrhiza 7498 genome databases. To mine all the potential GARP proteins that harbored B-motif in S. polyrhiza, the Hidden Markov Model (HMM) profile of B-motif was built using HMMER 3.3.2 based on the identified SpGARPs [44], and used for further search of SpGARP proteins. A total of 35 GARP genes were identified in S. polyrhiza, including 7 SpARR-B and 28 SpGLK genes. The detailed information of SpGARP genes was listed in Table 1, including gene ID, genomic position, gene length, exon number, protein length, molecular weight (MW), isoelectric point (pI), instability index, grand average of hydropathicity (GRAVY) and subcellular localization. The GARP genes were named as SpARR-B1-7  (Table S1).

Phylogenetic analysis of GARP superfamily genes in five species
The full-length amino acid sequences of SpGARPs, AtGARPs (14 AtARR-B and 41 AtGLK genes), OsGARPs (10 OsARR-B and 47 OsGLK genes), CeGARPs and WaGARPs were subjected to multiple sequence alignments using ClustalW [45]. The unrooted neighborjoining (NJ) phylogenetic tree was divided into eight branches belonging to three big groups (I, IIa-f and III). As shown in Fig. 1, the majority of ARR-B proteins were clustered into the group I, besides that, six members (At4G12020.1, CeARR-B8, OsGLK27, SpARR-B7, WaARR-B3, and WaARR-B4) were clustered into clade IIa. The orthologous genes of ZmGOLDEN2 in five species were clustered into the subgroup IIb. Both of Arabidopsis and rice possess two members, C. esculenta and S. polyrhiza have a single copy in their genome, and three orthologous genes were present in the smallest and simplest flowering plant W. australiana. The findings indicate the unusual evolutionary pathways of ZmGOLDEN2 orthologous genes in Araceae plants, especially in the duckweed plants. A total of 20 GLK genes (5 AtGLKs, 5 CeGLKs, 4 OsGLKs, 2 SpGLKs and 4 WaGLKs) were clustered into subgroup IIc. The members in subgroup IIc were involved in circadian oscillation, includ- genes in Arabidopsis and rice [41]. The members of 25 NIGT1/HRS1/HHO (also named NIGT1) TFs were clustered into subgroup IId. It has been found that NIGT1 TFs coordinated N and P responses in Arabidopsis [38,39,[46][47][48]. Six SpGLK genes belonged to NIGT1 subfamily and were named based on their topological locations in the phylogenetic tree: , and SpHHO6 (SpGLK28). Subgroup IIe harbored 23 KANADI (KAN) genes which play key roles in organ positioning, cell type patterning and organ morphogenesis of shoot apical meristem (SAM) [49,50]. Subgroup IIf contained 35 members which functioned to leaf development [41], including 7 AtGLKs, 11 OsGLKs, 6 CeGLKs, 3 SpGLKs and 4 WaGLKs. All the Phosphate Starvation Response (PHR)/PHR1-like (PHL) TFs were clustered into group III, SpGLK3 (SpPHR1) and SpGLK12 (SpPHR2) exhibited closer relationships to AtPHR1 and OsPHR2, which were the center regulator of phosphate starvation response (PSR) in Arabidopsis and rice, respectively [13,51,52].

Chromosomal locations, duplication analysis, and synteny analysis of SpGARP genes
As shown in Fig. 2, a total of 35 SpGARP genes distributed in 15 chromosomes and one contig (tig00010334_1, SpGLK28), most chromosomes contain 1-3 GARP genes except chr10 (4 SpGARP genes, SpGLK15/16 and SpARR-B5/6) and chr11 (6 SpGARP genes, SpGLK17-21 and SpARR-B7). The duplication events of SpGARP genes were analyzed using McScanX [53]. Four segmental duplication events (SpGLK3-SpGLK12, SpGLK5-SpGLK14, SpGLK6-SpGLK13 and SpGLK17-SpGLK22) and two tandem duplication events (SpARR-B5-SpARR-B6, SpGLK17-SpGLK18) were identified, all the duplicated gene pairs belong to the same group. A total of 10 SpGARP paralogs gene pairs were identified based on the stringent amino acid sequence homology analysis (Table S2). We also identified the orthologous relationship of SpGARP genes with A. thaliana  (Table S2). The K a /K s ratios of the paralogs and ortholog gene pairs were analyzed to estimate the evolutionary pressure on SpGARPs. The K a /K s ratio ranged from 0. 15 3B and Table S2). The syntenic relationship of SpGARP orthologs in C. esculenta (Fig. 3C) and O. sativa (Fig. 3D) showed that group IId (green) and group III (brown) share the majority of orthologs gene pairs.

Conserved motifs, gene structure and phylogenetic analysis of SpGARPs in S. polyrhiza
An unrooted phylogenetic tree was constructed to visualize the evolutionary relationships between GARP members, using 35 SpGARP protein sequences (Fig. 4A). As shown in Fig. 4B, 15 conserved motifs were identified using MEME, motif 1 and 2 were  [45]. The B-motif only had one tryptophan (W) residue while the MYB domain had three and MYB-like domain had two W residues [2]. As shown in Fig. 4C, there was only one W residue in the B-motif and a consensus sequence ((A/K)SHLQ(K/M)) in the third helix in SpGARP TFs, as well as the B-motifs in Arabidopsis, rice and tobacco [2,41,42].
The structure of exons/introns were determined by aligning the genomic DNA sequences and the fulllength cDNA of SpGARP genes (Fig. 5A). There were 1 (SpGLK2/8) to 11 (SpARR-B7) exons in SpGARP genes, and most of the SpGARP genes contained 4-8 exons. The members belonging to the same clade always share the similar gene structure and conserved motifs, such as the members of clade I have 5-6 exons, SpGLK2 and SpGLK8 (group IIc) only have one exon, and the SpNIGT1s contain 4-5 exons.

cis-regulatory element analysis of the SpGARP gene family
GARP TFs perform diverse functions in plants, such as nutrient sensing, chloroplast development and circadian clock oscillation. The cis-regulatory elements of the promoter regions play important roles in the expression of GARP genes for the environmental signals, such as nutrient stress and abiotic stress [40][41][42]. As shown in Table  S3, 246 potential cis-acting elements were identified in the promoter regions of SpGARP genes. These cis-elements were involved in stress and hormone responses. As shown in Fig. 5B, most SpGARP genes contained two or more types of P/N starvation related cis-elements in their promoter regions. WRKY71OS was the most abundant P starvation related cis-elements of SpGARPs, as well as NODCON2GM in N starvation related cis-elements.
We also predicted the TFs that may be involved in regulating the expression of SpGARP genes. As shown in Fig. 6, most SpGARP genes interacted with multiple TFs, suggesting that they may be involved in many physiological processes. Among these predicted TFs, Apetala2  (Table S3 and Fig. 6). Several SpGARP genes were regulated by GARP superfamily TFs, six members (SpARR-B1/3 and SpGLK1/5/6/20) contained the ARR-B TF binding site

Transcriptional patterns of GARP genes in S. polyrhiza
The RNA-seq data of S. polyrhiza was analyzed to explore the expression patterns of SpGARP genes. As shown in  To further understand the function of SpGARP genes under different N and P conditions. The correlations between the expression patterns of SpGARP genes were analyzed (Fig. 8). A total of 69 SpGARP gene pairs showed correlation expression patterns, including 54 positive correlation and 15 negative correlation SpGARP gene pairs. The expressions of three IId members (SpGLK9, SpGLK25, and SpGLK27) were positive correlated, and both negative correlated with SpGLK6 and SpGLK21. Some SpGARP genes belonging to different groups presented correlated expression patterns, such as SpGLK7/10/14/17/18/20, SpARR-B1 and SpGLK2/6/11/21/28, SpARR-B3 and SpGLK14/23/24.

Expression analysis of PHR and NIGT1 genes in response to P or/and N deprivation by qRT-PCR
To further clarify the potential abilities of GARP genes responding to P and N stresses, the expression profiles of six NIGT1 subfamily genes (SpGLK6, SpGLK9, SpGLK13, SpGLK25, SpGLK27, SpGLK28) and two PHR genes (SpGLK3, SpGLK12) were verified using qRT-PCR. SpGLK13 has not been detected in the cDNA samples from both CG and nutrient stress treatments. As showed in Fig. 10, SpGLK3 was upregulated under N starvation

Discussion
GARP is a plant specific TF superfamily with diverse functions which have been identified in several species, including Arabidopsis, rice, maize, tobacco, cotton. However, it had not been reported in S. polyrhiza, which limited our understanding of the molecular mechanism of this plant to N/P responses. In curent research, we identified 35 SpGARPs belonging to three groups which were further classified into eight subgroups. The numbers of SpGARP (35) and WaGARP genes (28) were less than those of A. thaliana (56), C. esculenta (46), O. sativa (56), which was consistent with previous report of WRKY TF family in aquatic plants. The reason may be related to the minimal gene set of duckweed genomes (18,708 genes in S. polyrhiza, and about 15,000 genes in W. australiana) [25,54,55]. However, S. polyrhiza had more IId members than C. esculenta, O. sativa and W. australiana, suggesting that NIGT1 subfamily underwent expansion in the adapting to various aquatic environments. It should also be responsible for the high capability of nutrient uptake in S. polyrhiza. The SpGARPs that gathered in the same subfamily had the similar gene and protein structures, suggesting the strong evolutionary conservation of SpGARP genes [56]. ARR-B family was classified into group I and IIa. However, IIa and IIb showed closer phylogenetic relationship and the similar gene and protein structures, indicating the diverse origins of ARR-B proteins.
GARP members played important roles in many physiological processes as transcriptional regulators [2]. The earlier research found that GARP TFs are involved in Fig. 6 Regulation networks between SpGARP and potential transcription factors chloroplast development [2]. In recent years, their function on nutrient sensing and uptake had gained more attention in plant science [39,57]. PHRs were the center regulator in the response of Pi starvation, however, the transcriptional level expression of AtPHR1, OsPHR2 had not been improved under the Pi starvation condition [13,34,52]. SPX protein withhold PHRs in the cytoplasm under Pi sufficient condition to avoid the toxicity of high Pi, and release PHR proteins into nucleus so as to activate PSR under Pi deficient condition [58,59]. NIGT1 could  [38,46]. In the other hand, NIGT1 TFS acted as transcriptional repressors of N starvation response-related genes [47].
TFs regulated the expression of target genes at the transcriptional level through binding to the cis-element in the promoter region. Many types of cis-acting elements were present in the promoters of SpGARP genes, including abundant of P and N starvation related cis-elements, indicating that most SpGARPs are involved in hormone and stress responses, especially P and N starvation. The expression patterns of SpGARP genes were analyzed in this study, 18 out of 35 SpGARP genes were differential expressed between fronds and roots. Furthermore, most members of group IId and III were upregulated in roots. These results indicated that the root might play important roles in mineral sensing and uptake in S. polyrhiza. The N or/and P starvation altered the expressions of group I and IId genes, especially SpGARP6/13/9/25/27, which was consistent with previous researches in Arabidopsis [36,37]. Salinity stress is one of the major abiotic stresses limiting plant growth and productivity, it also suppresses the growth of S. polyrhiza [60]. There  [11], ZeGLK3 in Z. mays [43], 27 SiGLK genes in S. lycopersicum [61], and GhGLK55/120 in G. hirsutum [40]. In our study, we found that the expression of SpGLK25 was downregulated under salt stress, which helps to explain why the salt stress affects N/P uptake, and growth rate in plants [60]. The co-expression network showed the regulatory relationship between SpGARP and N, P response genes involving N and P sensing, uptake, and assimilation. From the results of qRT-PCR, most of the SpPHR and SpNIGT1 genes were differentially expressed under PS, NS, LP, and LN treatment. Interestingly, NIGT1s showed different expression patterns, for example, SpGLK6/28 were strong induced by NS while SpGLK9/25/27 were downregulated under NS. The different expression patterns between PS and LP, NS and LN suggested that these genes are involved in N/P sensing. N and P are crucial for plant growth and food production. Deficiency of N and P mineral nutrients is a huge problem for agriculture. Lots of botanists and agronomists focus their study and research on the adaptation mechanism of plants under the restriction of N or P nutrition. Nazir et al. (2016) found 25 N-deficiency induced proteins between the low-N tolerant and low-N sensitive maize genotypes [62]. Kunar et al. (2018) identified 37 N and P nutrition candidate genes in wheat, including 24 N-use efficiency (NUE) and 13 genes [63]. A total of 12 P-use efficiency (PUE) traits and 136 single nucleotide polymorphisms (SNPs) were identified among 144 diverse mungbean (Vigna radiata.) genotypes [64]. Meena et al. (2021) found that the relative expression of some P stress induced (PSI) genes in mungbean accession IC333090 (P-deficiency and drought stress tolerant accession) were significantly higher than that of sensitive accession IC488526, such as SPX1, sulfolipid sulfoquinovosyl diacylglycerol 1 (SQD1), Phosphate1 (PHO1), and purple acid phosphatase 1 (PAP1) [65]. Furthermore, PHT1, PHO1 and SPX genes were found to be response to P or/ and N deficiency in the aquatic crop S. polyrhiza [66,67].
S. polyrhiza is distributed throughout the world, and in the spotlight of plant science, environment remediation, biomass energy, and food security [17]. The discharge of agricultural, industrial and domestic wastewater led to the accumulation of N and P in water, resulting in dramatic water eutrophication and algae blooms. It also destroys the ecological balance of water bodies, lowers the content of dissolved oxygen, deteriorates water quality, and even causes the death of aquatic creatures. Phytoremediation has been recommended as an alternative solution to treat eutrophic wastewater due to its costeffective, environment friendly and sustainable characteristics. S. polyrhiza is widely used in phytoremediation for several outstanding characteristics: fast growth rate, adaptability to a wide range of environmental conditions, high nutrient uptake capacity, and higher proportion of carbohydrate [68]. S. polyrhiza can also adapt itself to various nutrient environments by balancing N and P uptake and assimilation. Our comprehensive analysis of GARP gene superfamily in S. polyrhiza showed that the members of SpNIGT1 subfamily would be the hub regulators of N and P sensing and acquisition.

Conclusion
In this study, we identified 35 GARP genes in S. polyrhiza genome, including 7 ARR-B subfamily genes and 28 GLK subfamily genes. The gene structures and phylogenetic analysis suggest a complex evolution history of this gene family in S. polyrhiza. As the central regulator in P and N nutrients, the PHR/PHL subfamily genes were transcriptional induced by NS and LN treatments, while the NIGT1 subfamily genes could respond to P and N stresses, especially SpGLK9/25/27 which were upregulated under PS treatment and downregulated under NS treatment. This study provides insight into the evolution and function of GARP superfamily in S. polyrhiza, and facilitates the further functional verification of SpGARP genes.

Gene duplication and K a /K s analysis
The information regarding chromosome length and gene locations of GARP family genes in S. polyrhiza was extracted from the Generic Feature Format (GFF) files. The duplication events were defined based on the collinearity analysis of candidate gene pairs using MCS-canX (http:// chibba. pgml. uga. edu/ mcsca n2/) [53]. The synteny analysis of SpGARP genes with CeGARPs and OsGARPs was performed using MCScanX and visualized using Circos (http:// mkweb. bcgsc. ca/ circos/ table viewer/) [76]. The non-synonymous substitution rate (K a ) and synonymous substitution rate (K s ) of the duplication and orthologous gene pairs were calculated using PAMLX (http:// abacus. gene. ucl. ac. uk/ softw are/ paml. html) [77].

Gene structure analysis and identification of conserved motifs
Multiple Expectation Maximizations for Motif Elicitation (MEME, http:// meme-suite. org/) was employed for the analysis of conserved motifs in GARP proteins with the following parameters: maximum number of motifs, 15; motif length, 6 to 50 amino acids [78]. The structure of GARP genes, including intron and exon information, was visualized using the online tools Gene Structure Display Server 2.0 (http:// gsds. gao-lab. org/ index. php/) [79].

Promoter analysis of SpGARP genes
Cis-acting elements in the promoter regions of GARP genes (2000 bp upstream of the start codon) were predicted and analyzed in New PLACE (https:// www. dna. affrc. go. jp/ PLACE/? action= newpl ace/) [80]. The subset of data representing P and N starvation related to cis elements was visualized using TBtools [81]. To discover the TFs involving in the regulatory expression of SpGARP genes, the online tool PlantRegMAP (http:// plant regmap. gao-lab. org/ bindi ng_ site_ predi ction. php) was used to predict the potential binding sites of TFs in the promoter regions of SpGARP genes [82]. Then, the regulatory networks between SpGARP genes and potential TFs were presented using Cytoscape 3.7.0 [83].

RNA-seq atlas analysis
The temporal and spatial expression profiles of SpGARPs in different tissues/organs (leaves, roots, and stipule, Bio-Project PRJNA557001), under various nutrient stresses (supply with P and N (+ P/ + N), N deprivation (+ P/ − N), P deprivation (− P/ + N), P and N deprivation (− P/ − N), and none nutrient supply (H 2 O) for seven days, PRJNA724886), under salt stress (salt-0 h, salt-6 h, salt-12 h, ssalt-24 h, PRJNA563960) were obtained using the publicly available transcriptome data from NCBI [25,60]. The heat map showing the correlations between the expression patterns of SpGARP genes were generated with the ggplot package in R (version 4.1.2, https:// www.R-proje ct. org/). The significant correlated gene pairs should satisfy the following criteria: the absolute of the correlation value > 0.85, p-value < 0.05. The genes whose coverage TPM > 1 were filtered and used for WGCNA. Co-expression network modules were identified with the WGCNA package 1.63 in R to generate the signed weighted correlation network, and the network of genes was visualized in Cytoscape 3.7.0 [83,84].

Plant materials and treatments
S. polyrhiza strain 7498, which was gifted from Duckweed Stock Cooperative (http:// www. ruduc kweed. org/ datab ase. html) and stored in National Aquatic Biological Resource Center (http:// www. nabrc. ihb. ac. cn/), was used as the source of plant materials in the study. S. polyrhiza was cultivated in liquid half-strength MS solution at pH 5.8, under the conditions of 16 h/8 h photoperiod (day/night), irradiance of 85 µmol photons•m −2 •s −1 , and temperature of 25 °C. Ten days later, duckweed was treated in half-strength MS solution without P (PS treatment), half-strength MS solution without N (NS treatment), half-strength MS solution with 1 μM KH 2 PO 4 (LP treatment), half-strength MS solution with 1 μM NH 4 NO 3 (LN treatment). The samples were harvested at varied time points (0 (Control group, CG), 1, 3, 5, 7 and 10 days) after the treatment, and immediately frozen in liquid nitrogen and stored at − 80 °C for further analyses. Three samples were collected for each treatment at each time point.

RNA isolation and qRT-PCR analysis
The total RNA of duckweed was extracted using Ominiplant RNA Kit (CoWin Biosciences, Beijing, China). From total RNA, first-strand cDNA was synthesized using a PrimeScript ™ RT reagent Kit (TaKaRa, Dalian, China). The primers of SpGARP genes were presented in Table S4. qRT-PCR program was performed using the same methods in previous study [85]. Each reaction was analyzed in triplicate and the 2 −△△CT method was used to analyze the data [86]. The qRT-PCR results were statistically analyzed using SPSS 25.0 software. Significance differences were determined by one-way ANOVA and a Fisher's LSD test at the p < 0.05.