Comprehensive transcriptomic analysis of two RIL parents with contrasting salt responsiveness identifies polyadenylated and non-polyadenylated flower lncRNAs in chickpea.

Salinity severely affects the yield of chickpea. Understanding the role of lncRNAs can shed light on chickpea salt tolerance mechanisms. However, because lncRNAs are encoded by multiple sites within the genome, their classification to reveal functional versatility at the transcriptional and the post-transcriptional levels is challenging. To address this, we deep sequenced 24 salt-challenged flower transcriptomes from two parental genotypes of a RIL population that significantly differ in salt tolerance ability. The transcriptomes for the first time included 12 polyadenylated and 12 non-polyadenylated RNA libraries to a sequencing depth of ~50 million reads. The ab initio transcriptome assembly comprised ~34 082 transcripts from three biological replicates of salt-tolerant (JG11) and salt-sensitive (ICCV2) flowers. A total of 9419 lncRNAs responding to salt stress were identified, 2345 of which were novel lncRNAs specific to chickpea. The expression of poly(A+) lncRNAs and naturally antisense transcribed RNAs suggest their role in post-transcriptional modification and gene silencing. Notably, 178 differentially expressed lncRNAs were induced in the tolerant genotype but repressed in the sensitive genotype. Co-expression network analysis revealed that the induced lncRNAs interacted with the FLOWERING LOCUS (FLC), chromatin remodelling and DNA methylation genes, thus inducing flowering during salt stress. Furthermore, 26 lncRNAs showed homology with reported lncRNAs such as COOLAIR, IPS1 and AT4, thus confirming the role of chickpea lncRNAs in controlling flowering time as a crucial salt tolerance mechanism in tolerant chickpea genotype. These robust set of differentially expressed lncRNAs provide a deeper insight into the regulatory mechanisms controlled by lncRNAs under salt stress.

Introduction

Soil salinization is a major global concern for sustainable agriculture production and food security (Ahmed et al., 2021; Venkataraman et al., 2021). Globally, 20% of cultivated and 33% of irrigated agricultural land are affected by salinity, and this proportion is increasing, thus causing severe crop yield losses (Kaashyap et al., 2017, 2018). Chickpea is the second most important food crop and Australia is the largest producer worldwide (Bandekar et al., 2022; Wood and Scott, 2021). Given its high nutritional value, low glycaemic index and hypocholesterolemic properties, chickpea is increasingly described as the ‘food crop of the future’ (Havemeier et al., 2017; Singh et al., 2021).

One potential approach for agricultural sustainability and meeting future generations' food requirements is improving the genetic potential of cultivars to make them salt tolerant. To date, breeding efforts have identified QTLs for salt tolerance traits; however, climate variability, the polygenic nature of salt stress and the confounding effects of abiotic stresses in field trials make these QTLs unstable and unable to efficiently reveal the complexity of interacting genes involved in salt tolerance (Atieno et al., 2021; Soren et al., 2020). Furthermore, marker‐assisted selection has limitations in capturing the real‐time expression level of genes regulating salt tolerance, according to the presence or absence of genes involved in salt tolerance mechanisms (Chen et al., 2021; Qin et al., 2020). Gene expression studies have enabled identification of candidate genes associated with salt tolerance; however, these genes are controlled by large and complex gene networks and identifying the genetic function of these networks remains challenging (Ben‐Amar et al., 2016). Additionally, the studies aimed at understanding stress tolerance mechanisms have focused mainly on the coding genome, thus providing an incomplete picture of the transcriptional landscape.

Progress in RNA‐sequencing technology has recently led to a deeper understanding of transcriptome and the identification of rare and weakly expressed noncoding transcriptional units located in the intergenic and overlapping coding regions (Khemka et al., 2016; Nejat and Mantri, 2018). Recent evidence has definitively indicated that the noncoding portion of the genome is widely expressed and is involved in controlling biogenesis and gene regulation (Yuan et al., 2018; Zhang et al., 2021). Long noncoding RNAs (lncRNAs) are emerging as key regulators of the genome that respond primarily to stress and plant development (Gelaw and Sanan‐Mishra, 2021; Huo et al., 2021; Jain et al., 2021). Thousands of lncRNAs have recently been discovered in plants including Arabidopsis (Wang et al., 2021), Oryza sativa (Jain et al., 2021), Glycine max (Zhang et al., 2021), Medicago truncatula (Zhao et al., 2020), Lycopersicum esculentum (Yang et al., 2020), Camellia sinensis (Baruah et al., 2021a; Varshney et al., 2019), Citrus limon (Bordoloi et al., 2022) and Capsicum annuum (Baruah et al., 2021b) through RNA‐sequencing technology. However, lncRNAs are a large collection of highly heterogeneous transcripts with multiple features depending on their site of origin in the genome, sequence, expression level and interaction with neighbouring genes; therefore, their functional characterization remains challenging. The lncRNAs are classified primarily as polyadenylated and non‐polyadenylated, and are further divided into three subcategories according to their genomic location of origin and mode of action: long intergenic RNA (lincRNA), intronic RNA and naturally antisense transcribed (NAT) RNA (Xu et al., 2020; Yin et al., 2021). To elucidate the pathways regulated by lncRNAs in response to stress, and to understand their mode of action and function, these RNAs must be extensively classified and their potential target genes studied.

The lncRNAs have diverse and essential roles in regulating many genes involved in biological processes, such as splicing, histone modification, cellular localization and mRNA processing (Lucero et al., 2021; Rigo et al., 2020). LncRNAs accomplish these functions through controlling the transcription of proximal genes or by acting distally from target loci through regulating signalling pathways (Statello et al., 2021). Only a few studies have performed functional characterization of selected lncRNAs, such as COOLAIR and COLDAIR, which regulate flowering in Arabidopsis through FLOWERING LOCUS C (FLC) gene repression (Marquardt et al., 2014; Sanbonmatsu, 2019; Zhao et al., 2021). Another important lncRNA, npc536, imparts salt tolerance in Arabidopsis and regulates root growth during salt stress (Lu et al., 2016). Furthermore, lncRNAs are crucial for chromatin remodelling, histone modification and epigenetic gene regulation (Bohmdorfer and Wierzbicki, 2015; Chakraborty et al., 2018; Liu et al., 2018; Zhang et al., 2019). Some lncRNAs regulate stress responses and plant development through post‐transcriptional mRNA processing, and interaction with transcription factors and gene promoters (Ulitsky and Bartel, 2013). However, given the versatile functionality and target specificity of lncRNAs, current knowledge represents a small fraction of their vast gene regulatory potential.

Since this area of research is relatively new, and little information is available regarding the functional role of thousands of lncRNAs, studying the genome‐wide and differential expression of lncRNAs in response to various stresses is crucial. Co‐expression gene network analysis allows the molecular interactions of these lncRNAs to be studied on the basis of their expression and k‐mer content, thus holistically revealing the landscape of both the cis‐ and trans‐acting gene regulatory potential of lncRNAs in response to stress (Kaashyap et al., 2022). To date, no report has identified lincRNAs at the post‐transcriptional level and elucidated the molecular pathways through expression network analysis in response to stress in chickpea.

The aim of the present study was to comprehensively classify the polyadenylated and non‐polyadenylated lncRNAs at the post‐transcriptional level to allow characterization of naturally antisense transcribed and intergenic lncRNAs, and elucidate their molecular interaction pathways that control the expression of salt‐tolerant genes in chickpea.

Results and discussion

Overview of mapping and transcriptome assembly

To comprehensively identify lncRNAs according to their sites of origin in the genome, we generated 12 salt‐tolerant and 12 salt‐sensitive flower transcriptome datasets from two types of RNA: polyadenylated RNA [poly(A+)] and ribo‐depleted RNA [poly(A‐)]. RNA‐sequencing libraries from three biological replicates of each control and stress condition were deep sequenced to ~50 million paired‐end reads with insert sizes greater than 300 bp. In total, 1167 million reads after filtering of rRNA reads and quality trimming were mapped to the improved CDC frontier Kabuli v2.6.3 reference genome (Edwards, 2016) (http://doi.org/10.7946/P2G596).

Overall, 864 million reads were mapped to the chickpea genome, with an average 87% concordant pair alignment in the poly(A+) reads and 65% concordant pair alignment in the poly(A‐) reads. An ab initio transcriptome assembly comprised 67 064 transcripts from the poly(A+) RNA (32 982 transcripts) and poly(A‐) RNA (34 082 transcripts). Interestingly, the poly(A‐) assembly had 3% more assembled reads than the poly(A+) assembly.

A total of 17 642 transcripts (9362 transcript loci from the poly(A+) dataset and 8280 transcripts from the poly(A‐) dataset) were mapped to intergenic regions. Of these, 15 044 transcripts longer than 200 bp with coding potential value (CPC ≤ 0) were further scanned with INTERPRO to identify whether any protein signatures were present. Finally, 9419 transcripts with no protein signatures and no small noncoding RNAs, such as snoRNAs, microRNAs (miRNAs) or tRNAs, were identified as lncRNAs.

Distribution of lncRNAs in chickpea chromosomes

The eukaryotic genome comprises higher proportion of noncoding RNA than coding RNA (Palazzo and Lee, 2015). The location of coding genes and the identified lncRNAs were mapped on different chickpea chromosomes. The chromosome distribution of lncRNAs was compared with that of mRNAs. Among eight chromosomes and a scaffold, chromosome 1, chromosome 3 and the scaffold had more lncRNAs than the coding portion of the genome (Figure 1a). Furthermore, chromosome 5 had a nearly comparable number of lncRNAs to the coding part.

Figure 1

(a) Chromosome distribution of lncRNAs compared with mRNAs. Ca1, Ca3 and the scaffold had greater amounts of noncoding RNA than the coding genes. (b) Characteristic features of lincRNAs. The number of exons was primarily one or two for lncRNAs.

Characterization of the basic features of lncRNAs

The transcripts represented both coding and noncoding regions as identified on the basis of several characteristics. The primary characteristic feature of lncRNAs is their transcript length, which differs from that of protein‐coding transcripts (Cabili et al., 2011). The lncRNAs ranged from larger (~6331 bp) to smaller (200 bp) lengths. In contrast, protein‐coding genes had a maximum length of ~14 kb. The mean length of lncRNAs was 701 bp, which was shorter than that of the coding genes (5735 bp). This finding was consistent with previous reports indicating that lncRNAs are shorter than coding genes. For example, the lncRNAs identified in rice have an average length of 800 bp compared with the coding genes, which have an average length greater than 1.5 kb (Zhang et al., 2014a). Another important characteristic of lncRNAs is that they have fewer exons than coding genes (Cabili et al., 2011). For example, the number of exons varied from one to eight: 76% of lncRNAs had a single exon, 0.09% had two exons, 0.08% had three exons and only 0.003% had eight exons (Figure 1b).

These findings are similar to previous studies in humans (Cabili et al., 2011), fish (Pauli et al., 2012), C. elegans (Nam and Bartel, 2012) and plants (Wierzbicki et al., 2021). Furthermore, compared with the poly(A+) lncRNAs, poly(A−) lncRNAs had fewer exons and were more abundantly expressed across the conditions in the two genotypes. LncRNAs may act as enhancers and affect the expression of their neighbouring coding genes (Ørom et al., 2010). Enhancers are mostly non‐polyadenylated, and lncRNAs may have similar modes of action and function in the genome as enhancers (Sun et al., 2020). Insufficient studies have been performed on enhancer characterization in plants, but further studies on poly(A−) lncRNAs may provide elucidation in the future.

Classification and identification of lncRNAs

The lncRNAs were broadly classified into polyadenylated and non‐polyadenylated groups to understand their essential role in mRNA splicing and post‐transcriptional modification. On the basis of the fragments per kilobase of transcripts per million mapped reads (FPKM) values of poly(A+) and poly(A−) lncRNAs, principal component analysis revealed a significant variance of 64% between the salt‐tolerant and salt‐sensitive genotypes under stress conditions. In addition, a variance of 24% was observed between the expression of poly(A+) and poly(A−) lncRNAs in the tolerant versus sensitive genotype under control and salt stress conditions (Figure 2a). This finding suggested that these lncRNAs have specific expression patterns revealing large genetic variations between genotypes and are uniquely expressed in response to salt stress.

Figure 2

(a) Principal component analysis showing the variance between the polyadenylated and non‐polyadenylated lncRNAs expressed across the tolerant and sensitive genotype under salt stress conditions, according to FPKM values. (b) Total numbers of lincRNA and NATs observed after filtering of the transcripts through a stringent pipeline. These lncRNAs were longer than 200 bp, had a CPC score ≤ 0 and did not match the PFAM or RFAM databases. (c) Significantly enriched GO categories for differentially expressed lncRNAs in response to salt stress in chickpea.

In total, we identified 9419 lncRNAs expressed specifically in one genotype × treatment, 2345 of which were novel lncRNAs specific to chickpea in response to salt stress. Among these, 2588 polyadenylated lncRNAs (1008 lncRNAs + 1550 NATs) were expressed in the tolerant genotype, whereas 2931 polyadenylated lncRNAs (1219 lncRNAs + 1712 NATs) were expressed in the sensitive genotype (Figure 2b).

A total of 1901 non‐polyadenylated lncRNAs (695 lncRNAs + 1206 NATs) were expressed in the tolerant genotype, whereas 1999 (767 lncRNAs + 1232 NATs) were expressed in the sensitive genotype. Poly(A−) lincRNAs were weakly expressed, and 93.3% of them had FPKM values between 0.14 and >20, whereas only 73% of poly(A+) lincRNAs had values between 0.2 and >20. Similar results have been reported in maize and cucumber (Hao et al., 2015; Li et al., 2014a). Furthermore, 538 poly(A+) lncRNAs, 790 poly(A+) NATs, 386 poly(A−) lncRNAs and 606 poly(A−) NATs were specifically expressed in the tolerant genotype during stress. In contrast, 565 poly(A+) lncRNAs, 787 poly(A+) NATs, 399 poly(A−) lncRNAs and 622 poly(A−) NATs were uniquely expressed in the sensitive genotype during salt stress.

Notably, 7% more lncRNAs were expressed in the tolerant genotype during stress than control conditions, whereas 14% less lncRNAs were expressed in the sensitive genotype during stress compared with the control condition. The FPKM values for lincRNAs were higher in the stress condition than the control condition in the tolerant genotype. In contrast, lincRNAs were more highly expressed in the control condition than the stress condition in the sensitive genotype. In addition, more NATs than lncRNAs were present in both genotypes. Their abundance in the tolerant genotype compared with sensitive genotype suggested their crucial role in gene silencing during stress. NATs are a diverse and rare class of lncRNAs that regulate several critical biological processes, such as differentiation and development (Villegas and Zaphiropoulos, 2015). Furthermore, important functionally classified lncRNAs, such as COOLAIR, have been found to transcribe from the antisense strand of the genome. These results are consistent with those from a previous study in Oryza sativa, in which more lincRNAs than NATs were identified (Zhang et al., 2014b).

More polyadenylated lncRNAs and NATs were expressed in the stress condition than the control condition in the tolerant genotype. In contrast, fewer polyadenylated lncRNAs and NATs were expressed in stress condition than the control condition in the sensitive genotype. The greater number of poly(A+) RNA than poly(A−) lncRNAs suggested that these lincRNA might have potential role in post‐transcriptional modifications. However, both genotypes showed an increase in the numbers of non‐polyadenylated lncRNAs and NATs under stress. Post‐transcriptional modification is an important mechanism regulating the molecular mechanism of salt tolerance, as further confirmed with Gene Set Enrichment Analysis, in which post‐transcriptional gene silencing (GO:0035194), gene silencing by miRNA (GO:0035195), regulation of gene expression, epigenetic regulation (GO:0040029) and negative regulation of gene expression (GO:0010629) were among the significantly enriched GO categories (Figure 2c).

Differential expression of lncRNAs in response to salt stress

DESeq2 analysis with three biological replicates

Three biological replicates of the control and stress samples from both genotypes were subjected to DESeq2 analysis. The resulting DEGs (FDR < 0.05) were run through a stringent pipeline to filter differentially expressed lncRNAs (DE‐lncRNAs) responsive to salt stress. A total of 178 significant DE‐lncRNAs were identified with a transcript length greater than 200 bp and no PFAM or RFAM match. Of the 178 DE‐lncRNAs, the tolerant genotype had more DE‐lncRNAs (110), whereas the sensitive genotype had fewer DE‐lncRNAs (68). Notably, more lncRNAs were up‐regulated in the tolerant genotype, whereas more lncRNAs were down‐regulated in the sensitive genotype. LncRNAs control the expression of their neighbouring genes (Ulitsky, 2016). Therefore, the up‐regulation of lncRNAs also resulted in the up‐regulation of essential salt tolerance candidate genes in the tolerant genotype, whereas their down‐regulation resulted in the down‐regulation of the salt tolerance candidate genes in the sensitive genotype.

The top induced lncRNAs in the tolerant genotype included XLOC_018822 (FC: 14.13 ↑), XLOC_025093 (FC:11.75 ↑) and XLOC_020890 (FC: 6.23 ↑), and the top repressed lncRNAs in the sensitive genotype were XLOC_028665 (FC: −3.67 ↓), XLOC_028304 (FC: −2.97 ↓) and XLOC_026105 (FC: −1.89 ↓; Figure 3a).

Figure 3

(a) Heatmap showing the differential expression of lncRNAs between the tolerant and the sensitive genotype in response to salt stress. Genotype‐specific differential expression and abundance of (b): poly(A+) and (c): poly(A−) lncRNAs in chickpea.

The lncRNAs were 14‐fold up‐regulated in the tolerant genotype but only 2‐fold up‐regulated in the sensitive genotype. For example, the lncRNA XLOC_018822 was significantly up‐regulated in the tolerant genotype (FC: 14.13 ↑) but down‐regulated in the sensitive genotype (FC: −1.12 ↓; Figure 3b, c; Table 1).

Table 1

List of differentially expressed lncRNAs between the tolerant and the sensitive genotypes in response to salt stress

LncRNAs	Fold change in tolerant	Fold change in sensitive	FDR values	Chromosome	Start	Stop	Origin	Coding/Noncoding	Coding potential
LncRNAs up‐regulated in tolerant but down‐regulated in sensitive genotype
XLOC_018822	14.13	−1.12	6.97E‐28	Ca5	50 597 016	50 597 635	Novel isoform	Noncoding	−1.23
XLOC_025093	11.75	−1.05	1.45E‐23	Ca6	6 408 627	6 409 596	Novel isoform	Noncoding	−1.22
XLOC_020890	6.23	−1.15	6.31E‐12	Ca5	35 720 301	35 721 327	Novel isoform	Noncoding	−1.28
XLOC_020891	2.87	−1.06	0.000266	Ca5	35 711 332	35 712 213	intron	Noncoding	−1.01
XLOC_026939	2.56	−1.06	0.000119	Ca6	58 524 199	58 524 432	Novel isoform	Noncoding	−1.32
XLOC_028665	2.52	−3.67	0.000326	Ca7	23 262 265	23 263 410	Novel isoform	Noncoding	−0.88
XLOC_028304	2.44	−2.97	5.10E‐07	Ca7	12 253 634	12 256 713	Novel isoform	Noncoding	−1.00
XLOC_026105	2.42	−1.89	0.001797	Ca6	32 431 931	32 432 740	Intron	Noncoding	−1.17
XLOC_031230	2.39	−2.00	0.007747	Ca7	38 176 187	38 176 786	Intron	Noncoding	−0.83
XLOC_025223	2.20	−1.04	0.021398	Ca6	9 075 736	9 076 881	Intron	Noncoding	−1.04
XLOC_018777	2.12	−1.16	0.000615	Ca5	48 662 201	48 665 725	Intron	Noncoding	−1.11
XLOC_021514	2.06	−1.02	0.003199	Ca5	60 157 717	60 158 355	Intron	Noncoding	−1.24
LncRNAs up‐regulated in sensitive but down‐regulated in tolerant genotype
XLOC_018173	−2.49	1.00	0.004294	Ca5	27 917 783	27 918 246	Intron	Noncoding	−0.42
XLOC_030070	−2.48	1.73	0.00512	Ca7	4 531 299	4 532 519	Intron	Noncoding	−1.03
XLOC_030033	−2.23	1.05	0.000409	Ca7	3 870 897	3 872 850	Intron	Noncoding	−1.04
XLOC_000376	−2.20	1.88	0.014099	Ca1	6 169 102	6 169 755	Intron	Noncoding	−1.21
XLOC_001991	−2.17	2.82	0.034844	Ca1	2 694 912	2 695 478	Intron	Noncoding	−1.11
XLOC_010817	−2.03	1.26	0.005343	Ca3	47 875 594	47 876 073	Novel isoform	Noncoding	−0.94
XLOC_031002	−1.73	1.49	0.013768	Ca7	31 660 242	31 661 313	Novel isoform	Noncoding	−0.72

These fold changes were confirmed with qRT‐PCR, and the values correlated well with RNA‐sequencing data (Figure S3). Although differential expression of lncRNAs has been reported, this is the first study showing the comparative differential expression patterns of lncRNAs between parental genotypes of an RIL mapping population (salt‐tolerant JG 11 and salt‐sensitive ICCV 2), which segregate for the salt tolerance trait.

Co‐expression network analysis predicting lncRNA functions

Cell wall biogenesis and flower development

Weighted gene co‐expression network analysis (WGCNA) was performed to identify the role of lncRNAs in gene regulatory pathways. The WGCNA co‐expression analysis assembled the lncRNAs and potential gene targets into 41 modules (Figure S1).

Gene modules M22, M23 and M41 in the tolerant genotype, and gene modules M4, M10 and M29 in the sensitive genotype, were significantly up‐regulated during salt stress. The modules in the tolerant genotype were abundant in lncRNAs co‐expressed with genes involved in cell signalling and cell wall biogenesis, and genes‐encoding transcription factors and transporters (Figure 4a). In contrast, modules expressed in the sensitive genotypes had less abundant lncRNAs and consequently less abundant genes involved in cell signalling, cell wall biogenesis and transport.

Figure 4

(a) WGCNA gene modules up‐regulated during salt stress in tolerant and sensitive genotypes. LncRNAs were more abundant in the tolerant genotype than the sensitive genotype; hence, salt‐responsive genes are induced during stress. (b) Gene regulatory network showing the targets of lncRNAs. The red dots show the source genes, and the red arrows show the target genes. ABA induces the expression of lncRNAs controlling the expression of salt‐responsive genes. (c) Synergistic mode of lncRNA co‐expression with the cytochrome P450 gene. The up‐regulation of lncRNAs induces the expression of cytochrome P450 genes in the tolerant genotype while repressing their expression in the sensitive genotype during salt stress.

Interestingly, DE‐lncRNAs were significantly co‐expressed with important genes involved in flowering (FLC), cell wall biogenesis (expansin), cell signalling (cationic peroxidase), transcription (MYB and ERF) and cell transport (Na+/K+ transporters) in response to salt stress. Among the highly expressed genes in these modules, the hub genes were a cation/H+ antiporter, pollen receptor kinase, cytochrome P450, abscisic acid (ABA), flowering locus (FLC) and the transcription factors MYB and ERF. The hub genes of module 22 were the lncRNA Ca27835 (kME: 0.96), peroxidase (Ca31840, kME: 0.99), cation/H+ antiporter (Ca20179; kME: 0.99), pollen receptor‐like kinase (Ca18394, kME: 0.99), agamous‐like MADS box (Ca20058; kME: 0.96) and cytochrome P450 (Ca09424; kME: 0.96). These genes have crucial roles in pollen development, and their role as hub genes in the module up‐regulated in the tolerant genotype suggested that the pathway involved in pollen tube development is crucial for salt tolerance. The co‐expression network analysis revealed that intergenic lncRNAs control the expression of these genes. Notably, these genes are the targets of lncRNAs, but the expression of these lncRNAs are instigated by the abscisic acid signalling gene (Figure 4b). The transcription factors MYB and ERF are the targets of lncRNAs that control their downstream gene expression. However, in most cases, lncRNAs had salt‐responsive genes as their direct targets; therefore, the mode of lncRNA action can be direct or indirect. Transcription factors such as MYB and ERF are targets of lncRNAs and ABA genes. After being triggered by lncRNAs, MYB/ERF genes control the expression of cationic peroxidase, cytochrome P450, expansin and FLC. ERF is an important transcription factor that regulates salt tolerance. Interestingly, ERF also targets MYB gene expression, but the opposite is not true. ABA targets all the important salt tolerance genes in the network, whereas all genes in the regulatory network target cationic peroxidase. Notably, FLC and expansin target cytochrome P450, thus suggesting the importance of cytochrome P450 in regulating pollen development through cell biogenesis mechanisms in response to salt stress.

LncRNAs were more abundant in the tolerant genotype than the sensitive genotype; hence, the expression of the target salt‐responsive genes was up‐regulated. Importantly, lncRNA (Ca26840) was synergistically co‐expressed with cytochrome P450 (Ca09486) gene, which means with the induction of lncRNA, the gene was also induced. Interestingly, 1.5‐fold up‐regulation of this lncRNA induced 100‐fold expression of the cytochrome P450 gene in the tolerant genotype during salt stress (Figure 4c). In contrast, up‐regulation of this lncRNA did not affect the expression of cytochrome P450 in the sensitive genotype in response to salt stress. LncRNAs are highly tissue and genotype specific, and the tolerant genotype might have a different noncoding RNA mode of action from that of the sensitive genotype. This finding may be crucial for understanding the regulation of salt tolerance in chickpea.

LncRNAs control flowering through DNA methylation and chromatin modification

Salinity severely affects the reproductive phase and therefore identifying the genes in flower development and reproductive success is imperative to understand salt tolerance in chickpea. Physiological studies have indicated that the tolerant genotype maintains more flowers to combat salt stress. To establish the role of lncRNAs and identify their potential flowering gene targets during salt stress, the chromatin modification genes that were tightly co‐expressed with the FLC genes responsible for flowering in chickpea were analysed. Chromatin modifications are instrumental in controlling the essential flower development genes in plants. These modifications include DNA methylation and histone modification, which regulate chromatin and flowering gene expression during stress. We identified a strong correlation between the expression values of lncRNAs (Ca15486, Ca32852 and Ca16678) and histone deacetylase and histone‐lysine N‐methyltransferase, thus suggesting the involvement of these lncRNAs in chromatin remodelling during flowering (Figure 5a).

Figure 5

(a) Correlogram showing co‐expression of the lncRNAs with genes involved in chromatin remodelling. Long noncoding RNAs show high Pearson correlation coefficient values (red) with histone modification and DNA methylation genes. (b) Synergistic mode of lncRNA co‐expression with DNA methylation genes. The expression of lncRNAs was higher in the tolerant genotype than the sensitive genotype during stress. (c) Gene co‐expression network showing that Ca10465 (DNA methylation) (red node) is a target of a lncRNA (Ca32852) (blue node top left). Blue nodes are the genes’ neighbours closely co‐expressed in a network. The arrows in red show the targets of the genes, and green dots show the source genes. (d) Synergistic mode of expression of chromatin remodelling genes with the FLC gene. The expression of the chromatin remodelling gene Ca11940 and FLC gene (Ca14012) was twofold higher in the tolerant genotype than the sensitive genotype.

The lncRNA Ca32852 triggers a cascade of DNA (cytosine‐5) methyltransferase (Ca10465) gene signalling in response to salt stress. The DNA (cytosine‐5)‐methyltransferase gene showed unique tissue specificity and differential expression between the tolerant genotype (Ca07939; FC: 13.9 ↑) and the sensitive genotype (Ca07939; FC: 1.32 ↓; Figure 5b). The DNA methylation gene signals the symplekin gene (Ca05304), a critical histone modification gene. Other essential genes co‐expressed with DNA methyltransferase genes included those encoding splicing factors (Ca07011), an enhancer‐binding protein (Ca16836), mechanosensitive ion channels (Ca21004) and cleavage and polyadenylation specificity factor (Ca06873). Although the expression of lncRNAs and DNA methylation genes was twofold higher in the tolerant genotype than in the sensitive genotype, these genes were repressed in the tolerant genotype during stress (Figure 5c). These findings suggest that epigenetic modification may be caused by lncRNAs and their co‐expression with flowering genes.

Chromatin remodelling induces the expression of FLC genes, which suppress flowering by controlling the timing of flower initiation in plants. Interestingly, these genes were repressed in the tolerant genotype but induced in the sensitive genotype during stress, thus suggesting an important salt tolerance mechanism allowing the tolerant genotype to produce and maintain more flowers during salt stress (Figure 5d). Concomitantly, the histone‐lysine N‐methyltransferase gene (Ca03197) was closely co‐expressed with the chromatin and FLC genes. Furthermore, histone‐lysine N‐methyltransferase genes mediating chromatin remodelling are involved in repressing the flowering gene (Kim and Sung, 2017; Kim et al., 2017; Shea et al., 2019). Because salinity delays flowering and leads to flower abortion, thus resulting in severe crop yield losses, these master regulators may be crucial in establishing salt tolerance in chickpea cultivars.

Role of lncRNAs in mRNA splicing

LncRNAs modulate gene expression in response to stress through active involvement in events such as mRNA cleavage and splicing (Rigo et al., 2020; Song et al., 2021). To identify the interactions between lncRNAs and splicing factors, we identified lncRNAs in the co‐expression network. We searched for closely connected neighbour genes to understand the functions of lncRNAs in mRNA splicing. Most genes in the network were associated with splicing and post‐transcriptional/post‐translational modifications. For instance, the lncRNA Ca32852 triggers the expression of genes involved in mRNA splicing (Figure 6a). In contrast, the gene pre‐mRNA splicing factor Ca07011 was the hub gene of this important signalling cascade network. This gene is a target of an ATPase (Ca21731), ATP‐dependent helicase (Ca12310), cullin protein (Ca27354), G protein signalling modulator (Ca21728), intronic splice facilitator (Ca01366) and cleavage/polyadenylation specificity factor (Ca06873) (Figure 6b). Furthermore, the lncRNA‐controlled expression of splicing factors targets important salt tolerance genes, such as a calcium‐transporting ATPase (Ca13196), mechanosensitive ion channel (Ca21004) and zinc metalloprotease. Intriguingly, a pre‐mRNA splicing factor (Ca07011) targets the tRNA (cytosine(34)‐C(5))‐methyltransferase (Ca03738), in agreement with the finding that lncRNAs affects chromatin remodelling by interacting with target genes' specific splicing factor sites (Bardou et al., 2015, Ariel et al., 2015). Furthermore, the expression of different isoforms of genes, such as those encoding the ion channels that span the plasma membrane, was not only mediated by lncRNAs but also specifically regulated in response to stress conditions.

Figure 6

(a) Gene co‐expression network showing that Ca07011 (pre‐mRNA spicing factor) is a target of lncRNA (Ca32852). Blue nodes are the gene neighbours closely co‐expressed in a network. The arrows in red show the targets of the genes, and green dots show the source genes. (b) Gene co‐expression network showing that Ca07011 (pre‐mNA spicing factor) is a hub gene that controls the expression of genes encoding an ATPase, intronic splice factor and mechanosensitive ion channel protein under stress conditions. (c) IGV visualization of lncRNAs’ origins from intergenic regions of the chickpea genome. Differential expression of lncRNAs regulates the expression of a glutathione/chloride channel gene. Ca_v2.6.3_gene.gff3: Cufflink transcript files from tolerant control and stress; Intergenic.bed: intergenic regions of chickpea genome; Tolerant DEGs: differentially expressed genes in the tolerant genotype; Sensitive DEGs: differentially expressed genes in the sensitive genotype.

We further visualized the genomic locations of lncRNAs and coding genes in Integrative Genomics Viewer (IGV_2.3.98). The XLOC_022526 lncRNA spanned the intergenic region of the chickpea genome and had a longer transcript assembled in the tolerant genotype during stress (CUFF.17726) than the control condition (CUFF.15067.1) (Figure 6c). Interestingly, this lncRNA is located ~100 kb downstream from the glutathione/chloride channel gene (Ca20075). The differential expression of lncRNAs instigates the expression of glutathione/chloride channel genes located upstream. Therefore, this gene was up‐regulated in the tolerant genotype but was not expressed in the sensitive genotype during stress. Chloride channel genes have essential role in the efflux of Cl‐ ions. Therefore, lncRNA‐mediated regulation of this gene should be investigated as a critical factor in response to salt stress.

Functional annotation of lncRNAs through orthologous inference

To further confirm the biogenesis and functional role of lncRNAs in the regulation of stress tolerance genes, we performed BLAST analysis of the FASTA sequences of the lncRNAs obtained in this study against the functionally validated lncRNAs from databases (GreeNC, lncRNAdb) of six plant species: (i) Arabidopsis, (ii) Medicago, (iii) Glycine max, (iv) Phaseolus vulgaris, (v) Oryza sativa and (vi) Vitis vinifera. We found strong matches between chickpea lncRNAs and previously functionally validated lncRNAs from Arabidopsis and Medicago species. A total of 1130 poly(A+) lncRNAs were specific to chickpea, whereas the remaining lncRNAs showed at least a 32‐bit score match with other species. Similarly, 1215 poly(A‐) lncRNAs were specific to chickpea, whereas the remainder showed at least a 32‐bit score match with other species. Of the total lncRNAs, 24% matched with Medicago species and 15% matched with Glycine species, whereas only 6% showed BLAST hits with Arabidopsis species (Figure 7). Importantly, these poly(A+) and poly(A−) lncRNAs had specific sequences based on their transcription from the different RNAs and did not share any commonalities.

Figure 7

Comparative analysis of chickpea lncRNA sequences against six model plant species: (1) Arabidopsis; (2) Medicago; (3) Glycine; (4) Phaseolus; (5) Oryza and (6) Vitis.

The lncRNAs showed the greatest synteny with the model legume species, Medicago truncatula and Glycine max. However, interestingly, nearly 11% showed synteny with Vitis vinifera, which is very distantly related to chickpea, thus suggesting that some lncRNAs are evolutionarily conserved. This finding may aid in comparative studies to functionally annotate lncRNAs (Paytuvi‐Gallart et al., 2019). According to lncRNAdb (http://www.lncrnadb.org), we identified only seven functionally annotated lncRNAs in Arabidopsis (IPS1, At4, COOLAIR, COLDAIR, npc536, npc48 and TERRA; Franco‐Zorrilla et al., 2007; Shin et al., 2006; Swiezewski et al., 2009); two in Glycine max (IPS1, alias: TPSI/Mt4 family; At4, alias: TPSI family; Martín et al., 2000); three in Medicago truncatula (IPS1, At4 and ENOD40) (Girard et al., 2003); two in Oryza sativa (IPS1 and ENOD40); one in Vitis vinifera (IPS1) and none functionally characterized in Phaseolus vulgaris. Of those functionally characterized, the chickpea lncRNAs showed strong homology with only three lncRNAs: COOLAIR, IPS1 and At4. These lncRNAs have been described on the basis of overexpression of these transcripts in mutant plants or the addition of an antisense promoter to study the gene silencing mechanisms regulated by these lncRNAs (Huang et al., 2011; Zhao et al., 2021).

Twenty‐one lncRNAs (11 poly(A+) lncRNAs and 10 poly(A−) lncRNAs) showed significant homology with COOLAIR, AT4 and IPS1, which have been well characterized in Arabidopsis species. Two poly(A+) lncRNAs (XLOC_016687 and XLOC_032977) and four poly(A−) lncRNAs (XLOC_004010, XLOC_033371, XLOC_038053 and XLOC_012076) were homologous with COOLAIR lncRNA. COOLAIR is an antisense transcribed lncRNA that initiates at the terminator and terminates at the promoter of a gene, and is known to silence the FLC gene in response to cold stress (Kim et al., 2017; Zhao et al., 2021). Interestingly, XLOC_032977 was 92‐fold up‐regulated in the tolerant genotype compared with the sensitive genotype, but was only 2‐fold up‐regulated in response to salt stress. Similarly, the poly(A‐) lncRNAs were up‐regulated in the tolerant genotype compared with the sensitive genotype (Table 2).

Table 2

Functional annotation details of chickpea lncRNAs on the basis of syntenic relationships with model plant species

LncRNAs	Species	Cicer arietinum	Function	Expression in response to salt stress
IPS1	Oryza sativa	XLOC_016933	Sequester mir‐399	Up‐regulated (FC: 2.1 ↑)
	Brassica rapa
	Medicago truncatula
	Populus tremula
	Medicago sativa
	Glycine max
	Vitis vinifera
	Lycopersicon esculentum
COOLAIR	Arabidopsis thaliana	XLOC_016687	Cold‐induced silencing of FLC gene	Up‐regulated (FC: 5.5 ↑)
At4	Arabidopsis thaliana	XLOC_007084	Phosphate‐induced plant growth	Up‐regulated (FC: 8.0 ↑)
	Glycine max
	Medicago truncatula
TCONS_00046739	Medicago truncatula	XLOC_011795	Salt stress	Up‐regulated (FC: 2.5.0 ↑)

The lncRNA XLOC_016687, a NAT, was highly induced (5.5‐fold) in the tolerant genotype but significantly repressed (28.7‐fold) in the sensitive genotype during salt stress. Previous studies have shown similar results, wherein up‐regulation of COOLAIR lncRNAs results in regulation of FLC (Rosa et al., 2016). Plants undergo several phase changes to adapt to environmental cues (Zhao et al., 2017); the most crucial phase change is from the vegetative to the reproductive stage.

Consequently, along with the up‐regulation of the antisense transcribed lncRNA XLOC_016687, the FLOWERING LOCUS T gene was induced in the tolerant genotype (Ca31297; FC: 8.69 ↑) during the salt response, thus suggesting that antisense transcribed lncRNAs span the coding gene on the sense strand and may be involved in effectively controlling gene expression and/or silencing in response to salt stress.

The miRNAs are key regulators of developmental and physiological processes in plants (Qu et al., 2021; Secic et al., 2021). The lncRNAs At4 and INDUCED BY PHOSPHATE STARVATION1 (IPS1) regulate the phosphate content in root and shoot biomass thereby affecting plant growth (Franco‐Zorrilla et al., 2007). These lncRNAs regulate the post‐transcriptional modification in coding genes through miRNA‐mediated mRNA cleavage (Rojo Arias and Busskamp, 2019). Similarly, the lncRNA IPS1 is expressed in species including Vitis, Glycine, Phaseolus, Oryza and Arabidopsis. The IPS1 lncRNA interacts with mir‐399 and decreases the content of inorganic phosphate (Pi) in the shoots. Both IPS1 and At4 lncRNAs regulate phosphate homeostasis and thus plant growth. Eight poly(A+) lncRNAs and five poly(A−) lncRNAs were found to have significant homology to At4 lncRNA. Importantly, XLOC_007084 was up‐regulated 17.6‐fold in the tolerant genotype but was 5.52‐fold down‐regulated in the sensitive genotype in response to salt stress. In addition, only one lncRNA, XLOC_016933, showed homology with lncRNAs IPS1; this lncRNA was up‐regulated 2.09‐fold in the tolerant genotype during flower development in response to salt stress, thus suggesting that up‐regulation of IPS and At4 lncRNAs may increase the expression of mir‐399, which controls metal homeostasis in response to salt stress in chickpea.

Another important lncRNA, TCONS_00046739, plays a major role during root development in response to salt stress in Medicago truncatula (Wang et al., 2015). Three poly(A+) lncRNAs and two poly(A−) lncRNAs showed significant homology with TCONS_00046739 lncRNA. XLOC_011795 lncRNAs was 2.5‐fold up‐regulated in flowers of the tolerant genotype under salt stress. During flower development, this lncRNA may play a different biological role in inducing successful reproduction events under salt stress. Interestingly, of five lncRNAs showing homology to TCONS_00046739, three lncRNAs, XLOC_037004, XLOC_020713 and XLOC_037004, were found to be NATs. Co‐expression of the TCONS_00046739 lncRNA with the cytochrome P450 gene has been demonstrated to play an essential role in the formation of tapetum walls during pollen development (Pinot and Beisson, 2011; Xu et al., 2020).

Furthermore, the up‐regulation of lncRNA induced the silencing of genes such as FLC, which plays an essential role in flowering time (He, 2012; Menger and Rizvi, 2021; Waseem et al., 2020). Therefore, the up‐regulation of the TCONS_00046739 lncRNA may have a critical role in the overexpression of the cytochrome P450 gene in flower development in response to salt stress. Cytochrome P450 was highly induced in the tolerant genotype (Ca09486; FC: 98.3 ↑) but was repressed in the flowers of the sensitive genotype (Ca00600; FC: −330.8 ↓) in response to salt stress, thus suggesting that the XLOC_037004 lncRNA is involved in overexpression of the cytochrome P450 gene, thus regulating pollen development during salt stress. These results confirmed the K‐means clustering‐based co‐expression analysis of lncRNAs and their respective coding genes. According to the K‐means clustering, cytochrome P450 was closely co‐expressed with XLOC_024019, a lncRNA that was induced in the tolerant genotype (FC: 1.74 ↑) but was highly repressed in the sensitive genotype (FC: −2.28 ↓). These genes require further validation by measurement of lncRNA expression after coding gene overexpression or knockout. Nonetheless, our results indicate the active involvement of these lncRNAs in regulating gene expression.

Materials and methods

Plant material

Two improved chickpea cultivars, JG11 (salt tolerant) and ICCV2 (salt sensitive), were subjected to salt stress in a randomly complete block design at RMIT University, Australia, glasshouse facility. These genotypes are parents of a RIL mapping population and segregate for salt tolerance traits (Khan et al., 2015). Three biological replicates of each genotype were subject to control and salt stress conditions. The seeds were surface sterilized with 70% ethanol, and then rinsed several times with MilliQ water. The seeds were germinated in the dark and moist conditions until radicle emergence.

Three germinated seedlings were sown per 10.5‐inch‐diameter pots filled with 9.5 kg of pasteurized soil and later thinned to one healthy plant. Two adaptive doses of 40 mmNaCl (~1.17 g per kg of soil) were added twice during the life cycle. The first salt dose was given 1 week before sowing, and another dose was given 10 days after sowing and before the first flower stage (Kaashyap et al., 2018). The pots were sealed with sturdy tape to prevent salt leakage. To avoid water retention in pots, the soil field capacity was measured as the amount of water drained per hour per kilogram of soil. Pots were watered and maintained to 80% field capacity throughout the experiment. To monitor salt concentration, we assessed the soil's electrical conductivity and maintained it below ~1 dS/m as chickpea is an intrinsically salt‐sensitive crop (Pushpavalli et al., 2015). The first flower date was recorded, and fully opened flowers were collected when both genotypes in the control and stress conditions reached the first flowering stage. The flowers were snap frozen in liquid nitrogen and stored at −80 °C until RNA extraction.

Isolation of poly (A+) RNA

Flower tissues were ground to a fine powder in liquid nitrogen, and total RNA was isolated with a Qiagen RNeasy kit (GmbH, Germany). From 1 µg of total RNA, poly(A+) RNA was isolated with a Dynabeads mRNA purification kit (Thermo Fisher Scientific). To ensure minimum carryover of the rRNA, we performed extraction of poly(A+) RNA with oligo (dT) beads twice.

Isolation of poly(A−) RNA (RNA devoid of mRNA and rRNA)

After removal of poly(A+) RNA as described above, the total RNA was cleaned with AmpureRNA clean‐up beads and eluted in RNase‐free water. Next, the ribosomal RNA was depleted from this RNA sample with a TruSeq Stranded Total RNA kit with Plant Ribo Zero (Illumina, Inc.). The ribo‐depletion step was repeated twice to ensure no carryover of rRNA in the sample. This process led to the isolation of poly(A−) RNA devoid of poly(A+) RNA and rRNA.

RNA‐sequencing library preparation

The RNA‐sequencing libraries were prepared from the poly(A+) and poly(A−) RNA samples of the two genotypes. These included 12 poly(A+) and 12 poly(A−) RNA samples, each from two genotypes, two conditions (control and stressed) and three biological replicates. A TruSeq stranded library kit (Illumina Inc.) was used for constructing libraries from 100 ng of poly(A+) and poly(A−) RNA samples. The RNA samples were randomly primed and fragmented according to the standardized fragmentation time described to obtain a large cDNA insert size with a median size of 300 bp. The first‐strand libraries were generated with dUTP incorporation in the second strand, which terminated the second‐strand synthesis. The cDNA molecules were uniquely indexed, and six samples were sequenced per HiSeq 3000 (Illumina Inc.) lane, thus generating more than 50 million reads per sample (Figure S2).

Data processing

The RNA sequencing library quality was verified with FastQC. The libraries were screened for several quality parameters: read count, Phred score, length distribution, adapter contamination and RNA contamination. Subsequently, rRNA read contamination was assessed and removed with the sortMeRNA (sortmerna‐intel/2.1) tool (Kopylova et al., 2012) and inbuilt silva rRNA databases. Following this, the reads were trimmed and adapters were removed with the trimmomatic (trimmomatic/0.36) tool (Bolger et al., 2014), which maintained the paired‐ness of the reads.

Mapping and transcriptome assembly

The clean reads were mapped to the improved CDC frontier Kabuli v2.6.3 reference genome (http://doi.org/10.7946/P2G596) with the splice junction aligner tophat (tophat‐gcc/2.0.13) (Trapnell et al., 2009) with default mate‐pair distance parameters. The resulting bam file showing the accepted hits was input into the Tuxedo pipeline, and an ab initio transcriptome assembly was generated with Cufflinks (cufflinks‐gcc/2.2.1). Subsequently, individual Cufflink assemblies from 12 poly(A+) and 12 poly(A−) samples, including stress and control from two genotypes, were merged with Cuffmerge. Cuffmerge assigns class codes to the transcripts to the genome’s gtf file. Finally, the intergenic transcripts denoted with class code ‘u’ were isolated and input into a stringent filter pipeline to identify the lncRNAs.

Filter pipeline for lncRNA

Intergenic transcripts longer than 200 bp were analysed with BLASTx against the SwissProt and UniProt databases. The transcripts were then screened for their coding potential with coding potential calculator (CPC) software (Kong et al., 2007), and only CPC scores ≤ 0 were considered to be noncoding. Next, the transcripts were searched and scanned with InterPro scan to determine whether any protein resemblance or match to protein signatures was present. The transcripts showing PANTHER or PFAM matches were removed. The filtered unique transcripts were matched against RNA families (RFAM) to identify the lncRNAs. Thus, a stringent set of unique lncRNAs was obtained on the basis of the strand information, and characterized as NATs and lincRNAs.

Weighted gene co‐expression network analysis

Co‐expression network modules were obtained with raw gene counts from three biological replicates of each genotype, stress condition and RNA type (poly(A+)/poly(A−)). Raw gene counts obtained from HTSeq analysis were log‐transformed and used as input into the WGCNA package (v1.51) in R (Langfelder and Horvath, 2008). Genes with low counts (<1.00) and low correlation coefficients (<1.00) after log‐transformation were filtered, and the remaining genes were used to construct an adjacency matrix. The co‐expression network modules were attained with the blockwiseModules function and default steps prescribed in the package. A soft threshold power = 16; TOMtype = signed; mergeCutHeight = 0.25; and minModuleSize = 30 were chosen to indicate significance. The expression profile of each module was calculated based on the eigengene value (ME) and gene connectivity (kME) value indicating the correlation strength of an individual gene in each module, as calculated with the signedKME function. Genes with a kME value >0.90 were identified as central hub genes. Finally, a hypergeometric test was performed to assign FDR‐corrected values with the phyper function in the R program (Kaashyap et al., 2022).

Gene regulatory networks

The gene regulatory networks were created with the R Bioconductor packages WGNCA, knitr, limma, ggplots and reshape2. Genes with low counts (<1) were filtered, and the remaining genes were normalized with the function log2 (raw counts +1). A correlation distance matrix was constructed with the function cordist and adjacency.fromSimilarity with a power of 12 and the type selected as signed. A weighted network was developed with a threshold of 0.999. Genes with edges lower than the threshold value or with no edges were filtered out. The non‐positive and negative edges were identified and rescaled to 0 and 1. The adjacency matrix was converted to graphml format with the R package igraph (Kaashyap et al., 2022). The graph consisting of genes with a correlation value represented by the edges of the network was exported as network.graphml and visualized in Cytoscape v 3.8.2 with Prefuse Force Directed layout. Genes showing a high connectivity coefficient were denoted hub genes in the gene network. To identify the source and target genes, edge weights with a significant cut‐off were used.

Validation of differential expression of lncRNAs with qRT‐PCR

The qRT‐PCR primer pairs were designed from 10 significant DE‐lncRNAs among the tolerant and sensitive genotypes. cDNA was synthesized from the flower tissues with random hexamer primers and Superscript III RTase enzyme. Primer pairs specific to the DE‐lncRNA templates were used to analyse the expression patterns of these DE‐lncRNAs across the three biological replicates. In addition, the housekeeping gene EF1 was used as a positive control and for normalization (Figure S3).

Conclusions

The study comprehensively analysed lncRNAs according to their site of origin in genome from two contrasting chickpea genotypes in response to salt stress. More polyadenylated lncRNAs and NATs were differentially expressed in the tolerant genotype than the sensitive genotype during stress. The DE‐lncRNAs were co‐expressed with important flowering genes involved in functions such as chromatin remodelling, DNA methylation and flowering, thus suggesting their role in regulating flower development during salt stress. These genes showed synteny with functionally validated lncRNAs from legume crop species, thus confirming their role in regulating salt tolerance in chickpea. The comprehensive set of lncRNAs identified in this study will benefit the understanding of the complex molecular mechanisms underlying abiotic stress tolerance and in engineering salt stress‐tolerant crops.

Funding

This work was supported by Australia–India Strategic Research Fund (AISRF) Grant number GCF010013, Commonwealth of Australia.

Conflicts of interest

The author(s) declare no competing interests.

Author contributions

Conceptualization, N.M., K.H.M.S. and R.K.V; Data curation, M.K., S.K.; Formal analysis, M.K.; Funding acquisition, N.M., R.F., D.E., K.H.M.S. and R.K.V; Investigation, M.K., R.F. and N.M.; Project administration, N.M. and R.F.; Resources, N.M., S.K. and D.E.; Supervision, N.M., R.F.; Writing – original draft, M.K.; Writing – review & editing, M.K., N.M. All authors have read and agreed to the published version of the manuscript.

Figure S1 Weighted gene co‐expression network analysis (WGCNA) clustered salt‐responsive genes differentially expressed in tolerant and sensitive genotypes. The expression profile of each module was calculated on the basis of the eigengene value (ME) and the gene connectivity (kME) value indicating the correlation strength of an individual gene in each module.

Click here for additional data file.

Figure S2 Mapping statistics of poly(A+) and poly(A−) RNA‐sequencing reads aligned to the chickpea genome. A maximum of 87% concordant pair alignment was observed.

Click here for additional data file.

Figure S3 qRT‐PCR confirmation of the differential expression of lncRNAs observed in RNA‐sequencing data.

Click here for additional data file.