The nine citrus rootstock varieties, including Poncirus trifoliata (Pt026, Pt030, Pt034, Pt038), Citrus wilsonii Tanaka ‘Zhique’, C. limonia Pasquale ‘Volkamer’, C. limonia Osbeck ‘Canton Lemon’, C. reshni Hort. Ex Tan ‘Cleopetra’, and C. reticulata ‘Zhuju’, were used in this study. The plant materials were collected from National Citrus Germplasm Repository (Chongqing, China).

The seedlings of ‘Volkamer’ were used for analysis of genes expression in different tissues and organs, and under different hormone treatments and stresses. The seedlings of all varieties were used to study the correlation between root indexes and gene expression.

The genome and protein information of citrus plant (C. clementina) was downloaded from the Phytozome database (http://phytozome.jgi.doe.gov/pz/portal.html). Twenty-four response regulators’ protein sequences of Arabidopsis thaliana (http://www.arabidopsis.org/) (ARRs) were used as reference sequences to identify the members of the RR family in C. clementina (CcRRs) with protein basic local alignment search tool (BLASTP) (Mahram and Herbordt, 2015). The top E-value was less than 1 × 10⁻¹⁰. All the retrieved sequences were submitted to the HMMER database (https://www.ebi.ac.uk/Tools/hmmer/; Potter et al., 2018), Pfam database (http://pfam.xfam.org/; Punta et al., 2011) and NCBI CDD database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi; Marchler-Bauer et al., 2015) to confirm whether they were the members of RR family. The redundant members were removed manually. ExPASy website (http://web.expasy.org/) was used to predict molecular weight (MW) and isoelectric points (pI) of all proteins of the identified RR family members. CELLO website (http://cello.life.nctu.edu.tw/) was used to predict the subcellular localization of RR family members (Yu et al., 2006).

The chromosomal location of each RR gene was derived from the genome annotation information. Chromosomal distributions of these genes were generated using MapChart software (Voorrips, 2002).

MEGA7.0 software (Kumar et al., 2016) was used to construct the phylogenetic tree with protein sequences of RRs from different plant species based on the Maximum Likelihood Method (MJ). The bootstrap value was repeated 1000 times.

TBtools (Chen et al., 2020) was used to display the distribution of gene structure of RRs. The analysis of the protein conserved motif was conducted through the online MEME website (https://meme-suite.org/meme/tools/meme; Bailey and Elkan, 1995) to study the differences among RR family members. The maximum number of motifs was set to 5. 2000 bp upstream sequence of each RR gene was extracted from the citrus genome. The cis-elements in the promoter region were predicted using PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

After removal of seed coats, the embryos were cultured in petri dish at 28°C under moisturized conditions for 3-5 days for germination. The seedlings of Volkamer were grown in hydroponics for 25 days. The roots, stems and leaves were collected for gene expression analysis of CcRRs. Meanwhile, the remaining seedlings were treated with hormones and abiotic stresses. The root samples were collected at different time (0, 1, 3, 6, 12 h) after 50 μM 6-BA, 100 μM ABA, and 20% polyethylene glycol (PEG₆₀₀₀) treatments for RNA extraction, respectively. The roots were also sampled from the seedlings treated at 4°C and 200 mM NaCl for 0, 3, 6, 12 and 24 h, respectively. In addition, the seedlings with lateral root primordium just emerged (as shown in Supplementary Figure 1 ) were used to explore the expression of type A and type B RRs in different zones of root (RT: the meristematic/elongation zone; RM: the root elongation/differentiation and lateral root initiation zone; RC: lateral root growth zone). This division was based on the reports of Hwang et al. (2012), Goh (2019), Santos Teixeira and ten Tusscher (2019).

The seeds of nine citrus rootstocks under same state were selected and sown in the pots (30 × 30 cm) filled with the mix of vermiculite and perlite. Different zones of root were sampled at 25, 50 and 75 days after sowing. The samples were immediately frozen with liquid nitrogen and stored at -80°C for RNA extraction to explore the expression of the three RR genes. At the same time, the length of primary root and the number of lateral roots in RC zone of the root were counted. All the seedlings were grown at 28°C in a greenhouse under 16 h day length.

Total RNA was isolated from samples using an EASYspin Plus Plant RNA Rapid Extraction Kit RN38 (Aidlab, Beijing, China) following the manufacturer’s instructions. Subsequently, the RNA was reverse transcribed using the Thermo Scientific RevertAid MM (Thermo Fisher, Lithuania). qRT-PCR was performed with 10 μL reaction mixture containing 2 μL cDNA sample, 5 μL universal SYBR^® Green Supermix (2 × iTaq™, Bio-Rad, CA, USA), 0.2 μL of each primer (10 μM) and 2.6 μL sterilized water. The qRT-PCR program was as follows: 95°C for 60 s, followed by 39 cycles of 95°C for 20 s, 60°C for 30s. The primers used in this study are listed in Supplementary Table 1 . Each treatment included three replications. The expression levels of the genes were normalized with citrus Actin gene, and fold changes were calculated using 2^-ΔΔCt method (Zhang et al., 2020).

In this study, more than 100 RRs and RR like proteins were searched out from the C. clementina genome. The RRs with sequence redundancy or alternative splicing forms were removed. Furthermore, the RRs containing other domains, or lacking conserved motifs were also eliminated. Finally, a total of 29 potential RR proteins (including 8 pseudo RRs) were identified, and named as CcRR1-21 and CcPRR1-8 according to their positions on the scaffolds. The details of CcRRs, including their names, accession numbers, number of amino acids, molecular weight, isoelectric point (pI), grand average of hydropathicity (GRAVY) and putative subcellular localization were listed in Table 1 . The number of amino acids of CcRR proteins varied from 140 (CcRR12) to 873 (CcPRR7), with molecular weight (MW) from 15.29 (CcRR12) to 95.44 kDa (CcPRR7) and pI from 4.88 (CcRR19) to 8.92 (CcRR21). The prediction of subcellular localization showed that most of CcRRs were located in the nucleus, and a few of them were located in the plasma membrane, mitochondria or chloroplast.

To uncover the phylogenetic relationship of CcRRs, a phylogenetic tree was generated with 95 RR proteins, of which 33 were from Arabidopsis (24 ARR and 9 APRR), 33 from poplars (22 PtRR and 11 PtPRR) (Ramírez-Carvajal et al., 2008), and 29 from citrus (21 CcRR and 8 CcPRR) ( Figure 1 ). It showed that the 95 RRs were divided into 9 clusters. Based on the feature of protein structure, the 9 clusters could be classified into 4 types. The group in green represented type A CcRRs with only one D-D-K domain, including CcRR5-6, CcRR8, CcRR11-13 and CcRR19. Type C CcRR (pink) was mostly close to type A with a similar D-D-K domain. It is a small group with only four members in citrus (CcRR3, CcRR7 and CcRR17-18). CcRR1-2, CcRR4, CcRR9-10, CcRR14-16 and CcRR20-21 were placed in type B CcRRs (blue) with D-D-K domain and MYB-like DNA-binding domain. The group marked in orange color contained pseudo RRs with D-D-K domain and CCT domain. Chromosome mapping results showed that most of CcRR members were clustered on scaffold 3, 4 and 9, and a few of them were on scaffold 1, 2, 5, 7 and 8 ( Figure 2 ).

The diversity of exon-intron structure and motifs may reveal the evolution of gene family (Wang et al., 2017; Sun et al., 2022). Intron-exon analysis of CcRR genes ( Figure 3 ) indicated that the CcRRs classified into type A and B possessed 3 to 7 intronic regions, while there was only one intron in the type C CcRRs. The numbers of introns for CcPRRs varied more, ranging from 5 to 10. Five potential conserved motifs were identified from 29 CcRRs by online MEME program. The details of each motif were shown in Supplementary Figure 2 . As shown in Figure 3 , motifs 1, 2 and 4 encoding the D-D-K domain presented in most of CcRRs, but the sequences of D-D-K domain in some of CcRRs were incomplete. CcRR12 of type A, CcRR9 and CcRR21 of type B, CcPRR5 and CcPRR6 of pseudo RRs, and all type C CcRRs were missing part of the D-D-K domain sequence. All the type B CcRRs and CcPRR8 also possessed motif 3, which was identified as a conserved MYB-binding domain. Motif 5 was deemed as a conserved CCT motif and it existed in almost all CcPRRs except CcPRR4 and CcPRR6. These results suggested that the four types of CcRRs may play different roles in citrus.

Cis-element in promoter region of the gene plays an important role in the regulation of gene transcription and usually combines with TFs to determine the level of gene expression (Wang et al., 2017; Dai et al., 2022). Results of cis-element prediction showed that the promoters of CcRR family members contain a large number of light responsiveness elements ( Figure 4 ). Drought-induced, low-temperature response, defense and stress response elements were also found in promoter regions of CcRRs. In addition, there are many hormone response elements, such as abscisic acid, auxin, methyl jasmonate, salicylic acid, gibberellin response elements and zein metabolism regulation element. Some specific cis-acting elements were also identified, such as anaerobic induction, meristem development, and endosperm expression elements. The analysis of the predicted cis-elements in the promoter indicated that the transcriptional regulation of CcRR gene family members was probably related to the processes of growth and development of citrus plants, and the responses to various stresses.

Levels of gene expression in the roots, stems and leaves were analyzed by qRT-PCR to characterize the expression of type A and type B CcRR genes in different tissues of citrus ( Figure 5 ). CcRR12, CcRR13, CcRR20 and CcRR21 expressed at very low level in roots, suggesting that they are not active in this tissue, therefore, they were excluded from further work of this study. The remaining 13 CcRR genes showed various expression patterns in different tissues, and most of them had the highest expression in leaves, such as CcRR4-6, CcRR10, CcRR11, CcRR15 and CcRR16. Among them, CcRR4, CcRR5, CcRR6 and CcRR16 also highly expressed in roots.

Phytohormone and abiotic stress treatments could induce the expression of CcRR genes. Analysis of RR genes’ expression under hormone treatments and abiotic stresses showed that the transcriptional levels of CcRR4, CcRR5, CcRR6, CcRR16 were high and changed significantly under the conditions of PEG, NaCl, cold, ABA or 6-BA treatment, indicating that those CcRRs dynamically responded to the multiple stresses ( Figure 6 ).

Under 6-BA treatment ( Figure 6A ), the expressions of all the type A CcRRs (CcRR5, CcRR6, CcRR8, CcRR11, CcRR19) and some of type B CcRRs (CcRR2, CcRR15, CcRR16) were up-regulated at the early stage of the treatment and then declined, while the expressions of CcRR1, CcRR4, CcRR9, CcRR10 and CcRR14 were generally suppressed by the treatment. For ABA treatment ( Figure 6B ), the expression patterns of seven CcRRs (CcRR2, CcRR4, CcRR5, CcRR11, CcRR14, CcRR15, CcRR16) were basically the same. Their expressions up-regulated at the beginning of treatment and then declined, while the expression of CcRR19 was increasing throughout the treatment. However, CcRR1 and CcRR10 showed the opposite trend of decreased expression, while the expressions of CcRR6, CcRR8 and CcRR9 were only down-regulated at the early stage, but increased afterward.

Under PEG treatment ( Figure 6C ), the expressions of four genes (CcRR2, CcRR5, CcRR8 and CcRR19) were all up-regulated at 1h or 3 h after treatment, but down-regulated from 3h or 6h after treatment, while the expressions of CcRR4 and CcRR9 were up-regulated nearly in the whole period of treatment. On the opposite, the expression of CcRR10 was gradually decreased during the treatment, while the expressions of five genes (CcRR1, CcRR6, CcRR11, CcRR15, CcRR16) were down-regulated at the early stage of treatment and then recovered at late stage. Under low temperature condition ( Figure 6D ), the expressions of CcRR2, CcRR14 and CcRR16 increased in the entire period of treatment, while the expression level of type A CcRRs (CcRR5, CcRR6, CcRR8, CcRR11, CcRR19) and CcRR4 showed a trend of increasing first and then decreasing. However, the expressions of CcRR1 and CcRR10 were consistently decreasing, while the expressions of CcRR9 and CcRR15 were down-regulated at the early stage of treatment and then slight increased. For NaCl treatment ( Figure 6E ), the similar expression patterns were observed in the most of genes (CcRR2, CcRR4, CcRR5, CcRR6, CcRR8, CcRR9, CcRR14, CcRR15, CcRR16 and CcRR19). Their expressions were also up-regulated at the beginning of treatment and then declined at varying degree. The expressions of CcRR1 and CcRR11 declined during the whole period of treatment, while the expression of CcRR10 was down-regulated at the early stage of treatment and then up-expressed. The results of above observation indicated that most of CcRR genes involved in the processes responding to the various abiotic stresses and exogenous hormone treatments.

The eight CcRRs highly expressed in roots, also drastically respond to the abiotic stresses and exogenous phytohormone treatments in above study were selected to explore their functions involved in root development with Volkamer seedlings. The levels of gene expression were investigated in RT, RM, and RC zones of the roots as shown in Figure 7 . CcRR2, CcRR5, CcRR6 and CcRR19 predominately expressed in RM, among which the expression of CcRR5 was the most prominent. The expressions of CcRR4, CcRR16, especially CcRR10, were higher in RC than that in RT. However, expression of CcRR14 was significantly higher in both RT and RC than in RM. The results suggested that those CcRRs probably play different roles in root morphogenesis. Therefore, the three CcRRs (CcRR5, CcRR10 and CcRR14) that distinguishingly expressed in different zones of root were further investigated to determine their roles in root morphogenesis of 9 citrus rootstocks, which varied in the length of primary roots and the numbers of lateral roots. Results of analysis indicated that the expression of CcRR5 in RT zone was significantly correlated to the length of primary root, while the expression of CcRR10 in RC zone was significantly correlated to the number of lateral roots, and the expressions of CcRR14 were significantly correlated to the two root indexes ( Figure 8 ; Supplementary Table 2 ).