A NIL FLW29 containing yellow rust resistance gene Yr16 introgressed from a wheat variety ‘Cappelle-Desprez’ was crossed with wheat cultivar PBW343. The recipient parent FLW29 (resistant) and cultivar PBW343 (susceptible) were used for identifying the lncRNAs in response to Pst infection. Two wheat genotypes were inoculated with Pst pathotype 46S119 such that for each cultivar a total of 10 pots (five mock-inoculated and five inoculated) with three biological replicates at three different time durations (12, 48, and 72 hpi) were used for the experiment. Seedlings of both the cultivar were grown in plastic pots (20 × 20 × 20 cm³) at a distance of 1.5 cm (seed to seed) at the Indian Council of Agriculture Research-Indian Institute of Wheat Barley and Research (ICAR-IIWBR), Regional Station, Flowerdale, Shimla. When the plants attained the two-leaf stage (~15 days after sowing), fresh urediniospores of pathotype 46S119 were harvested from the infected wheat plants and suspended in sterile distilled water containing Soltrol (20 mg/100 ml), which helps in pathogen adherence to leaves. The spore suspension was sprayed on seedlings using an atomizer and later kept in the dark at 10°C for 16-h light/8-h dark to maintain the relative high humidity. Leaf samples were collected at three different time periods, i.e., at 12, 48, and 72 h post-inoculation. The collected samples were immediately placed in RNAlater^® and stored at −20°C until use. Four samples each at three time points with a total of 12 sample combinations, viz., resistant inoculated (FLW_T12, FLW_T48, and FLW_T72), resistant control (FLW_C12, FLW_C48, and FLW_C72), susceptible inoculated (PBW_T12, PBW_T48, and PBW_T72), and susceptible control (PBW_C12, PBW_C48, and PBW_C72), each with three biological replications, were pooled to increase the detection accuracy of transcriptome analysis.

Total RNA was isolated from Pst-treated and mock-inoculated leaf samples using the Qiagen RNeasy Mini kit (Qiagen Inc., Valencia, CA, USA) according to the manufacturer’s instructions, including the recommended treatment with DNase. The RNA quality was verified using RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA) on 2100 Bioanalyzer (Agilent Technologies, USA). RNA concentrations were determined with a NanoDrop ND-8000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA; Thermo Scientific, Wilmington, DE, USA). cDNA libraries were constructed using an Illumina-TrueSeq RNA library preparation kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s recommended protocol and sequencing was carried out on single HiSeq 4000 lane using 150-bp paired-end chemistry. The library preparation and sequencing were performed by commercial service providers (NxGenBio Life Sciences, New Delhi, India). Briefly, total RNA was used to purify poly(A) messenger RNA (mRNA) using oligo-dT beads to capture polyA tails and attached magnetic beads were used for two rounds of purification. During the second elution, enriched mRNA was fragmented into 200–500-bp pieces using divalent cations at an elevated temperature (94°C) for 5 min. With the use of SuperScriptII reverse transcriptase (Life Technologies, Inc., Carlsbad, CA, USA) and random primers, the cleaved RNA fragments were transcribed into first-strand cDNA. Fragments were end-repaired and A-tailed after second-strand cDNA synthesis and indexed adapters were ligated. To construct the final cDNA library, the products were purified and enriched by PCR. Libraries were sequenced using the paired-end (100 bp at each end) module of the Illumina HiSeq platform. After sequencing, the samples were demultiplexed and the indexed adapter sequences were trimmed using the CASAVA v1.8.2 software (Illumina Inc.).

The bioinformatics pipeline that was followed to identify the T. aestivum lncRNA transcripts is depicted in Figure 1 . Each fastq file was aligned with the T. aestivum reference genome (RefSeqv1.0; IWGSC, 2018) using Tophat (Trapnell et al., 2009) and the Cufflinks package was used to build the aligned reads (Trapnell et al., 2010). Cuffmerge was used to merge all transcript files into a single non-redundant transcriptome (Trapnell et al., 2010). The FASTA file was extracted from the combined GTF file using the GFF reads module of the Cufflinks package. Transcripts having lengths of <200 bp were removed by length filters using Bioperl code. These filtered transcripts were assessed using Coding Potential Calculator2 (CPC2) (Kang et al., 2017) and PLEK (predictor of long non-coding RNAs and messenger RNAs based on an improved k-mer scheme) (Li et al., 2014a) to calculate the coding potential, and only non-coding transcripts were retained for further analysis. The obtained transcripts were searched against tRNA, rRNA, snRNA, and snoRNA databases to filter out housekeeping genes using BLASTN. Further, in order to eliminate the coding transcripts, BLASTX was performed against the wheat protein downloaded from the UniProt database.

Identified lncRNAs were checked for their differential expression pattern between NILs under control and Pst-infected plants at three different time points. The calculation of values of expression in the form of fragments per kilo of transcript per million reads mapped (FPKM) was performed using the Cufflinks tool. Further, these expression values were compared between PBW343 and FLW29 lines for both control and treated conditions at each time point by the Cuffdiff tool (Trapnell et al., 2012) in the form of log(fold change) value. The lncRNAs that met the criteria of a threshold value of absolute log₂ fold change ≥2 (upregulated) and ≤−2 (downregulated) and p-value ≤0.001 were considered as differentially expressed lncRNAs.

Simple sequence repeats (SSRs) are microsatellite markers that are crucial for molecular characterization and provide useful insights into plant genetic diversity (Misganaw and Abera, 2017). They are co-dominant, highly variable, and uniformly dispersed over the whole genome (Pinto et al., 2006). For the identified putative lncRNA, the Krait tool (Du et al., 2018) was used to determine the frequency and distribution of SSRs (mono, di, tri, tetra, penta, and hexa). The default parameters of the Krait tool with respect to frequency of repeats were employed for predicting SSRs, viz., 10, 7, 5, 4, 4, and 43 for mono, di, tri, tetra, penta, and hexa nucleotide repeats, respectively.

To identify the lncRNA functions as a precursor of miRNA, the 122 precursor sequences of known 119 miRNAs were downloaded from the miRBase database (http://microrna.sanger.ac.uk/) (Griffiths-Jones et al., 2006, 2008) and aligned with identified lncRNAs. A lncRNA harboring a miRNA precursor sequence with 100% query coverage and similarity was considered a precursor of that miRNA. The hairpin loop formation in lncRNAs was analyzed using the miRNAFold server (https://evryrna.ibisc.univ-evry.fr/miRNAFold) (Tav et al., 2016), and the secondary structure was plotted using the Vienna RNAfold web server (http://rna.tbi.univie.ac.at/) (Gruber et al., 2015).

A total of 119 previously reported mature miRNAs of T. aestivum downloaded from the miRBase database (http://microrna.sanger.ac.uk/) (Griffiths-Jones et al., 2008) were used for interaction analyses with lncRNAs using psRNATarget (https://www.zhaolab.org/psRNATarget//) (Dai and Zhao, 2011). The lncRNAs that showed interaction with miRNAs were then considered for interaction analyses with mRNAs (coding sequence or CDS) of T. aestivum downloaded from National Center for Biotechnology Information (NCBI)) with an expected value threshold of 3.0 and target accessibility set at a maximum of 25. The interaction network of lncRNAs, miRNAs, and mRNAs was visualized using Cytoscape (http://cytoscapeweb.cytoscape.org/) (Shannon et al., 2003) version 3.9.1.

To understand the functional characteristics of the target genes of the identified lncRNAs, blastx was performed against the nr database. This was performed to identify the proteins that had the highest sequence similarity with the given transcripts to retrieve their functional annotations, and a typical cutoff e-value < e−10 was set. BLAST2GO (http://www.blast2go.com/b2ghome) (Conesa et al., 2005) program was used to obtain Gene Ontology (GO) annotations of the differentially expressed lncRNAs (DELs) for describing biological processes, molecular functions, and cellular components. The GO graph of the targeted genes of the identified lncRNAs was generated using the WEGO program (Ye et al., 2018). The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed on the KOBAS website (Bu et al., 2021).

To find the homology of T. aestivum lncRNAs with the previously known and reported lncRNA sequences of other plant species, BLASTn analysis was performed using an e-value of 0.001. Other known lncRNAs of different plant species were downloaded from different lncRNA databases. LncRNAs of Hordium vulgare were downloaded from CANTATAdb (http://cantata.amu.edu.pl/index.html), while O. sativa, Zea mays, and Sorghum bicolor were downloaded from GREENC database (http://greenc.sciencedesigners.com/). Also, to check the homology of the identified lncRNA to the wheat CDS, BLASTn was performed with percent identity >85%, query coverage >80%, and e-value < e⁻⁵⁰.

Wet lab expression analysis of a few Pst-induced differentially expressed lncRNAs was performed using semiquantitative PCR. For this analysis, total RNA from Pst-treated and mock-inoculated leaf samples was isolated at two time points (12 and 72 hpi) following the procedure described in the section, and it was purified using TURBO DNA-free™ Kit (Thermo Fisher Scientific, Waltham, MA, USA) to eliminate any chances of genomic DNA contamination. High-quality RNA was used to synthesize cDNA using the Revert Aid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, USA) according to the manufacturer’s instructions (Thermo Scientific, USA). Primer3 program was used to design gene-specific primers for a total of eight lncRNAs that were found to be differentially expressed between rust susceptible and resistant lines using RNA-seq analysis ( Table S15 ). Additionally, a wheat actin gene segment was amplified as a positive control using the primer pair 5′ CCAAGGGCTGTTTTCCCTAG 3′ and 5′ CTCAAGTACCCGATTGAGCA 3′. Amplification reactions were set up in 20-μl volume containing, 1 U of Taq DNA polymerase, 1× PCR buffer, 250 mM of dNTPs, 100 ng of cDNA, and primers at the concentration of 0.2 μM. The PCR program was as follows: 95°C for 5 min, followed by 35 cycles each consisting of 94°C for 1 min, 55°C–58°C for 1 min and 72°C for 1.5 min, and finally at 72°C for 8 min. Amplified products were separated on 2.0% agarose gel at constant 100 V with 1× Tris acetate EDTA (TAE) buffer (pH −8.0) and visualized in a gel documentation system (AlphaImager, American Instrument Exchange, Haverhill, MA, USA).

The total assembled 164,095 transcripts generated from Cuffmerge was filtered by length using a perl script where 340 transcripts having lengths less than 200 nt were discarded, and a total of 163,755 transcripts were retained. These transcripts were assessed using the CPC tool based on open reading frame (ORF) integrity and p-value. All coding labeled parts were filtered, which had ORF lengths less than 100 and p-value less than 0.5, and after that, 40,316 transcripts were obtained. Likewise, the Plek tool was also used to find out the non-coding RNAs. After filtering the coding RNAs, the remaining 57,200 transcripts were classified as non-coding. Further, the common transcripts IDs from CPC and PLEK results were obtained, and a total of 27,628 transcripts IDs were found to be common in both CPC and PLEK filters. By performing BLASTn against the non-coding RNA databases downloaded from the RNA Central (https://rnacentral.org/), genes like rRNA, tRNA, snoRNA, and snRNA were removed from the remaining transcripts. There were 359, 26, 55, and 2 hits found against rRNA, tRNA, snoRNA, and snRNA databases, respectively. After gradually removing all these hits, 27,186 transcripts were retained. Further, the BLASTX program was run against the UniProt database to remove protein parts from these reads. As a result of BLASTX, 7,031 hits were found and discarded, and a total of 20,155 transcripts were obtained. These transcripts were further filtered for their exon count level, and 13,347 mono-exonic transcripts were removed from subsequent analysis. The remaining 6,807 transcripts having more than one exon were identified as the putative lncRNAs. A total of 977, 696, 329,779, 560, and 540 lncRNAs were identified from samples FLW_C12, FLW_C48, FLW_C72, FLW_T12, FLW_T48, and FLW_T72, respectively. Similarly, in the case of PBW343, a total of 307, 386, 442, 541, 404, and 302 lncRNAs were identified from samples PBW_C12, PBW_C48, PBW_C72, PBW_T12, PBW_T48, and PBW_T72, respectively ( Table 1 ). It was observed that 923 (36.8%), 728 (29.1%), 274 (10.9%), and 503 (20.1%) lncRNAs were uniquely present under FLW control, FLW-treated, PBW control, and PBW-treated conditions, respectively, at 12 hpi. Further, 6 (0.2%) lncRNAs were found to be common in all the conditions at 12 hpi. At 48 hpi, 662 (33.5%), 525 (26.5%), 363 (18.3%), and 377 (19.1%) lncRNAs were unique for FLW control, FLW-treated, PBW control, and PBW-treated conditions, respectively. There were four (0.2%) lncRNAs common in all the conditions at 48 hpi. Similarly, 289 (18.9%), 501 (32.7%), 403 (26.3%), and 272 (17.8%) lncRNAs were unique for FLW control, FLW-treated, PBW control, and PBW-treated conditions, respectively, at 72 hpi ( Figure 1C ). It was also observed that eight (0.5%) lncRNAs were common in all the conditions at 72 hpi. Further, comparing all the common lncRNA at 12, 48, and 72 hpi, it was observed that only four lncRNAs were common in both the lines of both control and treated samples in all the different time points, which indicates that these lncRNAs are expressed under all the conditions.

The length distribution of the lncRNAs showed that the average nucleotide length of the lncRNAs of the FLW29 line was 1007, 973, 986, 966, 938, and 926 bp for C12, T12, C48, T48, C72, and T72, respectively. In the case of PBW343, the average nucleotide length found was 911, 968, 911, 984, 874, and 902 in C12, T12, C48, T48, C72, and T72 conditions, respectively. It was observed that the lengths of lncRNAs ranged from 203 to 5428 bp, with the vast majority having lengths between 600 and 900 bp under different conditions ( Figure S1 ). The identified lncRNAs of both lines showed that the exon count was between 2 and 14, and the maximum number of lncRNAs was bi-exonic followed by exon numbers 3, 4, 5, and so on ( Figure S1 ). In this study, 2,765 lncRNAs had two exons, 1,616 had three exons, 1,097 had four exons, and the rest had between 5 and 14 exons. Chromosomal distribution of the identified lncRNA was visualized using the Circos software, depicted in Figure 1B .

To check the conservation level of identified wheat lncRNAs with other species, blastn of lncRNAs was performed against lncRNAs of other plant species with the parameters, e-value < e⁻¹⁰, coverage > 50%, and percent identity >35%. A total of 217 significant blast hits were found out of which 99 lncRNAs were unique. Therefore, out of 6,807 lncRNAs, 99 were found homologous with the other cereal crops. LncRNAs of H. vulgare were found closer to the wheat lncRNA followed by the lncRNAs of O. sativa, Z. mays, and S. bicolor. The first 15 hits of lncRNA were listed with percent identity, coverage, and e-value ( Table S9 ).

To test the conservation of identified lncRNAs with the protein-coding genes, lncRNAs were tested for homology to wheat CDSs by BLAST with percent of identity >85%, query coverage >80%, and e-value < e⁻⁵⁰. The result showed that a total of 46 lncRNAs were homologous with 178 wheat CDS, among which one lncRNA (TCONS_00009015) was found homologous to 11 CDSs related to a receptor-like kinase (RLK) or more specifically serine/threonine kinase function. The detailed information of these lncRNAs when mapped with the CDSs of bread wheat is presented in Table S8 . This finding is consistent with the literature, which suggests that lncRNAs have very low conservation levels when compared to protein-coding mRNAs (CDS) and are species- and tissue-specific.

The FPKM values generated by Cufflinks were used further for obtaining DELs between the two lines—FLW29 and PBW343—under different conditions using the Cuffdiff software. A total of 10 significant differentially expressed lncRNAs were found, out of which five were upregulated and five were downregulated. In these five upregulated DE-lncRNAs, one lncRNA (TCONS_00163170) was common in C48, C72, T48, and T72, and four DELs were uniquely present, out of which three were from control condition (TCONS_00076516, TCONS_00093548, and TCONS_00100461) and one was from Pst-infected condition (TCONS_00073476). In the case of five downregulated DE-lncRNAs, one was from control conditions and four were from Pst-infected conditions. In these five downregulated DE-lncRNAs, all were uniquely expressed, out of which one was expressed in control (TCONS_00025410) and four in Pst-infected conditions (TCONS_00040012, TCONS_00053873, and TCONS_00066365) ( Table S7 ).

LncRNAs are long RNAs found in the nucleus, nucleolus, and/or cytoplasm that can serve as a precursor for smaller ncRNAs such as snRNAs, snoRNAs, and miRNAs. We explored the lncRNAs acting as precursors of known miRNAs in T. aestivum using the miRNAFold server. A total of 13 lncRNAs were predicted as a precursor of 10 miRNAs, which suggests that those lncRNAs can give rise to the mature miRNAs after being acted upon by nuclease enzymes like dicer and/or drosha ( Table S10 ). Figure 2 shows lncRNA TCONS_00052144 containing the precursor and mature sequences of miRNA tae-MIR9775 ( Figure 2A ) and another lncRNA, TCONS_00075437, acting as a precursor of two miRNAs: tae-MIR1128 and tae-MIR5175 ( Figure 2B ).

MiRNAs are short ncRNAs (18–23 nt) that regulate mRNA expression by binding to the 3′ UTR of protein-coding mRNAs. Depending on the complementarity of the miRNA–mRNA interaction, expression is inhibited or silenced. Full complementarity causes mRNA degradation, which silences the genes, whereas partial complementarity reduces mRNA expression, which downregulates the genes. LncRNAs can sometimes interfere in this process by acting as a miRNA sponge, preventing miRNA–mRNA binding. Identified lncRNAs and known miRNAs of wheat available at the psRNAtarget server were taken for this analysis. The identified wheat lncRNAs were uploaded as target sequences against the available wheat miRNAs to the psRNAtarget server and executed with the parameters of max UPE 25 and expectation ≤ 3. A total of 233 interactions with 199 unique lncRNAs and 65 unique miRNAs of T. aestivum were found ( Table S11 ). Target mRNAs of the identified miRNAs were also found using psRNAtarget by submitting wheat CDS as a target and previously identified 65 miRNAs as small RNAs. A total of 902 miRNA–mRNA interactions were found with 757 distinct mRNA transcripts involved in various functions. Individual lncRNA–miRNA and miRNA–mRNA networks were integrated and visualized using the Cytoscape software ( Figure 3A ), revealing an interaction network of lncRNAs, miRNAs, and mRNAs involving miRNA and target mRNAs that can possibly be interfered with by lncRNA, hence affecting the normal gene regulatory process. Multiple sets of interactions were detected such as a single lncRNA (TCON00045608) interacting with three miRNAs and a single miRNA interacting with several lncRNAs and mRNAs ( Figures 3B–D , Table S11 ). A total of 50 interconnected clusters were present in the network. Further, hub genes were identified by using the application cytoHubba, which is itself a part of the Cytoscape tool. The top 1 hub gene was identified on the basis of 11 different algorithms, viz., MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, and Clustering Coefficient. A total of four hub genes (tae-miR1133, tae-miR1122b-3p, tae-miR1127b-3p, and tae-miR167b) were identified, and their targeted lncRNAs and mRNA are given in Figures 4A–D , respectively.

To investigate the functions of lncRNAs, we analyzed the potential targets of lncRNAs by finding target miRNAs and their interactions with known protein-coding genes or mRNAs of wheat. In our study, 757 target genes of lncRNAs were annotated, which are involved in a variety of metabolic processes like biotic stress tolerance, disease resistance, regulation of cell cycle, and cell morphogenesis. Specifically, we observed many lncRNAs targeting biotic stress response-related genes like serine/threonine-protein kinase, cytochrome P450, NBS-LRR-like protein, rust resistance kinase Lr10-like, putative disease resistance protein RPM1, RPP13, At1g58400, and transcription factors like bZIP, NAC, and MYB, which signify their roles in biotic stress response ( Table 2 ). The functional annotation revealed a predominance of different GO categories for the analyzed target genes. It revealed that a total of 20 biological processes, 13 molecular functions, and two cell components were significantly altered in response to stripe rust. The most significantly enriched biological process includes the GO terms like regulation of cellular (GO:0009987), metabolic (GO:0008152), biological regulation (GO:0065007), response to stimulus (GO:0050896), localization (GO:0051179), and immune system process (GO:0002376). Among molecular functions, GO terms like binding (GO:0005488), ATP-dependent activity (GO:0140657), catalytic activity (GO:0003824), antioxidant activity (GO:0016209), and transporter activity (GO:0005215) were significantly enriched. Further, for the cellular component, the GO terms like cellular anatomical entity process (GO:0110165) and protein-containing complex (GO:0032991) were regulated ( Figure 5 , Table S12 ). With the use of the Kobas tool, the KEGG analysis revealed 27 significantly enriched pathways (p ≤ 0.05) in the targeted genes of the identified lncRNAs, some of which were related to the metabolism of ascorbate and aldarate, carbon fixation in photosynthetic organisms, metabolic pathways, arginine biosynthesis, purine metabolism, and the biosynthesis of secondary metabolites, having key significance in plants in disease response ( Figure 5 , Table S13 ). Moreover, 11 domains with significant importance to biotic stress response have been found in the targeted genes of lncRNA such as NAD-binding domain (IPR006140), serine/threonine-specific protein phosphatase (IPR006186), serine-threonine/tyrosine-protein kinase catalytic domain (IPR001245), ubiquinol–cytochrome c reductase hinge domain (IPR023184), NAD-dependent epimerase/dehydratase (IPR001509), reverse transcriptase zinc-binding domain (IPR026960), protein kinase domain (IPR000719), WRKY domain (IPR003657), and NAC domain (IPR003441) ( Figure S2 ).

All of the lncRNA transcripts found in this research were utilized to find potential microsatellites using Krait v1.1.0, a powerful and fast tool with a user-friendly graphic interface for identifying microsatellites across the genome (Du et al., 2018). A total of 635 lncRNAs were identified as having SSRs out of a total of 6,807 putative lncRNAs. The number distribution of SSRs detected from lncRNAs as mono, di, tri, tetra, penta, and hexa were 35, 180, 314, 75, 21, and 9, respectively ( Figure S3 , Table S14 ). Among the microsatellites, the di-nucleotide (CT) motif was the most frequent (8%), followed by the tri-nucleotide (CCG) pattern (5%).

Of the eight lncRNAs that we analyzed, the expression of only two (TCONS_00025410 and TCONS_00073476) could be detected, and their expression, in general, varied between susceptible and resistant genotypes and over two time points (12 and 72 hpi), although only slightly. At 12 hpi, both mock-inoculated and Pst-inoculated plants of sensitive genotypes showed upregulation of the lncRNA TCONS_00025410. However, this lncRNA was not expressed in mock-inoculated plants but showed upregulation in Pst-inoculated plants at the same time point (12 hpi) in the resistant genotype. At 72 hpi, TCONS_00025410 was downregulated in both sensitive and resistant genotypes and in both mock-inoculated and Pst-inoculated plants. The expression of another lncRNA, TCONS_00073476, was downregulated in mock-inoculated plants of sensitive genotype at 12 hpi. However, it was upregulated in the Pst-inoculated plants of sensitive genotype and both the mock and Pst-inoculated plants of resistant genotypes at the same time point. However, at 72 hpi, this RNA showed downregulation in mock as well as Pst-inoculated plants in both sensitive and resistant genotypes ( Figure 6 ).