Genome-wide identification
The whole-genome screening analysis identified a total of 12,822 NBS-encoding genes across 34 species. The species-based genome-wide identification demonstrated that the number of genes present in each genome was independent of its genome size. The genome size to the number of NBS gene regression (R2) values was observed as below the cutoff value (R2 = 0.015). Interestingly, we observed a correlation of genome size to the number of NBS genes in Gossypium sp. (R2 = 0.669). For instance, G. hirsutum with a genome size of 2.4 Gb has the highest number of NBS genes (708 NBS genes), followed by G. barbadense (622 NBS genes) with a genome size of 2.3 Gb, and G. arboreum, a diploid genome with a 1.7 Gb genome, has 365 NBS genes. Similarly, G. raimondii with a 0.7 Gb genome has 323 NBS genes (Fig. 1, Tables S1 and S2).
Genome-wide identification and classification of NBS genes in Land plants. The length of the bar represents the total number of genes and different colors demonstrate the top 20 subclasses, based on the domain architecture.
Classification of NBS genes in land plants
The Pfam domain analysis identified several decoy domains in addition to the NBS domain. The deep analysis of domain architecture reported several functional protein domains associated with NBS domains (Table S2). So, based on the conserved domains, motifs, and their sites in the primary sequences of NBS proteins, all 12,822 genes were divided into 168 classes concerning the presence /absence and copy number of associated domains (Table S3). Of these classes, the 33rd class (only the NBS domain) has 8591 NBS genes, showing the largest class of NBS superfamily and most of the classes were highly specific with few numbers of NBS genes. So, based on the observation, it can be estimated that the main superfamily of NBS proteins contains only the NBS domain. The species-wise distribution and diversity of identified NBS classes were also interesting. As the lower plants possessed only either one or two specific classes e.g., Selaginella moellendorffii belong to a lycophyte, and possessed only 33th classes (only NBS single domain) with 16 NBS genes. Similarly, Spirodela polyrhiza is a species of duckweed known by the common name common duckmeat and has only two classes 33rd and 60th (NBS and NLR). In contrast, the evolutionary development and adaptation process causes drastic changes in the genome of higher land plants. As moved from simple to complex plants, the NBS classes became more complex. For instance, Hordeum vulgare, barley, a member of the grass family, is a major cereal grain grown in temperate climates globally and has the highest number of NBS classes and several unique NBS domain architectures. Similarly, Vitis vinifera, the common grape vine, a species of flowering plant, native to the Mediterranean region, also has twenty NBS classes. Among all land plants, the most common domain architectures were detected as NBS, NLR, TIR-NBS, TIR-NBS—LRR, and RPW8—NBS domain architectures. However, we also observed highly species-specific classes like Methyltransf_11—NBS, NBS-Glyco_transf_8, TIR-NBS-Lectin_legB-Pkinase, TIR-NBS-RHD3, and NB-LRR-NB-LRR-NBS-Retrotran_gag_2. Overall, out of 168 classes, 94 classes were specific-specific, some classes were family-specific, and were genus-specific like the Gossypium species showed a close relationship regarding the classes. The tetraploid species (G. barbadense and G. hirsutum) of Gossypium sp. have more diversity rather than other diploid species (Fig. 1, Tables S3 and S4).
Orthologs analysis: orthogrouping, duplication, and overlapping,
The orthologs groping analysis of 12,822 NBS in 33 land plants demonstrated that 91.8% of total NBS proteins were assigned orthogroups. However, only 8.2% of NBS proteins were highly unique and could not fit into any orthogroups. A total of 664 orthogroups were predicted and out of these, 347 orthogroups (holding 1135 NBS proteins) were species-specific orthogroups and the mean and median orthogroup size was 18 and 3 genes, respectively (Table S5). The species-based comparative genomics overexposed the genetic makeup of different land plants regarding NBS genes. We identified many common/core genes as well as species-specific unique genes based on sequence divergence. For instance, in the case of Arabidopsis layrata, a total of 209 NBS genes were predicted and categorized into 37 orthogroups (with 2 species-specific orthogroups) covering 206 NBS genes. The remaining three genes did not fit any orthogroups due to their clear sequence divergence as compared to other NBS genes. Similar results were also observed for other land plants. For instance, in the Gossypium sp. the Gar, Gba, Ghi, Ghe, Gth, and Gra did not show any species-specific unique orthogroup or ortholog. But observed some unique NBS genes were categorized as “unassigned orthogroup genes” (Table S6). The number of genes in orthogroups of each species was also assessed and identified that the “OG000” orthogroup possessed the highest number with 2535 NBS proteins and shared all species except some lower plants including Gbi, Mpo, and Smo. In addition, some orthogroups were highly specific to species like OG400 was only present in Al with ALNBS134 and ALNBS135 genes. Similarly, the OG024 was specific to Mtr with more than 15 NBS genes and the OG454 (GhirNBS444 and GhirNBS445), and OG455 (GhirNBS626 and GhirNBS627) were specific to G. hirsutum. Similar unique genes were also identified in other land plants (Tables S7, S8, S9 and S10). The orthogroups duplication event analysis conformed several duplication events during the evolutionary process e.g., the OG000 duplicated 1932 times with 100% confidence and 1600 times with a 50% confidence level. Similarly, other orthogroups were also passed by some duplication events (Tables S11 and S12).
To find the relationships among different species the orthogroups and orthologs overlapping were also assessed and compiled into clusters using dendrogram. The result demonstrated that Osa, Hvu, Zmy, Bra, and Sbi formed one clade, and Sly, Stu, and Sme formed another clade. Similarly, the Gossypium sp. Including Gth, Gtu, Gra, Gba, Ghi, Gar, and Ghe formed one clade. In addition, the Gar (A-genome) and Ghe (A2-genome) shared sister branches, and the Ghi (AD-genome), and Gba (AD2-genome) formed another sister branches showing the highest genome similarity index (Figs. S1, S2, and Table S13).
Phylogenetic tree analysis: duplication events at nodes and gene-based phylogenetic tree
The duplication events at internal and terminal nodes indicated the evolutionary relationship of NBS genes in land plants. The duplication events at terminal nodes demonstrated that Ginko biloba was highly conserved and did not demonstrate any duplication at the terminals as compared to other land plants. Entering the Gossypium genus, the N24 node further duplicated more than 210 times and separated other Gossypium sp. from G. raimondii (Table S14). In conclusion, gene duplication events at internal and terminal nodes demonstrated the expansion of the NBS gene family from lower plants to higher plants. The duplication events also cause the gain and loss of NBS-associated terminal domains, as we observed in the domain architecture analysis. This diversity in NBS proteins might be an adaptation of land plants (Fig. 2).
Evolutionary study of NBS genes in land plants. (A) Basic statistics of orthogroups with abbreviated species names. (B) Gene duplication events at terminal nodes of the phylogenetic tree, and (C) gene duplication events at internal nodes of a phylogenetic tree.
The gene-based phylogenetic tree was divided into 282 sub-clusters and each cluster was correspondence to the orthogroups. The highest number of genes was observed in cluster OG000, which gradually decreased to and decreased to OG282. In addition, some clades were highly species-specific, and some were genus while very few were family specific. Due to the large number of genes in top orthogroups only OG2 was present in the phylogenetic tree (Fig. S3).
RNA-based expression profiling
For a deep understanding of NBS genes in different tissues under different biotic and abiotic stresses, we have selected three plant species including A. thaliana, Z. mays, G. arboreum, and G. hirsutum. Based on the RNA-seq expression data most important putative genes were further filtered for detail study. The tissue-specific expression profiling of NBS genes in A. thaliana generally demonstrated that most of the genes are differentially expressed in leaf, shoot, seedling, flower, and silique. However, very little or negligible expression is observed in pollen, endosperm, embryo, and seed. At the orthogroups (OGs) level, some OGs were highly specific to some tissues like OG6 (At4G33300OG6 and At5G04720 OG6), OG15 (AT5G45490 OG15) had the highest expression in root tissue (Table S15). Similarly, the OG2 (AT1G61300OG2, AT5G63020 OG2, and AT1G61190 OG2) showed similar patterns in leaf, shoot, and seedling. The cladogram among tissues and genes also demonstrated a co-expression network of genes. In the case of Zmy tissues, three major co-expression pattern networks were observed among NBS genes (Fig. 3, Tables S16, S19, S22 and S25).
RNA-seq-based expression profiling of NBS gene in different tissues. The OGs associated with gene IDs represent the orthogroups. (A) A. thaliana, (B) G. hirsutum, (C) Z. mays and (D) G. arboreum.
The RNA-seq data of abiotic stresses included oxidative stress, UV, nutrient deficiency, dark, salt, cold, drought, heat, ozone, and wounding stresses. Generally, it was observed that the OG members mostly co-expressed under various abiotic stresses. In A. thaliana, the NBS genes are most commonly expressed in all stresses except oxidative and UV stresses. Among OGs, the OG6, OG11, and OG1 showed similar expression patterns with the highest expression under most of the stresses. Similarly, the OG6, OG0, and OG351 formed one co-expression clade in the Z. mays plant under various abiotic stresses (Fig. S4, Tables S17, S20 and S23).
The biotic stress expression profiling included different pathogens from aphids to viruses (aphids, nematodes, fungi, bacteria, and viruses). The analyses identified several NBS gene responses under different pathogenic stresses in Arabidopsis, Z. mays, G. arboreum, and G. hirsutum. For instance, the OG1 (AT1G66090OG1, AT4G19520OG1, AT5G41750OG1, and AT1G72900OG1) significantly upregulated under Slerotinia sclerotiorum (Fig. 4A). Similarly, under viral disease, the OG1, OG5, and OG6 showed putative responses during viral infection in Arabidopsis and Z. mays (Fig. 4B). The G. arboreum is naturally resistant to several viral diseases, so we have taken the CLCuD grafted RNA-seq for the identification of the NBS gene’s role in the presence of viruses. We identified several differentially expressed NBS genes in G. arboreum under grafted-CLCuD. The OG2 (Gar06G25610OG2, GarUnG04050OG2), OG6 (Gar11G21680OG6), OG4 (Gar11G229450OG4) and OG115 (Gar10G07690OG115) showed significantly upregulation under CLCuD (Fig. 4C). The Gar06G25610OG2 was further validated through gene silencing approaches in this study (will discuss later) (Fig. S5). As we are so interested in the identification of most putative NBS genes and their role in G. hirsutum under cotton leaf curl disease, which is one of the major challenges in Pakistan, we presented a two-contrast agronomic trait associated accession (Mac7 (a highly CLCuD tolerant G. hirsutum accession developed by USDA, but low production) and Coker312 (highly susceptible, most regenerative accession, commonly used for tissues culturing, also posseted important agronomic traits) RNA-seq data. The results were surprising, we found differential regulation of NBS genes in tolerant accession (Fig. 5, Tables S18, S21, S24 and S26). In summary, we found that the OG0, OG2, OG5, OG6, and OG15 have a highly responsive role in different tissues and under various biotic and abiotic stresses in Arabidopsis, Z. mays, G. arboreum and G. hirsutum.
RNA-seq based expression profiling of NBS gene under different biotic stresses. The OGs associated with gene IDs represent the orthogroups. (A) A. thaliana, (B) Z. mays and (C) G. arboreum.
Expression profiling of G. hirsutum NBS genes in different tissues under various biotic stresses. (A) Cotyledon tissues, (B) leaf tissues, (C) expression of NBS gene in two contrasting accessions, 1; Mac7 (tolerant to CLCuD), 2; Coker312 (highly susceptible to CLCuD) under CLCuD infection in leaf. (D) under various root-associated pathogens in root tissue.
Variants detection in tolerant and susceptible accessions
The genome-wide genetic diversity of Mac7 (G. hirsutum tolerant) and Coker312 (G. hirsutum susceptible) with reference to TM-1 G. hirsutum reference genome, identified several SNPs and InDels in NBS genomics regions of Mac7 (InDels: 2989, SNPs: 3594) and Coker 312 (InDels: 2646, SNPs: 2527). The identified variants were characterized into four impact levels i.e., impact as high (affecting splice-sites, stop and start codons), moderate (non-synonymous), low (synonymous coding/start/stop, start gained), and modifier (upstream, downstream, intergenic, UTR). A comparative study of high-impact SNPs and InDels associated genes between Mac7 and Coker identified several unique and common variants e.g., Coker 312 (7 SNPs, 10 InDels) and Mac 7 (22 SNPs, 26 InDels) unique variants. While 4 genes were common between the two accessions. Similarly, the modifier and moderate and low-impact variants were also observed in the two accessions (Fig. 6, Tables S27, S28 and S29).
Genetic variations of NBS genes between Mac7 and Coker 312. (A) variants distribution region Based on impact level of variants on gene functions. (A) HI; High Impact, (B) MFI; Modifier Impact, (C) MI; Moderate Impact, (D) LI; Low Impact, comparison at SNPs (single nucleotide polymorphisms) and InDels (Insertion and Deletion).
Physiochemistry, gene ontology, and metabolic pathways
The amino acid (AA) sequence analysis of Gossypium hirsutum showed a length ranging from 500 to 1500 amino acid (AA), molecular weight (MW) ranged from 50 to 150 kDa, isoelectric point (pI) value ranged from 2 to 6, charge ranged from + 55.5 to − 44.5 and the Grand Average of Hydropathy ranged from 0 to − 0.2. The gene ontology comprises three components; biological process, molecular function, and cellular components. In the case of molecular function, the results demonstrated that the NBS is mainly involved in the Adenosine diphosphate (ADP) and Adenosine triphosphate (ATP) binding activity with signal transduction. The KEGG pathway analysis indicated the role of the NBS gene in plant-pathogen interaction (PPI) and signaling pathways (SP). The promoter analysis identified several stress-responsive elements in the promoter region of NSB genes. A total of 72 cis-regulatory elements were identified some of the responsive elements were present in all NBS genes like TATA-box, CAAT-box, Box, ARE, G-box, GT1-motif, ABRE, CGTCA-motif, TGACG-motif, TCT-motif, TCA- element, MBS, GATA-motif, CAT-box, and O2-site. Most of the cis-regulatory elements belonged to stress responses (Fig. S5, Table S30).
Molecular docking, and protein–protein interaction
For the proteins modeling, we selected genes based on upregulated in G. arboreum (naturally resistant to CLCuD) under grafted-CLCuD and whitefly induced in Mac7 (G. hirsutum, tolerant to CLCuD) covering OG0, OG2, OG15, and OG43. The selected genes’ translated proteins were used for 3D structure prediction using the I-TASSER server. The PDB database proteins 6J5T-C,6S2P-N, 7CRB-A,4TZH-A, and 4U09-A were used as template sequences. The 3D molecules of ATP and ADP were downloaded from the cheminformatics database. The docking results demonstrated a stable interaction of NBS genes with ATP and ADP with a range of − 7 kcal/mol to 8.2 kcal/mol, except Gar09G25760_OG43 and Ghir_A12G021600_OG15, which showed below − 6.8 kcal/mol and − 6.6 kcal/mol, respectively. The binding affinity of Gar12G23120_OG0 showed that it has a more stable interaction with ATP as compared to the ADP molecule. Similarly, the OG2 genes (Gar01G01860 OG2, Gar06G24920 OG2) showed more affinity to the ATP molecule as compared to the ADP molecule. In contrast to this, the OG0 (Gar11G29700OG0, Ghir_D13G021900OG0) showed more stable interaction with the ADP molecule. The interacting residues of NBS proteins varied from ATP to ADP e.g. the ATP binds with Lys179, His240, Trp139, Val261, Gln 264, Glu142, and Val 284, whereas the ADP binds with Asp143, Glu142, Asn289 and Lys286 in Gar12G23120_OG0 protein. Similarly, other binding results also demonstrated similar patterns (Fig. S6, Table S31). For further detailed molecular interaction of OG2 genes (Gar06G24920_OG2) with CLCuD viral proteins (AC1, AC2, AC3, AC4, AV1, and AV2), we assessed the interaction level by Gibs free energy. The Gar06G24920_OG2 protein demonstrated the highest stability value (− 1352 kcal/mol) with AC1 viral proteins followed by AV2 (− 1103.8), AV1 (− 1026.4), and so on. Similarly, another member of OG2 group protein (Gar06G25610/Gohir.A06G19220/Ghir_A06G020580), also demonstrated high interaction value with AC1 (− 1506.6 kcal/mol) and AV1(− 1255.8 kcal/mol), AC3 (− 1230.8 kcal/mol). The active residues were also changed with the viral protein complexes. (Fig. S7, Table S32).
Silencing of OG2 (G2) genes enhanced virus titer in G. arboreum
TRV-based VIGS assay in cotton is a well-established approach to conducting functional studies of OG2 (Gar06G24920_OG2)genes. FDH-228 plants were inoculated with VIGS constructs to silence the subjected genes, 15 days post infiltration, a completely bleached phenotype was visible on TRV-GrCLA1 inoculated cotton plants, at this time, RT-qPCR based gene silencing was performed (Figs. S8, S9). The results demonstrated a considerable reduction in the expression level of the G2 gene in FDH-228 silenced plants, compared to the TRV:00 inoculated control plants (Fig. 7A). The graft-mediated CLCuD inoculation has been well demonstrated by Ullah et al.58 in G. arboreum plants. Since this breakthrough, it has become possible to identify the resistance imparting genes of G. arboreum against CLCuD. After assessing the gene silencing, an equal number of silenced and control plants were inoculated with CLCuD by grafting CLCuD harboring scions, and exposing them to viruliferous whitefly. 25 days post-CLCuD exposure, enhanced virus concentration was witnessed in G2 (Gar06G24920_OG2) silenced plants in comparison to TRV:00 inoculated controls both under graft and viruliferous whitefly exposure (Fig. 7B). Minor disease symptoms were evident only on graft-inoculated G2 (Gar06G24920_OG2) silenced plants (Fig. 7C). These results suggested the likely involvement of the G2 (Gar06G24920_OG2) gene in CLCuD resistance in FDH-228 plants.
qPCR-based estimation of gene expression and CLCuD Symptoms development on VIGS plants. Panel (A), Orange represents a decreased expression of the G2 gene in G. arboreum plants. Green shows showing expression of G2 in empty vector inoculated control plants. Typical CLCuD symptoms of lower grade appeared on G2 VIGS plants. Panel (B) shows qPCR-based estimation of virus titer in both graft and whitefly exposed VIGS plants of G. arboreum. ** is for a significant difference. Panels (C) is showing minor vein thickening on G2 VIGS plants and TRV:00 inoculated FDH-228 plant leaf with no symptoms. The plants were picturized using a microlens-assisted camera of the Apple iPhone 11 pro max model a2161.
- SEO Powered Content & PR Distribution. Get Amplified Today.
- PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
- PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
- PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
- PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
- Source: https://www.nature.com/articles/s41598-024-62876-5