Pangenome-wide analysis of cyclic nucleotide-gated channel (CNGC) gene family in citrusSpp. Revealed their intraspecies diversity and potential roles in abiotic stress tolerance

1 Introduction

Calcium is an important macronutrient for plant growth and development and is involved in signaling pathways as a secondary messenger. It also plays a key role in the defense mechanism of plants against abiotic stress (Lecourieux et al., 2006;Kudla et al., 2018）.Calcium sensor proteins belong to three main families including calmodulin (CaM) and calmodulin-like proteins (CMLs) (Yang and Poovaiah, 2003;Bender and Snedden, 2013), calcineurin-B-like proteins (CBLs) (Luan, 2009), calcium dependent protein kinases (CPKs) and calcium and calmodulin dependent protein kinase (CCaMK) (Cheng et al., 2002;Wang et al., 2015）.Calcium binding to these calcium sensors induces a conformational change that triggers either a particular target protein or directly stimulates kinase activity by taking into account CPKs (Ranty et al., 2016）.In contrast, several families of ion channels regulate the uptake of calcium including Cyclic nucleotide-gated channels (CNGCs), two pore channel 1 (TCP1), ionotropic glutamate receptors, and several other channels (Demidchik et al., 2018）.

CNGCs belong to the nonselective cation channels that are found in both animals and plants. Plant CNGCs was first discovered in 1998 while scanning calmodulin-conjugated transporters (HvCBT1) in barley (Mäser et al., 2001）.CNGCs are ligand-gated channels that are calcium permeable and involved in the interaction of cyclic nucleotides and calcium dependent signaling pathways (Talke et al., 2003）.CNGCs are calcium sensors in eukaryotes while calcium is important for plant growth, development, light signaling, drought and salt stress, and pathogen tolerance (Ranty et al., 2016）.CNGCs get activated by the binding of cyclic nucleotides (cNMP) and their activity gets inhibited by Ca²⁺/CaM binding (Trudeau and Zagotta, 2002）.Calcium is very helpful in regulating plant growth under stress conditions. There are 6 TM domains (S1-S6) and a pore region in CNGCs, fifth and sixth domains along with the Cyclic nucleotide-binding domain (CNBD) and CaM binding domains are present at C-terminal. CNBD comprises a phosphate binding cassette (PBC) and a hinge region (Duszyn et al., 2019）.中国人民银行结合糖和磷酸根of cyclic nucleotide binding (CNB) ligand and the hinge region contributes to the efficacy of ligand binding and selectivity (Li et al., 2019）.CNGCs are also involved in plants responses to various abiotic and biotic stress conditions. (Jha et al., 2016）.

CNGC gene family has been reported inArabidopsis thaliana(Mäser et al., 2001),Brassica oleracea(Kakar et al., 2017),Zea mays(Hao and Qiao, 2018),Ziziphus jujube Mill.(Wang et al., 2020),Nicotiana tobacum L.(Nawaz et al., 2018),Triticum aestivum L.(Guo et al., 2018), Oryza sativa(Nawaz et al., 2014), Brassica rapa(Li et al., 2019), Pyrus bretschneideri Rehd(Chen et al., 2015）.andSolanum lycopersicum(Saand et al., 2015）.On the basis of the phylogenetic classification in the aforementioned plants, this gene family is classified into four major groups and the fourth group is further divided into two sub-groups namely as; I, II, III, IV-A, IV-B. A single reference genome is not enough to capture diversity present among the members of a species (Golicz et al., 2016）.Thus, it brings a bias to study gene family members in plants solely based on a single genome. Therefore, it is suggested to conduct pangenome-wide analysis for gene family characterization (Tahir Ul Qamar et al., 2019）.The first ever concept regarding pangenome was introduced when the pangenome ofStreptococcus agalacitaewas developed (Tettelin et al., 2005）.Pangenome of a species comprises core genes that are present in all members, accessory genes that are present in few but not in all members, and unique genes that are present only in specific members (Tahir ul Qamar et al., 2020;Ismail et al., 2022;Zanini et al., 2022）.

Citrus is an economically important fruit crop as it is widely used both as a fruit and as a juice (Liu et al., 2019）.它是多年生作物和主要种植在China, Brazil, India, United States, Mexico, Spain, and Italy (Liu et al., 2012）.Citrinae是一大群柑橘类的水果树that belong to the subfamily Aurantioideae and the family Rutaceae. Based on botanical features Citrinae is categorized into three types i.e., primitive citrus, near citrus, and true citrus (Wang et al., 2017）.The well-known Citrus varities include;Atlantia buxfolia(Chinese box orange),Citrus sinensis(sweet orange),Citrus grandis(pummelo),Citrus recticulata(mandarin),柠檬(lemon),Citrus paradisi(grapefruit) andPoncirus trifoliata(Trifoliate orange) (Liu et al., 2019）.Citrus varities widely influenced by drought stress as the productivity, growth, and yield of citrus get reduced after facing drought stress (Osakabe et al., 2014）.However, few drought resistant varities are also reported which can withstand against this stress, including navel orange and trifoliate orange (Bhusal et al., 2002;Koshita and Takahara, 2004;Pingping et al., 2017）.

In present study,C. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliatawere selected for pangenome-wide analysis of CNGCs gene family, as they have good quality assembled genomes and their annotations are available at chromosome level. The quality of genome assembly or sequencing directly affects the quality of results (Vaattovaara et al., 2019), therefore, the aforementioned species were preferred to reduce the biasness. CNGCs gene family has been studied in several plant species at single genome-wide level (Mäser et al., 2001;Nawaz et al., 2014,2018;Chen et al., 2015;Saand et al., 2015;Kakar et al., 2017;Guo et al., 2018;Hao and Qiao, 2018;Li et al., 2019;Wang et al., 2020), but no pan-genome-wide analysis has been performed before. Therefore, current study aims to provide a comprehensive pangenome-wide representation of CNGCs gene family in citrus species, which will serve as the foundation for future gene family researches.

2 Materials and methods

2.1 Identification of cyclic nucleotide-gated channel family genes inC. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliata

20 CNGC protein sequences ofA. thalianataken from TAIR database (https://www.arabidopsis.org/) (Rhee et al., 2003) were used as query and BLASTp search was performed on Citrus pan-genome to breeding database (CPBD;https://citrus.hzau.edu.cn/) (Liu et al., 2022) againstC. sinensisv2.0,A. buxfoliav2.0,P. trifoliatav1.0,C. recticulatav2.0, andCitrus grandis (L.) Osbeck.cv.Wanbaiyouv1.0. The resulting BLAST hits were manually processed to remove duplicates and isoforms and the final hits were used for further analyses.

To check the presence of specific domains, databases including SMART (https://smart.embl-heidelberg.de/) (Schultz et al., 2000), CDD (https://pfam.xfam.org/) (Marchler-bauer et al., 2011), and HMMER (https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan) (Potter et al., 2018) were used. This eliminated those sequences that didn’t have specific conserved domains required for CNGC protein function. Domain architecture was constructed using the HMMER database. Molecular weight (MW), Theoretical isoelectric point (PI), Instability index (II), Aliphatic index (AI), and Grand average of hydropathy (GRAVY) were determined by using the web-based tool ProtParam available at the EXPASY server (https://web.expasy.org/protparam) (Gasteiger et al., 2003）.Subcellular localization was determined using CELLO version 2.5 (https://cello.life.nctu.edu.tw/) (Yu et al., 2006）.

2.2 Multiple sequence alignment and phylogenetic analysis

To comprehend the phylogenetic relationships of identifiedCNGCs, multiple sequence alignment of identified CNGC protein sequences ofC. sinensis,A. buxfolia,C. recticulata,C. grandis,P. trifoliataalong with already reported protein sequences ofO. sativa(Nawaz et al., 2014),Z. jujuba(Wang et al., 2020), Z. mays(Hao and Qiao, 2018),A. thaliana(Köhler and Neuhaus, 2000) andP. bretschneideri(Chen et al., 2015) was done using ClustalW program and a phylogenetic tree was constructed by using online server IQ-tree (https://iqtree.cibiv.univie.ac.at/) (Nguyen et al., 2015) with Maximum Likelihood (ML) method and 1,000 replicates while other parameters were set to their default values. The tree was visualized and edited using the online server iTOL (https://itol.embl.de/) (Letunic and Bork, 2021）.

2.3 Chromosomal location, gene structure, and conserved motif analysis

The chromosomal location, start and end sites ofC. sinensis,A. buxfolia,C. recticulata,C. grandis,andP. trifoliatawere retrieved from the CPBD database and a genetic linkage was constructed by using TBtools (Chen et al., 2020）.The gene and CDS sequences ofC. sinensis,A. buxfolia,C. recticulata,C. grandis,andP. trifoliatawere retrieved from the sequence fetch option at the CPBD database (https://citrus.hzau.edu.cn/) (Liu et al., 2022）.The GSDS v2.0 (https://gsds.gao-lab.org/) (Hu et al., 2015) was used for the visualization of gene structures of CsCNGCs, AbuCNGCs, CreCNGCs, CgCNGCs, and PtCNGCs. Conserved motifs were identified by using MEME (Multiple EM for Motif Elicitation) suite 5.4.1 (https://meme-suite.org/meme/db/motifs) (Bailey et al., 2009）.All parameters were set to their default values except the number of motifs that were set to 10.

2.4 Gene duplication and promoter analysis

The location ofCNGCgenes inC. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliatawas retrieved from the CPBD database (https://citrus.hzau.edu.cn/) (Liu et al., 2022）.All genes possessing ≥70% sequence identity were considered duplicated genes (Hu et al., 2021）.DnaSP v6.0 (Librado and Rozas, 2009) offline tool was used to calculate the rate of Non-synonymous (Ka) and synonymous substitutions (Ks) of duplicated gene pairs. To calculate the selection pressure that assisted in the evolution of theCNGCgene family Ka/Ks ratio was used. The formula for calculating duplication time was the following:T=Ks/2x(wherexrepresents substitutions per synonymous site per year and is equal to 6.56 × 10⁻⁹for dicots) (He et al., 2016）.The cis-elements in 2000bp coding regions ofCsCNGCs,CreCNGCs,CgCNGCs,AbuCNGCs,andPtCNGCswere retrieved from the Citrus pan-genome to breeding database (CPBD,https://www.citrus.hzau.edu.cn/) (Liu et al., 2022）.While the types, numbers, and functions of these cis-elements were analyzed by using PlantCare web-based tool (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002）.

2.5 Putative miRNA target prediction, protein-protein interaction network, and gene ontology analysis of citrusSpp.

Plant microRNA Encyclopedia (PmiREN;https://pmiren.com) database was utilized to acquire mature miRNA sequences ofC. sinensis.For putative miRNA target prediction CDS sequences of the potential target,CsCNGCswere utilized and were submitted at the psRNATarget server (https://www.zhaolab.org/psRNATarget/home) (Dai et al., 2018) along with the respective mature miRNA sequences ofC. sinensiswith default considerations. The regulatory association between targetCsCNGCsand predicted miRNAs was visualized using Cytoscape software (Shannon et al., 1971）.The interaction among members of the CNGC protein family and other proteins from the citrus plant was predicted using the STRING database (https://string-db.org/）.32 CsCNGC protein sequences were uploaded to the STRING database with ‘Citrus sinensis’ being selected as reference species. The level of connection used was sixth and other parameters were kept by default. PPI network was visualized and edited using Cytoscape software (Shannon et al., 1971）.Citrus Pan-genome2breeding database (CPBD;http://citrus.hzau.edu.cn/) (Liu et al., 2022) was utilized to analyze gene ontology (GO) enrichment of CitrusSpp.using the gene IDs ofCNGC基因。

2.6 Expression profiling ofC. sinensisunder drought stress

To demonstrate the expression ofC. sinensisunder abiotic stress (drought) in leaves, RNA-seq data was downloaded from the NCBI-SRA database (https://www.ncbi.nlm.nih.gov/sra) (BioProject: PRJNA792482). Reference genome and GFF3 files were downloaded from the Citrus pan-genome to breeding database (CPBD,https://citrus.hzau.edu.cn/) (Liu et al., 2022）.To check the quality of paired-end data (in FASTQ format) FASTQC was utilized and Trimmomatic was used for trimming and improving the quality of reads. Then HISAT2 was used for the alignment of reads to the C.sinensis v2.0genome. To normalize gene expression in terms of Fragments per kilobase of transcripts per million mapped reads (FPKM) Cufflinks were used. The heatmap was constructed using pheatmap function of R-language (Ihaka and Gentleman, 1996）.

2.7 Drought stress treatment, ribonucleic acid isolation, and quantitative real-time reverse transcription–polymerase chain

柑橘类植物生长中的控制ental conditions in a growth chamber (having 60 ± 3% humidity, 27 ± 2°C temperature, and 5000 LUX light intensity) with recommended fertilizer and water treatment. Four months old citrus plants were subjected to drought stress and leaves were collected and 0, 10, and 20 days of drought stress. Control and drought-stressed leaves were harvested for RNA extraction. Zomanbio (Cat no. ZP401-2) total RNA-pure reagent (Lot#200F12F) was used to extract total RNA and the complementary DNA (cDNA) was synthesized by using Zomanbio (M-MLV, ZR102-3) reverse transcriptase kit (Beijing, ZOMAN Biotechnology Co., Ltd.) according to the manufacturer instructions. For quantitative real-time polymerase chain reaction (qRT-PCR) ChamQ universal master mix SYBR (Vazyme, Q711-02) and LongGene (Model: q2000b) fluorescence quantitative PCR instrument (Langji Scientific instrument Co., Ltd.; Hangzhou, China) were used whereas citrus actin gene was used as an internal reference. 2^-(ΔΔCt) method was applied to analyze the qRT-PCR expression data in Excel (Microsoft Corp., Redmond, WA, United States). Statistix 8.1 (Tallahassee Florida, United States) statistical software was used for analyzing all qRT-PCR data and the Excel program was used for graphs. The qPCR primer information is characterized (Supplementary Table S1）.

2.8 3D Structure prediction of cyclic nucleotide-gated channels in citrus spp.

Three-dimensional (3D) structures of 13 CNGC proteins were predicted, including 9 proteins fromC. sinensis, one fromA. buxfolia,and two fromP. trifoliata.Among these 13 proteins, 3D structures of 12 CNGC proteins were predicted by using Alphafold2 (https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb) (Jumper et al., 2021) Whereas, the 3D structure of CsCNGC1.4 was predicted by using trRosetta (https://yanglab.nankai.edu.cn/trRosetta/) (Du et al., 2021)由于其长度即1427 aa。蛋白质的方法es were visualized by using Pymol (Yuan et al., 2017）.For validation of these predicted structures SAVES server (https://saves.mbi.ucla.edu) was used.

3 Results

3.1 Identification of cyclic nucleotide-gated channel genes inC. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliata

A total of 32 putative genes inC. sinensis, 27 genes inC. recticulata, 30 genes inC. grandis, 31 genes inA. buxfolia,and 30 inP. trifoliatawere identified. The identifiedCNGCgenes were named based on their phylogenetic relationships withCNGCsinArabidopsis.Figure 1is showing the homologs ofArabidopsis CNGCgenes present in five species under study. Most of the identified members ofArabidopsis CNGCsare present in five species under study except forAtCNGC3,AtCNGC6,AtCNGC9,AtCNGC11,AtCNGC12,andAtCNGC20. All other members have a variable number of homologs present in five Citrus species. Further, two unique genes were identified:CNGC13andCNGC14.CNGC13is present in three plant species includingC. sinensis,A. buxfolia,andP. trifoliatawhile absent inC. grandisandC. reticulata.CNGC14is present in only one plant species,P. trifoliatawhile being absent in the other four species.

Conserved domains that were predicted inC. sinensis,A. buxfolia,C. recticulata,C. grandis,andP. trifoliatainclude Cyclic Nucleotide Binding Domain (CNBD or cNMP), Ion trans (IT), Cap family effector domain (CAP_ED) and other ion trans domains (Supplementary Table S2）.Ion trans and cNMP binding domains were the most conserved among all. Domain architecture was constructed according to the prediction results of the HMMER database. cNMP binding domain was not present in CsCNGC1.1, CsCNGC1.2, CsCNGC1.5, CsCNGC1.7, CsCNGC1.8, CsCNGC1.9, CsCNGC2.1, CsCNGC2.2, CsCNGC4.1, CsCNGC7, CsCNGC15.2 and CsCNGC19 according to prediction results of HMMER database but SMART database prediction confirms the presence of cNMP binding domain in these proteins. The domain architecture ofC. sinensisis given in (Figure 2）.

Ion trans domain was absent in CreCNGC1.2, and CreCNGC1.6 according to HMMER database prediction. While CDD prediction confirms the presence of the Ion trans domain in CreCNGC1.6. Results of the HMMER database demonstrate the absence of the cNMP binding domain in CreCNGC1.1, CreCNGC1.3, CreCNGC1.4, CreCNGC2.1, CreCNGC4.1, CreCNGC7, CreCNGC15.1, CreCNGC15.3, and CreCNGC19 according to prediction. Domains predicted by the SMART database indicated the presence of the cNMP binding domain in these proteins (Supplementary Figure S1）.

The following proteins ofC. grandisCgCNGC1.1, CgCNGC1.2, CgCNGC1.3, CgCNGC1.4, CgCNGC1.5, CgCNGC1.6, CgCNGC2.1, CgCNGC2.2, CgCNGC4.1, CgCNGC7, CgCNGC15.1, CgCNGC15.3, and CgCNGC19 didn’t have cNMP binding domain as per HMMER database prediction. Whereas, the cNMP binding domain was predicted to be present in all these proteins except CgCNGC1.3 as supported by SMART prediction (Supplementary Figure S2）.

AbuCNGC1.1、AbuCNGC1.3 AbuCNGC1.5 AbuCNGC1.7, AbuCNGC1.8, AbuCNGC1.9, AbuCNGC2.1, AbuCNGC4.1, AbuCNGC7, AbuCNGC15.1, AbuCNGC15.3, and AbuCNGC19 are those aforementioned proteins ofA. buxfoliathat have cNMP binding domain absent in them according to the prediction of HMMER database. But taking into account the domains predicted in these proteins by the SMART database the cNMP binding domain was present in all of them. Prediction results of the HMMER, SMART, and CDD database demonstrate the absence of the Ion trans domain in AbuCNGC1.4 (Supplementary Figure S3）.

PtCNGC proteins that have cNMP binding domain absent in them include PtCNGC1.1, PtCNGC1.3, PtCNGC1.4, PtCNGC1.5, PtCNGC2.1, PtCNGC2.2, PtCNGC2.3, PtCNGC2.4, PtCNGC4.1, PtCNGC7, PtCNGC14, PtCNGC15.1, PtCNGC15.3 and PtCNGC19 as predicted by HMMER database. But the cNMP binding domain was absent only in PtCNGC14 and present in all the other aforementioned PtCNGCs (Supplementary Figure S4）.Details of CNGCs reported in other plants are shown inTable 1.

3.2 Physiochemical properties and subcellular localization analysis of cyclic nucleotide-gated channels in citrusSpp.

The detailed physio-chemical properties of 150 CNGC proteins of five CitrusSpp.are shown in (Table 2）.C. sinensishad protein length ranging from 492–1553aa, molecular weight (MW) ranging from 56.13–177.77 (KDa), and Isoelectric point (PI) ranging from 6.38–9.52, Instability index (II) was above 40 for 25 proteins ofC. sinensis, indicating that most of the proteins were unstable. GRAVY values of 29 proteins ofC. sinensis为负表明大多数防eins were hydrophilic. Results of subcellular localization suggested that all putative CsCNGC proteins were present in the plasma membrane.

The protein length of CreCNGCs ranged from 371–1513aa, molecular weight (MW) ranged from 43.15–173.29 (KDa), Isoelectric point (PI) ranged from 6.33–9.44, 21 CreCNGCs proteins were unstable as they have II above 40. GRAVY values of 24 CreCNGCs were negative indicating that maximum proteins were hydrophilic. Only CreCNGC1.3 was localized in the plasma membrane and nuclear compartments while the rest were found to be localized in the plasma membrane.

CgCNGCs have protein lengths ranging from 299–1289aa, molecular weight (MW) ranging from 34.21–147.32 (KDa), Isoelectric point (PI) ranging from 6.06–9.52, Most of the proteins (22) ofC. grandiswere unstable as these proteins have II greater than 40. GRAVY valuesfor27 proteins ofC. grandiswere negative suggesting that these proteins were hydrophilic. All of the CgCNGCs were found to be localized in the plasma membrane.

A. buxfoliahad protein length ranging from 286–1335aa, molecular weight (MW) ranging from 33.62–766.01 (KDa), Isoelectric point (PI) ranging from 6.02–9.7, Instability index (II) was above 40 for 21 proteins ofA. buxfoliarevealing that most of the proteins were unstable in the test tube. 24 proteins ofA. buxfoliawere hydrophilic as their GRAVY values were negative while 7 proteins were hydrophobic as their GRAVY values were positive. Results of subcellular localization demonstrated that all AbuCNGC proteins were found to be present in the plasma membrane.

P. trifoliata’s protein length ranged from 575–1250aa, molecular weight (MW) ranged from 12.10–143.55 (KDa) forP. trifoliata, Isoelectric point (PI) ranged from 6.62–9.54, Instability index (II) was above 40 for 25 proteins ofP. trifoliatasuggesting that most proteins were unstable. 25 PtCNGC proteins were hydrophilic because GRAVY values for these proteins were negative. ForP. trifoliata PtCNGC2.5was localized in the plasma membrane as well as cytoplasmic and nuclear compartments while the rest were localized in the Plasma membrane. Hence, we can conclude that most of the proteins ofCitrus Spp.were basic, unstable, hydrophilic, and localized in the Plasma membrane. The Citrus CNGC proteins that were stable can be used as a biomarker for further studies.

3.3 Phylogenetic analysis

In total, 20AtCNGCs, 16 OsCNGCs, 12 ZmCNGCs, 21 PbrCNGCs, 15 ZjCNGCs, 32 CsCNGCs, 27 CreCNGCs, 30 CgCNGCs, 31 AbuCNGCs, and 30PtCNGCsgenes were classified into four groups and the fourth group was further classified into two sub-groups, I, II, III, IV-A, IV-B each containing the different number of members. The maximum number of members were present in Group IV (84 members) divided into the clade of Group IV-B with 71 members: two members fromA. thaliana(AtCNGC2and4), three fromO. sativa(OsCNGC2,4aand4b), three fromZ. mays(ZmCNGC10,11and12), five fromP. bretschneideri(PbrCNGC2,4,7,8and9) three fromZ. jujuba(ZjCNGC13,14and15),12fromC. sinensis(CsCNGC4.1-4.6andCsCNGC2.1-2.6), 10 fromC. recticulata(CreCNGC4.1-4.6andCreCNGC2.1-2.4), 12 fromC. grandis(CgCNGC4.1-4.6andCgCNGC2.1-2.6), 9 fromA. buxfolia(AbuCNGC4.1-4.6,AbuCNGC2.1-2.3) and 12 fromP. trifoliata(PtCNGC4.1-4.6andPtCNGC2.1-2.6) and Group IV-A with 13 members: two fromA. thaliana(AtCNGC19and20), two fromO. sativa(OsCNGC19aand19b), one fromZ. mays(ZmCNGC9), two fromP. bretschneideri(PbrCNGC19and20), one fromZ. jujuba(ZjCNGC12), one fromC. sinensis(CsCNGC19), one fromC. recticulata(CreCNGC19), one fromC. grandis(CgCNGC19), one fromA. buxfolia(AbuCNGC19) and one fromP. trifoliata(PtCNGC19）.The minimum number of members present in the clade of group II with 29 members two fromZ. mays(ZmCNGC4andZmCNGC5), three fromO. sativa(OsCNGC5a, OsCNGC5b,andOsCNGC5c), two fromZ. jujube(ZjCNGC4andZjCNGC5), two fromP. bretschneideri(PbrCNGC5andPbrCNGC6), five fromA. thaliana(AtCNGC5, AtCNGC6, AtCNGC7, AtCNGC8,andAtCNGC9), three fromC. sinensis(CsCNGC5, CsCNGC7,andCsCNGC8), three fromC. recticulata(CreCNGC5, CreCNGC7,andCreCNGC8), three fromC. grandis(CgCNGC5, CgCNGC7,andCgCNGC8), three fromA. buxfolia(AbuCNGC5,AbuCNGC7,andAbuCNGC8), three fromP. trifoliata(PtCNGC5, PtCNGC7,andPtCNGC8）.The number of members in other groups was also different as Group I had 66 members and Group III had 56 members in total (Figure 3）.

CNGCsfrom every group shared a clade withArabidopsis CNGCmembers that are a dicot, which demonstrates thatCNGCsemerged after the divergence of monocots and dicots. The close association of CNGC members in CitrusSpp.with AtCNGCs demonstrates that these are orthologs of CNGCs inArabidopsis. Members of the same group might have similar structures and functions. The results of phylogenetic analysis of CNGCs in CitrusSpp.were different than those inA. thaliana, O. sativa, T. aestivum, N. tobaccum, B. oleracea, B. rapa, P. bretschneideri, Z. jujubaas current analysis revealed that Group IV clade was largest with 84 members in total and the clade of group II was smallest with 29 members. The number of members in group IV was almost consistent with the previously reported number of members inZ. mays,在第四组的86名成员。最低麻木吗er of members present in the clade of group I was 25. Overall, the number of members was different in each group as compared to previously reported CNGC members in other plants.

3.4 Gene structure and conserved motif analysis

Gene structure analysis revealed that members from each subspecies are having their own set of exons and introns. Exons that belong to group I ofCsCNGCranged from 6 to 17 while exons that belong to group I ofAtCNGCsranged from 7 to 9. Exons that belong to group II ofCsCNGCwere 7 while exons that belong to group II ofAtCNGCsranged from 6 to 9. Exons that belong to group III ofCsCNGCandAtCNGCranged from 6 to 7. Exons that belong to group IV-A ofCsCNGCwere 12 while exons that belong to group IV-A ofAtCNGCranged from 10 to 11. Exons that belong to group IV-B ofCsCNGCranged from 7 to 14 while exons that belong to group IV-B ofAtCNGCranged from 8 to 9.

Ten motifs were identified in CsCNGCs and named motif 1 to motif 10. Motif 1 represents a combination of the Calmodulin binding motif (CaMB) and motif for the IQ domain. Motif 6 represents the hinge motif, while motif 9 represents the PBC motif. Both these motifs together constitute the cNMP/Cyclic nucleotide-binding domain (CNBD). Other motifs are responsible for unknown functions. The gene structure and logo of conserved motifs ofC. Sinensisare given in (Figure 4,Supplementary Figure S5）.

Group I ofCreCNGCscontained exons ranging from 4 to 17, 5 to 7 exons exist in group II ofCreCNGC, exons that exist in group III ofCreCNGCswere 6–8, and 12 exons were present in group IV-A ofCreCNGCsand 7 to 13 exons were present in group IV-B ofCreCNGCs. Motif 5 represents the PBC motif and motif 2 contains the hinge motif, CaMB motif, and motif for the IQ domain. Other motifs were representing motifs of unknown function (Supplementary Figures S6A,B）.

Group I ofCgCNGCcontained 4 to 17 exons, group II ofCgCNGCcontained 7 exons, group III ofCgCNGCcontained 6 to 8, and group IV-A ofCgCNGCcontained 12 exons and group IV-B ofCgCNGCcontained 7 to 14 exons. Motif 3 represents a combination of CaMB motif and motif for IQ-domain, motif 5 represents hinge region motif and motif 7 represents PBC motif (Supplementary Figures S7A,B）.

InAbuCNGCexons of group I were ranging from 4 to 18, exons of group II were ranging from 6 to 7, exons of group III were ranging from 4 to 7, exons of group IV-A 12, and exons of group IV-B were ranging from 7 to 14. Motif 3 represents a combination of CaMB motif and motif for IQ-domain, motif 8 represents PBC motif and motif 2 contains hinge motif (Supplementary Figures S8A,B）.

Exon number for group I ofPtCNGCranged from 6 to 16, exon number for group II ofPtCNGCranged from 7 to 12, exon number for group III of PtCNGCranged from 6 to 9, and exon number for group IV-A ofPtCNGCwere 12 and exon number for group IV-B of PtCNGCranged from 7 to 18. Motif 2 represents the CaMB motif and motif for the IQ domain, motif 3 represents the Cyclic nucleotide-binding domain that contains both PBC and hinge motif. The representation of motifs and logo of conserved motifs ofP. trfoliatais displayed (Supplementary Figures S9A,B）.

Hence, the PBC motif, hinge motif, CaMB motif, and motif for the IQ domain was conserved in 5 CitrusSpp.indicating that genes identified in the current study are trulyCNGC genes.

3.5 Chromosomal mapping

InC. sinensis32 genes were distributed unevenly on 8 out of 9 chromosomes. C. sinensishad maximum genes (10) at chromosome 9, minimum genes (2) at chromosomes 3, 6, and 8, and there was no gene on chromosome 7. The distribution ofCsCNGCon chromosomes is given in (Figure 5）.

27 genes were mapped unevenly at 8 out of 9 chromosomes inC. recticulata.InC. recticulatamaximum genes (8) were present at chromosome 9, minimum genes (1) were present at chromosome 3, and no gene was present at chromosome 7 (Supplementary Figure S10）.InC. grandischromosome 9 carried maximum genes (8), chromosome 3 carried minimum genes (1), and none of the genes was present on chromosome 7 (Supplementary Figure S11）.InA. buxfoliachromosome 9 contained maximum genes (8), chromosome 3 contained minimum genes (1) and none of the genes was present on chromosome 7 (Supplementary Figure S12）.InP. trifoliatathere were maximum genes (7) present at chromosome 9 and chromosome 1, there were minimum genes (2) present at chromosomes 2 and 8 and there was no gene present at chromosomes 4 and 5 (Supplementary Figure S13）.Thus, it can be inferred thatCNGCgenes were distributed unevenly at 8 out of 9 chromosomes in CitrusSpp.except forP. trifoliatain which genes were distributed at 7 out of 9 chromosomes.

3.6 Gene duplication events

The duplication pairs resulting from segmental duplication inC. sinensisincludeCsCNGC2.3/CsCNGC2.4,CsCNGC2.3/CsCNGC2.5, CsCNGC2.3/CsCNGC2.6, CsCNGC2.4/CsCNGC2.5, CsCNGC2.4/CsCNGC2.6, CsCNGC2.5/CSCNGC2.6. The gene pairs that were tandemly duplicated inC. sinensisincludeCsCNGC1.1/CsCNGC1.8, CsCNGC1.8/CsCNGC1.9, CsCNGC2.1/CsCNGC2.2, CsCNGC7/CsCNGC8. The gene pairs ofA. thalianaCNGCs that were tandemly duplicated includeAtCNGC3/AtCNGC11, AtCNGC7/AtCNGC8, AtCNGC11/AtCNGC12, AtCNGC19/AtCNGC20. The gene pairs ofA. thalianaCNGCs that were segmentally duplicated includeAtCNGC3/AtCNGC13, AtCNGC5/AtCNGC8, AtCNGC6/AtCNGC7, AtCNGC6/AtCNGC9, AtCNGC10/AtCNGC13, AtCNGC14/AtCNGC17. Genes were duplicated segmentally as well as tandemly in bothC. sinensisandA. thalianaindicating that both segmental and tandem duplications are involved in the expansion ofCsCNGC基因。此外,非同义substi的速率tutions (Ka), rate of synonymous substitutions (Ks), Ka/Ks, and duplication time (MYA) were calculated. The Ks of 6 segmental duplicates inC. sinensisranged from 0.0118 to 0.8374, also Ks of 4 tandem duplicates ranged from 0.2099 to 2.0937, and duplication time of both segmental and tandem duplicates ranged from 0.89 MYA to 159 MYA. The Ka/Ks value ofCsCNGC1.1/CsCNGC1.8, CsCNGC2.1/CsCNGC2.2, CsCNGC2.3/CsCNGC2.4,CsCNGC2.3/CsCNGC2.6, CsCNGC2.4/CsCNGC2.5, CsCNGC2.4/CsCNGC2.6, CsCNGC2.5/CSCNGC2.6was less than 1 indicating the occurrence of purifying selection in duplication of these genes. The Ka/Ks value ofCsCNGC1.8/CsCNGC1.9, CsCNGC2.3/CsCNGC2.5, CsCNGC7/CsCNGC8was greater than 1 indicating the role of positive selection in duplication of these genes. Similarly, the Ks of 6 segmental duplicates inA. thalianaranged from 0.2735 to 1.1835, and also the Ks value of 4 tandem duplicates ranged from 0.0504 to 0.9547, and the duplication time of both segmental and tandem duplicates ranged from 3.84 MYA to 90.20 million years ago (MYA) (Table 3）.

InC. recticulatagene pairs that were the product of segmental duplication includeCreCNGC2.2/CreCNGC2.4, CreCNGC2.3/CreCNGC2.4.The gene pairs that were the product of tandem duplication includeCreCNGC2.2/CreCNGC2.3, CreCNGC7/CreCNGC8, and CreCNGC15.1/CreCNGC15.2.Altogether 5 gene pairs were duplicated and among these 3 gene pairs were tandemly duplicated indicating the role of tandem duplication in the expansion ofCreCNGC基因。The Ks of 2 segmental duplicates inC. recticulatawere 0.50, also Ks of 3 tandem duplicates ranged from 0.05 to 0.56, and the duplication time of both segmental and tandem duplicates ranged from 4.23 MYA to 42.74 MYA. The Ka/Ks value ofCreCNGCswas less than 1 indicating that purifying selection has occurred in this duplication event (Supplementary Table S3）.

The segmentally duplicated gene pairs ofC. grandisincludeCgCNGC2.3/CgCNGC2.5, CgCNGC2.3/CgCNGC2.6, CgCNGC2.4/CgCNGC2.5, CgCNGC2.4/CgCNGC2.6, CgCNGC2.5/CgCNGC2.6.The tandemly duplicated gene pairs includeCgCNGC1.5/CgCNGC1.6, CgCNGC2.1/CgCNGC2.2, CgCNGC2.3/CgCNGC2.4, CgCNGC7/CgCNGC8, CgCNGC15.1/CgCNGC15.2.Overall, 10 gene pairs were duplicated and out of these 5 gene pairs were segmentally duplicated and 5 were tandemly duplicated indicating the equal contribution of both events in the expansion ofCgCNGC基因。The Ks of 5 segmental duplicates inC. grandis5的范围从0.03到0.74,也Ks串联duplicates ranged from 0.02 to 0.44, and the duplication time of both segmental and tandem duplicates ranged from 1.71 MYA to 69.20 MYA. The Ka/Ks value of five gene pairs was less than 1 indicating the role of purifying selection in the duplication of these genes. The Ka/Ks value of four gene pairs was less than 1 indicating the role of purifying selection in the duplication of these gene pairs and one gene pair (CgCNGC2.3/CgCNGC2.4) greater than 1 indicating the role of positive selection in the duplication of this gene pair (Supplementary Table S3）.

The gene pairs that were segmentally duplicated inA. buxfoliaincludeAbuCNGC1.2/CreCNGC1.7, AbuCNGC1.2/AbuCNGC1.9, AbuCNGC1.3/AbuCNGC1.7, AbuCNGC1.3/AbuCNGC1.9, AbuCNGC2.2/AbuCNGC2.3.Tandemly duplicated gene pairs includeAbuCNGC1.1/AbuCNGC1.8, AbuCNGC1.1/AbuCNGC10, AbuCNGC1.8/AbuCNGC1.9, AbuCNGC7/AbuCNGC8.In total 9 gene pairs were duplicated and among these 5 gene pairs were segmentally duplicated and 4 were tandemly duplicated indicating the role of segmental duplication in the expansion ofAbuCNGC基因。The Ks of 5 segmental duplicates inA. buxfoliaranged from 0.85 to 2.55, also Ks of 4 tandem duplicates ranged from 0.14 to 1.85, and the duplication time of both segmental and tandem duplicates ranged from 11.36 MYA to 195 MYA. The Ka/Ks value of eight gene pairs was less than 1 indicating the role of purifying selection in the duplication of these gene pairs. While the Ka/Ks value of only one gene pair (AbuCNGC7/AbuCNGC8) was greater than 1 indicating the role of positive selection in the duplication of this gene pair (Supplementary Table S3）.

The duplicated gene pairs that arise from segmental duplication inP. trifoliataincludePtCNGC2.2/PtCNGC2.3, PtCNGC2.2/PtCNGC2.4, PtCNGC2.5/PtCNGC2.6, PtCNGC5/PtCNGC8.The gene pairs that arise from tandem duplication includePtCNGC1.4/PtCNGC1.5, PtCNGC2.1/PtCNGC2.3, PtCNGC2.1/PtCNGC2.4, PtCNGC2.3/PtCNGC2.4, PtCNGC14/PtCNGC1, PtCNGC15.1/PtCNGC15.2.总共有10个基因对重复的和them, 4 gene pairs were segmentally duplicated and 6 were tandemly duplicated indicating the role of tandem duplication in the expansion ofPtCNGC基因。此外,非同义substi的速率tutions (Ka), rate of synonymous substitutions (Ks), Ka/Ks, and duplication time (MYA) were calculated. The Ks of 4 segmental duplicates inP. trifoliataranged from 1.09 to 1.83, also Ks of 6 tandem duplicates ranged from 0.06 to 1.45, and the duplication time of both segmental and tandem duplicates ranged from 4.81 MYA to 111.25 MYA. Mostly gene pairs have Ka/Ks value of less than 1 indicating the role of purifying selection in the duplication of these gene pairs. While the Ka/Ks value ofPtCNGC1.4/PtCNGC1.5was greater than 1 indicating the role of positive selection in the duplication of this gene pair (Supplementary Table S3）.

3.7 Cis-regulatory elements/promoter analysis of citrusSpp.

To clearly understand the role of cis-regulatory elements (CREs) inCsCNGCs, CreCNGCs, CgCNGCs, AbuCNGCs, and PtCNGCs,and the cis-elements in 2 kb upstream of TSS were identified. The results suggested that cis-elements of four types were identified namely, hormone-responsive, light-responsive, stress-related cis-elements, and plant development-related cis-elements inCsCNGCs, CreCNGCs, CgCNGCs, AbuCNGCs, PtCNGCs.

It was observed that cis-elements responsible for light responsiveness were present abundantly inCsCNGCs.Overall, 24 cis-elements responsible for light responsiveness were determined out of which Box 4 element was present in 31CsCNGCs,GT1-motif and G box elements were present in 22CsCNGCsand 23CsCNGCs很少和其他人在场CNGCs.Among 11 hormone-related cis-elements, the ABRE element was present in 22CsCNGCs,the CGTCA motif and TGACG motifs were present in 20CsCNGCs,and 21CsCNGCs,TCA element, and TATC box were present in 9CsCNGCsand11 CsCNGCsand other hormone-related elements were present in very fewCsCNGCs.Among 5 stress-related cis-elements, MBS element (drought inducible) was present in 16CsCNGCs,TC-rich repeats element (defense responsive) was present in 13CsCNGCsand LTR, GC motif, WUN motif was present in very fewCsCNGCs.Out of 8 development-related cis-elements, GCN4_motif and circadian were present in 7CsCNGCs andO2 site element was present in9 CsCNGCs很少和其他人在场CsCNGCs.The results demonstrate thatCsCNGCsare involved in plant growth, development, and response to abiotic stress. The graphical representation of the location and types of cis-elements present inCsCNGCsis given in (Figure 6）.

Cis-elements responsible for light responsiveness were present abundantly inCreCNGCs. Overall, 27 cis-elements responsible for light responsiveness were determined out of which Box4 was present in 25CreCNGCs,G box was present in 22CreCNGCs, andGT1 motif was present in 17CreCNGCs很少和其他人在场CreCNGCs. Among 10 hormone-related cis-elements, ABRE was present in 18CreCNGCs,TGACG motif was present in 17CreCNGCs,TCA element was present in 14CreCNGCs很少和其他人在场CreCNGCs.Among 4 stress-related cis-elements, LTR was present in 12CreCNGCs, MBS was present in 10CreCNGCs, TC-rich repeats element was present in 9CreCNGCs, andGCmotif was present in 4CreCNGCs.Out of 6 development-related cis-elements RY element, O2 site, and GCN4 motif were present in 6CreCNGCs很少和其他人在场CreCNGCs(Supplementary Figure S14）.

CgCNGCsalso contained a number of light-responsive CREs. Overall, 24 cis-elements responsible for light responsiveness were determined out of which Box4 was present in 28CgCNGCs,G box was present in 25CgCNGCs, and GT1 motif was present in 22CgCNGCs很少和其他人在场CgCNGCs.Among 9 hormone-related cis elements ABRE was present in 23CgCNGCs,CGTCA motif was present in 17CgCNGCs,TGACG motif was present in 18CgCNGCs很少和其他人在场CgCNGCs.Among 5 stress-related cis-elements, MBS was present in 18CgCNGCs, TC rich repeats element was present in 13CgCNGCs,LTR was present in 11CgCNGCs很少和其他人在场CgCNGCs.Out of 7 development-related cis-elements, circadian was present in 5CgCNGCs,GCN4 motif was present in 5CgCNGCs很少和其他人在场CgCNGCs(Supplementary Figure S15）.

Cis-elements responsible for light responsiveness were present abundantly inAbuCNGCs.Overall, 24 cis-elements responsible for light responsiveness were determined out of which G box was present in 26AbuCNGCs,箱4是总统ent in 26AbuCNGCs, andGT1 motif was present in 21AbuCNGCs,TCT motif was present in 21AbuCNGCs很少和其他人在场AbuCNGCs.Among 9 hormone-related cis elements ABRE was present in 25AbuCNGCs,TGACG motif and CGTCA motif were present in 16AbuCNGCs很少和其他人在场AbuCNGCs.Among 4 stress-related cis-elements, MBS was present in 17AbuCNGCs,TC-rich repeats element was present in 10AbuCNGCs,很少和其他人在场AbuCNGCs.Among 6 development-related cis-elements, circadian was present in 7AbuCNGCs,O2 site was present in5 AbuCNGCs很少和其他人在场AbuCNGCs(Supplementary Figure S16）.

Cis-elements responsible for light responsiveness were present abundantly inPtCNGCs.Overall, 24 cis-elements responsible for light responsiveness were determined out of which Box 4 was present in 25PtCNGCs,G box was present in 21PtCNGCs,and GT1 motif was present in 18PtCNGCs很少和其他人在场PtCNGCs. Among 10 hormone-related cis-elements, ABRE was present in 22PtCNGCs,CGTCA motif and TGACG motif were present in 17PtCNGCs很少和其他人在场PtCNGCs.Among 5 stress-related cis-elements, MBS was present in 17PtCNGCs,TC-rich repeats element was present in 11PtCNGCs,LTR was present in 7PtCNGCsand others are present in very fewPtCNGCs. Among 7 development-related cis-elements, the RY element was present in 5PtCNGCs,Circadian was present in 4PtCNGCs很少和其他人在场PtCNGCs(Supplementary Figure S17）.

3.8 Pan-genome wide investigation of miRNAs targetingCsCNGCgenes, protein-protein interaction, and gene ontology enrichment analysis

A total of 226 miRNAs were identified that targeted 32CsCNGCswith expectation values ranging from 3.5 to 5 (Figure 7A）.Only 1 miRNA was targetingCsCNGC17with an expectation value of 3.5, while 16 miRNAs were targetingCsCNGC7where all miRNAs have expectation value 5 except Csi-miRN925 with expectation value 4.5, 10 miRNAs were targetingCsCNGC1.1,CsCNGC2.3andCsCNGC2.4, 4 miRNAs were targetingCsCNGC1.2, 9 miRNAs were targetingCsCNGC1.3, CsCNGC10andCsCNGC13, 11 miRNAs were targetingCsCNGC1.4andCsCNGC1.8,5 miRNAs were targetingCsCNGC1.5,CsCNGC15.1andCsCNGC2.1, 7 miRNAs were targetingCsCNGC1.6andCsCNGC18, 13 miRNAs were targetingCsCNGC1.7, CsCNGC16andCsCNGC8, 3 miRNAs were targetingCsCNGC1.9, CsCNGC15.2, CsCNGC2.6, CsCNGC4.3andCsCNGC5, 2 miRNAs were targetingCsCNGC4.1, CsCNGC4.2andCsCNGC19,6 miRNAs were targetingCsCNGC2.2,8 miRNAs were targetingCsCNGC2.5, CsCNGC4.4andCsCNGC4.5, 6 miRNAs were targetingCsCNGC4.6. Detailed information related to these miRNAs regulatedCsCNGCsis given in (Supplementary Table S4）.Among these miRNAs, most of them were responsible for inhibiting the cleavage of target transcript while only a few were involved in inhibiting the translation of target genes.

The PPI network of CsCNGC proteins was constructed to reveal the interaction among these proteins and related proteins (Figure 7B) to understand their degree of connectivity and ultimately their functional relativity. It has been shown that the highest degree of connectivity was shown by syntaxin-121, a protein from theC. sinensisplant, which suggests that this protein may have some functional connectivity with the CNGC proteins. Similarly, other proteins including Membrin-11 and some vesicle-associated membrane proteins (acc: XP_006479311.1) also showed a higher degree of interaction. Among CNGC members, CsCNGC4.6, CsCNGC4.2, and CsCNGC4.3 had higher interactions with other CNGC members as well as other related proteins. CsCNGC2.3, CsCNGC4.1, CsCNGC4.4, CsCNGC16, CsCNGC18 and CsCNGC19 had relatively lesser interactions. This level of connectivity reveals that these proteins might be involved in similar pathways thus regulating particular reactions and performing similar functions.

GO enrichment analyses were carried out on 5 CitrusSpp.to increase our understanding of the dynamic roles of CNGCs genes at the molecular level (Figure 7C;Supplementary Table S5）.Based on GO analysis genes are classified into three major categories: biological process (BP), cellular component (CC), and molecular function (MF). Genes were mostly related to biological processes (4), molecular functions (4), and then cellular components (1). In the biological process, category 147 out of 150 genes were involved in ion transport (GO:0,006,811) and transmembrane transport (GO:0,055,085), 59 genes were involved in potassium ion transport (GO:0,006,813), and only 2 genes were involved in trehalose biosynthetic process (GO:0,005,992). In the cellular component category, 147 genes were mainly found in the membrane (GO:0,016,020) which is consistent with the subcellular localization prediction result. In the molecular function category, 147 genes were involved in ion channel activity (GO:0,005,216), 59 genes out of 150 are involved in voltage-gated potassium channel activity (GO:0,005,249), 32 genes in protein binding (GO:0,005,515), and only 2 genes were involved in the trehalose-phosphatase activity (GO:0,004,805).

3.9 Expression profiling ofC. sinensisunder drought stress

RNA-Seq data analysis was performed for leaves sample ofC. sinensisunder drought stress in two cultivars namely Newhall navel (NHE) orange, and Gannanzao (GNZ) navel orange at 0,10 and 20 days. The results suggest thatCsCNGC2.1andCsCNGC1.4were highly up-regulated in cultivar I (20 days) andCsCNGC4.2(10 days).CsCNGC1.3, CsCNGC15.1, CsCNGC15.2, CsCNGC16andCsCNGC18were slightly up-regulated in cultivar II (Figure 8）.

CsCNGC1.1was slightly up-regulated in cultivar I (20 days) and cultivar II.CsCNGC2.2was slightly up-regulated in both cultivars.CsCNGC1.2was slightly up-regulated in cultivar I (20 days).CsCNGC1.7was slightly up-regulated in cultivar I (10 days).CsCNGC2.2was slightly up-regulated in both cultivars.CsCNGC2.3was slightly up-regulated in cultivar I.CsCNGC2.4was slightly up-regulated in cultivar I (10 days). The increase in the expression level of these duplicated genes suggests that they not only evolved in number but in function also.CsCNGC4.3 and CsCNGC4.4were slightly up-regulated in cultivar I (10 days) and cultivar II.CsCNGC4.5was slightly up-regulated in cultivar II.CsCNGC7andCsCNGC8were slightly upregulated in cultivar II.CsCsCNGC4.6was highly down-regulated in cultivar I (20 days) and slightly down-regulated in cultivar I (10 days).CsCNGC1.4, CsCNGC1.5, CsCNGC1.9, CsCNGC2.1, CsCNGC2.3, CsCNGC2.4, CsCNGC2.5, CsCNGC2.6, CsCNGC4.2andCsCNGC5were slightly down-regulated in cultivar II.CsCNGC1.8was slightly down-regulated in cultivar I (20 days). Most of these genes have evolved through duplication which suggests that they are involved in the stress modulating process either by upregulating or downregulating their expression level. This change in their expression level may contribute to modulating stress response in drought stress as well.CsCNGC10was slightly down-regulated in cultivar I (10 days) and cultivar II.CsCNGC13was slightly down-regulated in both cultivars.CsCNGC15.1, CsCNGC15.2, CsCNGC16was slightly down-regulated in cultivar I.CsCNGC19was slightly downregulated in cultivar I (20 days).CsCNGC1.1andCsCNGC1.2in cultivar I (10 days),CsCNGC1.3andCsCNGC1.5in cultivar I,CsCNGC1.6in cultivar I (20 days) and cultivar II,CsCNGC1.7in cultivar I at (20 days) and cultivar II,CsCNGC1.8in cultivar I (10 days) and cultivar II,CsCNGC1.9in cultivar I,CsCNGC10in cultivar I (20 days),CsCNGC17in both cultivars, andCsCNGC18in cultivar I,CsCNGC19in cultivar I (10 days) and cultivar II,CsCNGC2.1in cultivar I (10 days),CsCNGC2.4in cultivar I (20 days),CsCNGC2.6in cultivar II,CsCNGC4.1, CsCNGC4.2, CsCNGC4.3 and CsCNGC4.4in cultivar I (20 days),CsCNGC7in cultivar I (10 days) andCsCNGC8in cultivar I were those genes that have no change in expression after providing stress condition (Figure 8）.

3.10 Expression validation of the citrus cyclic nucleotide-gated channel genes through quantitative reverse transcription-polymerase chain reaction

To explore the role and relationship betweenCNGCgenes and drought stress, the citrus plant was treated with drought stress under different conditions (Figure 9）.The results showed the expression level of different genes under no treatment and drought treatment at 10 and 20 days. According to the qRT-PCR results, theCsCNGC1.4gene had higher expression after 10 and 20 days of drought treatment compared to the expression level when no stress was applied (Figure 9）.The same pattern of gene expression was observed for other members includingCsCNGC2.1, CsCNGC2.3,CsCNGC2.4,CsCNGC4.2,CsCNGC4.3,CsCNGC4.4andCsCNGC5. The level of gene expression increased after 10 days of treatment and further increased after 20 days of treatment.CsCNGC4.6有不同的表达模式,水平of gene expression under controlled conditions was higher. Drought treatment for 10 days decreased the level of gene expression, while the level of gene expression was again increased after 20 days of drought stress but still lesser than the controlled condition. Two unique genesCsCNGC13andPtCNGC14had the same expression pattern being lesser expression under controlled conditions while increased after treatment with drought stress. Results suggest that these members of theCNGCgene family were sensitive to stress conditions, and thus are involved in stress regulation.

3.11 3D Structure prediction of CNGCs in citrus spp.

The protein structures that were predicted are having almost similar structures except for CsCNGC1.4 and PtCNGC14 which had unique structures (Figure 10）.Three-dimensional structures of solely thirteen CNGC proteins were predicted because these were differentially expressed proteins. Predicted structures of all CNGC proteins were visualized in the interactive 1 preset of Pymol (Yuan et al., 2017) where different colors are used to represent alpha helices and beta sheets. Each CNGC protein contained alpha helices and beta sheets. The long spirals were representing alpha helices while wide arrows were representing beta sheets. The templates used by tRrosetta for modeling the structure of CsCNGC1.4 were 5VA1, 7NP4, 5U6O, and 6UQF. CsCNGC1.4 had 55 alpha helices, CsCNGC2.3 had 38 alpha helices, CsCNGC2.4 had 37 alpha helices, and CsCNGC4.3 had 18 alpha helices, CsCNGC4.6 had 41 alpha helices, PtCNGC14 had 24 alpha helices and PtCNGC13 had 28 alpha helices. While CsCNGC2.1 and CsCNGC4.4 contained 27 alpha helices, CsCNGC13 and AbuCNGC13 contained 26 alpha helices. CsCNGC1.4 had 14 beta sheets, PtCNGC14 had 2 beta sheets, CsCNGC2.1, CsCNGC13, and AbuCNGC13 contained 8 beta sheets while the rest contained 10 beta sheets. The predicted structures of all these CNGC proteins were almost similar except for CsCNGC1.4 and PtCNGC14 suggesting that these proteins are potentially functionally similar too.

4 Discussion

TheCNGCfamily is characterized by the presence of a CNBD domain and 6 TM domains along with a pore region (Saand et al., 2015）.In the present study, theCNGCgene family is reported inC. sinensis, C. recticulata, C. grandis, A. buxfolia, and P. trifoliata. The presence of Ion trans and CNBD domain inC. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliataconfirm the genes identified are trueCNGC基因。Most of the proteins inB. oleracea(Kakar et al., 2017), O. sativa(Nawaz et al., 2014), Z. mays(Hao and Qiao, 2018), T. aestivum(Guo et al., 2018), Z. jujuba mill(Wang et al., 2020）.were basic, unstable, hydrophilic, and localized to the plasma membrane and similar results were found for citrus spp. in the present study. The localization of citrus CNGC proteins to the plasma membrane means that these are ion channel proteins and are involved in the uptake of calcium across the membrane. Pangenome-wide analysis provides a comprehensive overview of diversity at the genomic level involving multiple species, which may lead to the identification of unique genes which are present in specific species instead of being present in all genomes under study (Tahir ul Qamar et al., 2020）.Similarly, in this study two unique genes were identified includingCNGC13andCNGC14.The function of these members has not been yet identified inA. thaliana. Although, the function of these two members has been identified inO. sativa(Xu et al., 2017;Cui et al., 2020）.The number of members inC. grandisandP. trifoliatais the same as the number of members inB. rapawhile the number of members in CitrusSpp.is higher than that inZ. mays(12) (Hao and Qiao, 2018),Z. jujuba(15) (Wang et al., 2020),O. sativa(16) (Nawaz et al., 2014),S. lycopersicum(18) (Saand et al., 2015),A. thaliana(20) (Mäser et al., 2001),P. bretschneideri(21) (Chen et al., 2015),B. oleracea(26) (Kakar et al., 2017),and lower than that inN. tobacum(35) (Nawaz et al., 2018),T. aestivum(47) (Guo et al., 2018）.The phylogenetic analysis classified theCNGC家庭成员分为四大组和两个子任务groups, I, II, III, IV-A, and IV-B that were the same asA. thalianabut some members were missing in CitrusSpp.The members that belong to the same group could have similar structures and functions. Group members inC. sinensis, C. recticulata, C. grandis, A. buxfolia, and P. trifoliatawere named by the phylogenetic relationships withCNGCmembers ofA. thaliana. However,CNGC1.1-1.5andCNGC10were present in group I of CitrusSpp.WhileCNGC13which belongs to the same group was present only inC. sinensis, A. buxfolia,andP. trifoliata.CNGC5, CNGC7,andCNGC8were present in group II of CitrusSpp. CNGC15.1-15.2, CNGC17,andCNGC18were present in group III of CitrusSpp.WhileCNGC14belongs to the same groupandwas only present inP. trifoliata, CNGC15.3also belongs to the same group and was present inC. recticulata, C. grandis, A. buxfolia,andP. trifoliataexcept forC. sinensis.CNGC16also belongs to the same group and was present inC. sinensis, C. recticulata, C. grandis, andP. trifoliatawhileCNGC16.1andCNGC16.2were present inA. buxfolia.CNGC19was present in group IV-A of CitrusSpp.,whileCNGC2.1-2.3andCNGC4.1-4.6were present in Group IV-B of CitrusSpp. CNGC2.4, CNGC2.5,andCNGC2.6also belong to the same group whereCNGC2.4was present inC. sinensis, C. recticulata, C. grandis,andP. trifoliataexceptA. buxfolia, CNGC2.5andCNGC2.6were present inC. sinensis, C. grandisandP. trifoliataexceptC. recticulataandA. buxfolia. In the current study Group IV constituted the largest clade with 84 members while the clade of group II was the smallest with 29 members. While, inA. thaliana(Mäser et al., 2001) clade of group I was the largest and the clade of group IV was the smallest. InB. rapa(Li et al., 2019) group I constituted largest clade and clade of group IV-B was smallest. InZ. mays(Hao and Qiao, 2018) clade of group IV-B was largest and clade of group I was smallest. InB. oleracea(Kakar et al., 2017) clade of group IV was largest and clade of group II was smallest. InP. bretschneideri(Chen et al., 2015) clade of group I was largest and clade of group II and IV-A was smallest. InO. sativaclade of group III was largest and group II was smallest.

Results of chromosomal mapping suggested that most of the genes were present on chromosome 9 inC. sinensis, C. recticulata, C. grandis, A. buxfolia,and on chromosome 1 inP. trifoliata.Minimum genes were present in chromosomes 3, 6, and 8 inC. sinensis,chromosome 3 inC. recticulata,C. grandis, and A. buxfoliawhile chromosome 2 and chromosome 8 inP. trifoliata.The distribution ofCNGCgenes on chromosomes in CitrusSpp.was different as compared to other plants in which the gene family is already reported includingB. oleracea(Kakar et al., 2017) in which maximum genes were present on chromosome 1 and 5 and minimum genes were present on chromosome 7.B. rapa(Li et al., 2019) in which maximum genes were present on chromosome 1 and minimum genes were present on chromosomes 6, 7, and 9,P. bretschneideri(Chen et al., 2015) in which maximum genes were present on chromosomes 1, 8, and 15, and minimum genes were present at 2, 9, 13, 16 and 17,N. tobaccum(Saand et al., 2015) in which chromosome 1 and 8 carried maximum genes and minimum genes were present at chromosome 22 and 11. The gene structures ofC. sinensis, C. recticulata, C. grandis, A. buxfolia, and P. trifoliatawere somewhat similar toA. thalianaas the number of exons and introns of Citrus plants that are being studied were not exactly same asA. thaliana,O. sativa,and other plants. Conserved motif analysis suggested that motifs for IQ domain, CaM binding motif, and CNB motifs were present inC. sinensis, C. recticulata, C. grandis, A. buxfolia, P. trifoliataas reported inB. oleracea(Kakar et al., 2017) in which all the above-mentioned motifs were present.Z. jujube(Wang et al., 2020) also had a similar pattern of motifs. Others includeN. tobaccum(Nawaz et al., 2018) in which CNB motif CaM binding motif and motif for IQ domain were present, andT. aestivum(Guo et al., 2018) in which Cyclic nucleotide binding motif and motif for IQ domain were present. InZ. mays(Hao and Qiao, 2018) motif 3 was the combination of both CaMB and motif for the IQ domain, while motif 4 was the CNB domain and motifs 1, 2, 5, 8, 9, and 10 were transmembrane domains. The motifs were closely related to CNGC motifs inZ. mays. Cis-regulatory elements (CREs) that were present in promoter regions of CitrusSpp.were mainly of four types light responsive, stress-related, hormone-related, and development related. InZ. mays(Hao and Qiao, 2018) hormones, stress, and development-related cis-regulatory elements were present.O. sativa(Nawaz et al., 2014), Z. jujuba(Wang et al., 2020),andN. tobaccum(Nawaz et al., 2018) also contained all these stress-responsive elements. Cis-regulatory element analysis shows that theCNGCgene family is involved in plant response to light, hormone, and abiotic Gene duplication mainly contributes to the expansion of a gene family in plant species. InP. bretschneiderimainly segmental duplication has played role in the expansion of theCNGCgene family (Chen et al., 2015）.In ArabidopsisCNGCsboth segmental and tandem duplications contributed to the expansion of the CNGC gene family. Similarly, both segmental and tandem duplications played a role in the expansion of theCNGCgene family inC. sinensis, C. recticulata, C. grandis, A. buxfolia,andP. trifoliata. InO. sativa(Nawaz et al., 2014) three gene pairs were found to be segmentally duplicated includingOsCNGC1/OsCNGC2, OsCNGC10/OsCNGC11, OsCNGC15/OsCNGC16,and one gene pair was found to be tandemly duplicated includingOsCNGC2/OsCNGC3. Hence, both tandem and segmental duplications contributed to the expansion of theCNGCgene family inO. sativa(Nawaz et al., 2014）.InN. tobacumtheCNGCgene family was also considered to be expanded through both segmental and tandem duplications (Nawaz et al., 2018）.Most of theOsCNGCswere upregulated under abscisic acid treatment (ABA) i.e., 12 and indole acetic acid (IAA) treatment i.e., 11, and very few genes were upregulated under kinetin (KN) i.e., 2 and ethylene (ETH) treatment i.e., 6, where genes belonging to same groups showed similar expression patterns. Under cold stressOsCNGCsthat were present in phylogenetic groups I, II, and III were upregulated and those present in group IV were downregulated whereOsCNGC6exhibited the highest expression andOsCNGC16exhibited the lowest expression. Under pathogen stress where two phytopathogens were inoculated with 4 weeks old rice seedlings includingPseudomonas fuscovaginaeandXanthomonas oryzae pv. oryzae (Xoo)the expression patterns ofOsCNGCsdemonstrated that exceptOsCNGC5andOsCNGC6all otherOsCNGCswere up-regulated underXoowhile all the fourteenOsCNGCswere significantly up-regulated underP. fuscovagineinoculation. Thus, all the OsCNGCs that were duplicated were exhibiting similar expression patterns alongside relevance in their functions.OsCNGC1andOsCNGC2were duplicated genes and were also exhibiting similar expression patterns under abiotic stress i.e., Abscisic acid (ABA) and indole acetic acid (IAA) treatment and pathogenic stress that demonstrates that their functions were overlapping (Nawaz et al., 2014）.The 10 duplicated gene pairs inC. sinensisexhibit similar expression patterns exceptCsCNGC2.1/CsCNGC2.2whereCsCNGC2.1was highly up-regulated in cultivar I at 20 days drought stress and slightly down-regulated in cultivar II while that as not true forCsCNGC2.2. Among 10 duplicated gene pairsCsCNGC1.8was slightly up-regulated in cultivar I at 20 days drought stress,CsCNGC1.9was slightly down-regulated in cultivar II,CsCNGC2.3,CsCNGC2.4andCsCNGC2.5were slightly down-regulated in cultivar II,CsCNGC2.4was slightly up-regulated in cultivar I at 10 days drought stress,CsCNGC2.6was slightly down-regulated in cultivar I at 20 days drought stress whileCsCNGC7andCsCNGC8in both cultivars and aforementioned genes in remaining cultivars were having no change in expression. Thus, we can hypothesize that duplicated genes exhibit similar expression patterns and function overlapping in CitrusSpp.too. It seems that some evolutionary events such as duplication could affect the members of CNGC gene family. On the other hand, mutations in the structure, including upstream/downstream site and coding sequence site of members could change the expression levels of CNGC genes (Abdullah et al., 2021;Faraji et al., 2021;Heidari et al., 2021）.InT. aestivum(Guo et al., 2018),O. sativa(Nawaz et al., 2014), A. thaliana(Mäser et al., 2001), P. bretschneideri(Chen et al., 2015), Z. mays(Hao and Qiao, 2018), Z. jujuba(Wang et al., 2020),andS. lycopersicum(Saand et al., 2015) the CNGC family members were different indicating that gene duplications and gene losses have played an important role in the creation of new genes and functions. The increase in the number ofCNGCgene family members was an important event that contributed to the ability of these plants to adapt to changing environmental conditions.

The miRNAs are non-coding RNAs that regulate gene expression. In this study, a total of 226 putative miRNAs were identified that targeted 32CsCNGCs. Several miRNAs were targeting each gene exceptCsCNGC17which was targeted by a single miRNA andCsCNGC7was targeted by 16 miRNAs. InB. oleracea14 miRNAs were identified that targeted 17BoCNGCs(Kakar et al., 2017).After eliminating false positives based on a threshold value of 5 there remained 5 miRNAs that targeted 9BoCNGCs.Out of these miRNAs, bol-miR838days had five target genes while the rest of them were targeting only one gene. The majority of the miRNAs were related to cleavage while only two miRNAs were involved in the inhibition of translation of target genes. InN. tobacum162 tobacco miRNAs were identified that targeted 18NtabCNGCs(Nawaz et al., 2018).After eliminating false positives based on a threshold value of 4 there remained 79 miRNAs. While, after applying a threshold value of 3 there remained 6 miRNAs from 3 families that comprised 8NtabCNGCs. Most of the genes were having target sites for multiple miRNAs exceptNtabCNGC19which contained the target site of a single miRNA. Prior studies support the evidence that miRNAs are involved in stress response and adaptation including topping and wounding inN. tobacumand miRNAs are also involved in drought signaling in rice (Root, 2016）.The study done by (AAB et al., 2019) demonstrates a list of drought-tolerant plant crops with the involvement of genes of specific gene families and the role of their respective miRNAs. Hence, we can conclude that miRNAs in CsCNGCs will also be involved in their response to drought stress. PPI network analysis showed the interaction among citrus CNGC proteins as well as with the other citrus proteins. Higher connectivity was shown by CNGC and other genes which shows their involvement in pathways. The PPI results performed on BoCNGC proteins show that these proteins also have higher connectivity among themselves and with other proteins suggesting their integrated role in biotic, abiotic stress, and hyper-sensitivity resistance (Kakar et al., 2017）.In maize, the PPI network analysis was conducted based on interactions found on STRING. Similarly, the ZmCNGC proteins also showed connectivity within the CNGC members as well as with the homologous proteins fromArabidopsis(Hao and Qiao, 2018）.在棉花、功能交互分析demonstrated that most of the GhCNGC proteins were found to have higher connectivity with a receptor kinase present in the plasma membrane, FLS2 that activates immune signaling. Several other proteins were showing interactions with RSTK, MOL, and TAD3 which are involved in growth and developmental functions (Zhao et al., 2022）.These results regarding interactions of CNGC family members show their contribution of these genes to the functional as well as regulatory diversity in plants and might be helpful in future research to better understand the functions ofCNGC基因。As CNGCs are ion channels, so according to GO enrichment these genes are present in the plasma membrane, act as transmembrane ion transporters, and are involved in ion channel activity, potassium and calcium ion transport activity, and protein binding activity. InBrassica oleracea, according to biological processes, the BoCNGCs are associated with ion channel activity for transmembrane transport, negative regulation of defense responses, salicylic acid biosynthesis, responses to chitin, and plant-type hypersensitive responses. BoCNGCs are present in the plasma membrane and participate in cellular activities related to transduction, binding, and transport (Kakar et al., 2017）.

Expression patterns ofCsCNGCin leaves samples under drought stress at 10 and 20 days indicated that three genes namelyCsCNGC1.4, CsCNGC2.1,andCsCNGC4.2were highly up-regulated whileCsCNGC4.6was highly down-regulated. Out of two unique genes identified in this study, one is present inC. sinensis,CsCNGC13.The expression analysis of this gene in two cultivars is down-regulated under drought stress which shows some specialty in terms of abiotic stress regulation. These results were similar to the ones demonstrated by earlier studies such as expression patterns ofN. tobaccumshowed that 18CNGCgenes (NtabCNGC2, 3, 5–7, 14, 16–21,and29–34) were up-regulated under Calmodulin stress, 16CNGCgenes (NtabCNGC1, 3–7, 14, 16, 17, 26–28,and30–33) under drought stress and 10CNGCgenes (NtabCNGC2, 3, 5–7, 14, 16, 17, 19and20) under cold stress and some genes were downregulated in response to these stresses (Nawaz et al., 2018）.Expression patterns ofO. sativademonstrated that 10OsCNGCgenes were up-regulated under cold stress, and group IV members were down-regulated under cold stress (Nawaz et al., 2014）.InZ. jujuba ZjCNGC10,8,2,and 15 were downregulated under cold stress (24 h), andZjCNGC4and 12 were up-regulated under cold stress (1 h). The majority ofZjCNGCswere down-regulated after being treated with salt stress, particularly group III members, and the same was the case forZjCNGCsunder alkaline stress (Wang et al., 2020）.InB. oleracea13BoCNGCsgenes were up-regulated under cold stress. However, moreBoCNGCswere up-regulated under pathogen stress ofXanthomonas campestris pv. campestris(Xcc) as compared to those treated with cold stress (Kakar et al., 2017).Promoter and expression analysis revealed some genes that have variable expression under abiotic stress. It is hypothesized that several hormones and abiotic stress-related elements control the variable expression level ofCsCNGCsunder various abiotic stress conditions. As a result, this study confers that these genes can be used in future research due to their importance in abiotic stress response.

5 Conclusion

In this study, a total of 32 genes inC. sinensis, 27 genes inC. recticulata, 30 genes inC. grandis, 31 genes inA. buxfolia,and 30 inP. trifoliatawere identified as belonging to theCNGCsgene family.CNGC基因被确定基于CNGC-specific图案s and domains.CsCNGCs, CreCNGCs, CgCNGCs, AbuCNGCs, and PtCNGCshave diversity in their functions, protein lengths, and gene structures. Previously, Genome-wide studies have been done on theCNGCgene family in other plants but the present study is illustratating a pangenome-wide representation of theCNGCgene family among five CitrusSpp. To the best of our knowledge, this is the first research implementing the concept of pangenome-wide analysis and will be helpful for further pan-genome wide studies on other plants in the future. This analysis provided a detailed explanation regarding the pattern of evolution ofCNGCsin CitrusSpp.their intron-exon patterns, distribution ofCNGCgenes on chromosomes, prediction ofCNGCspecific motifs and domains, duplication type, along with promoter region analysis indicating which regulatory elements are more likely to influence the expression of particular genes. Phylogenetic analysis revealed thatCNGCsof these five citrus species were clustered into four major groups and two sub-groups. A fewCNGCsin the groups were missing or might be duplicated during evolution. CREs analysis reveals the association of gene families in response to abiotic stresses. The miRNAs also play a role in the response ofCNGCgenes to drought stress alongside regulating the expression of these genes. PPI network analysis also provided insights into their connectivity suggesting their involvement in functional regulation. GO enrichment was executed to understand the functions ofCNGCsat the molecular level. Expression profiling was done on tissue-specific data ofC. sinensisunder drought stress that demonstrates thatCsCNGC1.4,CsCNGC2.1,CsCNGC4.2were highly upregulated andCsCNGC4.6was highly downregulated under drought stress. Unique genesCsCNGC13andPtCNGC14also showed higher expression in drought stress. These genes can be used in further studies to develop stress-resistant crops. One can visualize and understand the genomic diversity among the Citrus species being examined. We have observed significant inter and intra-species diversity of the CNGC gene family members. The diversity observed could be due to differences in sequencing approaches. Therefore, further experiemnts are required to get deep insights.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

KZ and MR drafted the manuscript, prepared illustrations, and discussed the content with the MS and KF. MT conceived this research topic and revised the contents of the manuscript. FA revised the manuscript. AA and MA were involved in the final development of the manuscript and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

The authors are thankful to the Researchers Supporting Project number (RSP 2022R491), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at://www.gosselinpr.com/articles/10.3389/fgene.2022.1034921/full#supplementary-material

References