Co-Regulation of Long Non-Coding RNAs with Allele-Speciﬁc Genes in Wheat Responding to Powdery Mildew Infection

: Powdery mildew (caused by Blumeria graminis f. sp. tritici ; Bgt ) is an important fungal disease of wheat ( Triticum aestivum ) worldwide, and results in signiﬁcant crop damage in epidemic years. Understanding resistance mechanisms could have undoubted beneﬁts in controlling disease and minimizing crop losses. The recent explosion in genomic knowledge and the discovery of noncoding RNAs have led to the idea that long ncRNAs (lncRNAs) might be key regulators of protein-coding gene expression. However, in-depth functional analyses of lncRNAs in wheat remain limited. Here, we evaluated the possible role of lncRNAs in regulating functional genes in wheat responding to Bgt pathogen, using genome-wide transcriptome data and quantitative RT-PCR. Our results demonstrated that both long intron ncRNAs (linncRNA) and long intergenic ncRNAs (lincRNAs) play roles in regulating allele-speciﬁc genes, including transcription factors, both positively and negatively. The correlation of expression between lincRNAs and ﬂanking functional genes increased as the spacing distance decreased. Co-expression of microRNAs, their target lncRNA and target functional genes showed that lincRNA interacted competitively with functional genes via miRNA regulation. These results will be beneﬁcial for further dissecting molecular mechanisms of lncRNAs functions at the transcriptional and post-transcriptional levels in wheat.


Introduction
Long noncoding transcripts (lncRNA), defined as a group of RNA transcripts that exceed 200 nt in length with no apparent discernible coding potential, were previously seen as the 'junk' RNAs or 'dark matter' of the genome [1]. However, the recent explosion in genomic knowledge demonstrated that ncRNAs can play roles as key regulators of protein-coding gene expression, either directly or indirectly [2,3], such as competing with endogenous RNAs to regulate miRNA levels [4] and scaffolding ribonuclear protein complexes [5]. Compared with mammalian systems, the functional dissection of plant lncRNAs is still in its infancy. Initial identification of plant lncRNAs was based on bioinformatic searches in cDNA databases for RNAs with poor coding capacity [6,7]. Fortunately, high-resolution maintained by the College of Agronomy of NWAFU. Ten-day-old wheat seedlings were inoculated with Bgt conidia collected from sporulating seedlings of Shaanyou 225 pre-infected 20 days before.

Identifying Functional Genes Adjacent to Differentially Expressed Long Noncoding Transcripts (lncRNAs)
We identified lncRNAs of wheat line N9134, that were regulated in expression pattern after inoculation with Pst and Bgt separately, from the RNA-Seq database obtained in our previous study [14]. Briefly, after all annotated and pathway identified gene were removed, the lncRNAs were identified following four rigorous criteria (transcript length; encoding less than 50 aa; not any transposable elements (TEs); and no gap) as previous described. OrfPredictor was used to identify protein-coding regions in each strand. Differential gene expression analysis was performed with the bioconductor package DESeq, version 3.2. To identify genes neighboring lncRNAs, all 283 assembled DE-lncRNAs were mapped onto the reference genome [29] through alignment with BlastN at p-value < 1.0 × 10 −10 , and the adjacent functional genes were predicted according to the annotation in URGI (Unité de Recherche Génomique Info).

Real-Time Quantitative PCR Analysis
LncRNA, adjacent functional genes and miRNA expression profiles in the contrasting NILs were analyzed by SYBR green-based real-time quantitative PCR (Q-PCR) with cDNA, after cDNA synthesis and RNA extraction from infected leaves sampled at 0, 6,12,24,36,48,72,120,168, and 240 h post-inoculation (hpi). Three independent biological replications were performed for each time point. Q-PCR was performed on a QuantStudio™ 7 Flex Real-Time PCR System (Life Technologies Corporation, USA) with the FastKing RT kit (with gDNase) (TIANGEN, Beijing). Sequence-specific primers of relevant genes (Supplemental Table S1) and β-actin-F/R were designed using the Primer Premier 5 Design Program (Premier Biosoft International, Palo Alto, CA, USA) and were used to quantify the accumulation of transcripts, and to normalize the amounts of cDNA in samples, respectively. To ensure the specificity of PCR amplification or eliminate the interference of homologues in the other subgenomes, primers were selected with mismatched bases to specific homologues by mapping to genome sequences (EnsemblPlants, http://plants.ensembl.org/Triticum_aestivum/Info/Index). The reverse primer for miRNA was according to the instructions for the miRcute Plus qPCR kit (TIANGEN, Beijing). PCR was conducted in a 20-µL volume containing 10 µL 2 × SYBR Green PCR Master Mix (Takara, Dalian, China), 0.2 µM each primer and 2 µL template (6× diluted cDNA from leaf samples). The amplification program was as follows: 95 • C for 10 s; 40 cycles of 95 • C for 5 s, and 60 • C for 31 s. For each sample, reactions were carried out in triplicate and three non-template negative controls were included. Products were analyzed by melt curves obtained at the end of amplification, while the 2 −∆∆CT method was employed to quantify the relative gene expression. The correlation coefficients between co-expression genes were calculated with Pearson statistical method and t-test was used to test the statistical significance at the level of 0.05.

Identification of Transcription Factor Genes Adjacent to Differentially Expressed Long Non-Coding RNA in Wheat Responding to Pathogen Infection
Since lncRNAs play a regulatory role in the expression of nearby protein-coding genes and even gene clusters [15,22], we identified 461 functional genes close to 249 DE-lncRNAs, by alignment of the DE-lncRNA sequences with the reference sequences of Chinese Spring (IWGSC RefSeq v1.1) [29] (Supplemental Table S2). Among the functional genes, 27 transcription factors (TFs) were identified close to the interesting DE-lncRNAs, as listed in Table 1. These TFs could be classified into 15 types or families, including WRKY, NAC, MYB, C2H2, MADS, bHLH (Basic Helix-Loop-Helix), bZIP (basic leucine zipper), AP2/ERF (Activating Protein 2/ethylene responsive factor), CSD (the cold-shock domain), NF-X1 (nuclear transcription factor, X-box binding 1), B3 (plant-specific B3 superfamily), BES (BRI1-EMS suppressor), TUB (Tubby protein), GNAT (GCN5-related N-acetyltransferase), and mTERF (mitochondrial transcription termination factor). In addition, we similarly identified related resistance genes close to DE-lncRNAs (Table 2) (distance < 0.1 Mb) based on the annotation in IWGSC. The detailed information of locations was shown in the Supplemental Table S2. To further check the co-expression of DE-lncRNAs with these functional genes in response to pathogen stress, we used DESeq analysis, setting the threshold change at ≥2-fold and the false discovery rate (FDR) at 1.0%. The analysis identified 249 functional genes adjacent to 181 DE-lncRNAs as differentially expressed among pathogen-infected groups, compared with non-inoculated leaves as the control. This result suggested that about 75% of DE-lncRNAs might form DE-lncRNA-functional gene pairs. Note: The types of lncRNA were given in the symbols 'LincRNA, linncRNA, lpncRNA, and luncRNA' representing long intergenic ncRNAs, long intron ncRNAs, promoter lncRNAs, and untranslation region lncRNA. The star symbol mean that the functional genes were differential expressed in previously RNA-Seq profile.

Co-Expression of Long Non-Coding RNAs with Adjacent Functional Genes
By DESeq analysis of adjacent functional genes, we identified 6 DE-lncRNA-TF pairs using RNA-Seq data ( Table 2). Considering the hysteresis quality of the time-points in RNA-Seq data and the lack of comparable samples, we reassessed the relationship between the expression of DE-lncRNA and their nearby TFs. We randomly selected five predicted pairs of genes and tested the expression in resistant/susceptible NILs N9134R/S under powdery mildew pathogen stress. The gene expression pattern of lincRNA T10_unigene_BMK.12768 (hereafter abbreviated as T10.12768) was very similar to that of the adjacent gene TraesCS1B01G243100 (annotated as WRKY55-like), as shown in Figure 1.  Table S3, Value for significance at p = 0.05 and 8 df is 0.549). For TaNAC68L, T10.65297, TaNAC17L, T4.45663, the C2H2 gene and T16.1187, the maximum expression levels were detected at 168 hpi in N9134R. In the compatible line, co-expression was similarly detected in T10.12768 vs. WRKY55L, and T10.65297 vs. NAC68L, although the expression of the TFs showed stronger responses than their nearby lncRNAs at 24 and 120 hpi, respectively. However, the expression of NAC17L and TraesCS1B01G146800 (a C2H2 type TF) exhibited the opposite pattern in the susceptible line, especially at the early stage after Bgt inoculation. The inconformity of expression was also exhibited between T13.19448 and TraesCS4D01G172200.1 (homologous to WRKY64/70) in both compatible and incompatible lines. The expression of T13.19448 fluctuated only slightly in the resistant genotype. In contrast, the TF gene showed marked fluctuation: expression of WRKY64/70 was down-regulated at 6 hpi, then increased stepwise to a peak at 36 hpi, followed by a rapid decline back to a minimum at 120 hpi, but then increased again, reaching maximum accumulation at 240 hpi. Strikingly, the expression of lincRNA T13.19448 in the susceptible genotype increased progressively to a peak at 24 hpi, then decreased to a minimum at 72 hpi, followed by a slow increase to 240 hpi, but the expression of WRKY64/70 decreased rapidly at 6 hpi, followed by stable expression at the other time-points. This suggests that the lincRNA T13.19448 might negatively regulate the WRKY64/70 transcription factor at early stages. Note: The types of lncRNA were given in the symbols 'LincRNA, linncRNA, lpncRNA, and luncRNA' representing long intergenic ncRNAs, long intron ncRNAs, promoter lncRNAs, and untranslation region lncRNA. The star symbol mean that the functional genes were differential expressed in previously RNA-Seq profile.  Similarly, we compared the expression of DE-lncRNAs with that of several functional genes. The expression of linncRNA T10.79431 was nearly synchronized with the nearby gene TraesCS1B01G394900 (annotated as leaf rust 10 disease-resistance locus receptor-like protein kinase-like protein, R-RLK) in N9134R (Figure 2), while the expression of linncRNA T16.26398 and functional gene TraesCS5B01G565100 (a MAPK kinase) followed similar patterns in both resistant and susceptible lines. The correlation coefficient values reached to 0.998 and 0.968 in N9134R, respectively (Supplemental Table S3). The lncRNA T16.21355 and T13.23192 are long intergenic ncRNAs (lincRNAs). The former is located in the interval between TraesCS1B01G416400 (pre-mRNA-processing-splicing factor 8A, PRPF) and TraesCS1B01G416500 (ferredoxin-3 protein), while the latter is flanked by TraesCS1B01G384100 (TaDNAJ 10) and TraesCS1B01G384200 (unnamed protein). The physical distances from lincRNA T13.23192 to its two flanking functional genes are 31,632 and 50,821 bp, respectively. The distance between T16.21355 and TraesCS1B01G416400 is 11,232 bp, which is far less than the 232,840 bp distance from the lncRNA to TraesCS1B01G416500. These results showed that the expression patterns of lncRNAs were not totally coincident with functional genes, although similar expression patterns could be seen between them at several time-points ( Figure 2). For example, both T16.21355 and TraesCS1B01G416400 (splicing factor 8A) were upregulated at 36 and 120 hpi in the susceptible line, and at 168 hpi in the resistant line. T13.23192 and TaDNAJ 10 were induced at 168 hpi in the resistant N9134R, while T13.23192 and the unnamed protein gene TraesCS1B01G384200.1 were induced at 36 hpi in the susceptible N9134S. In addition, the unnamed protein gene expression peak at 36 hpi was accompanied by the other two peaks at 6 and 120 hpi, which showed an on-and-off pattern. These results showed that the lncRNAs are generally co-expressed with adjacent protein-coding genes, but not in all cases. It is noted that lncRNAs T10.65297, T4.45663, T10.79431, and T16.26398 are long intron ncRNAs, while T16.1187, T10.12768, T13.19448, T13.23192, and T16.21355 belong to long intergenic ncRNAs. The long intron ncRNAs (linncRNAs) seem to have closer co-expression relationships with their nearby functional genes. Taking the distance between lincRNAs with their adjacent functional genes into consideration, we noted that the distance from lincRNA T13.19448 to TraesCS4D01G172200 is 98,620 bp (r = 0.555), from T16.1187 to TraesCS1B01G146800 is 28,814 bp, and These results showed that the lncRNAs are generally co-expressed with adjacent protein-coding genes, but not in all cases. It is noted that lncRNAs T10.65297, T4.45663, T10.79431, and T16.26398 are long intron ncRNAs, while T16.1187, T10.12768, T13.19448, T13.23192, and T16.21355 belong to long intergenic ncRNAs. The long intron ncRNAs (linncRNAs) seem to have closer co-expression relationships with their nearby functional genes. Taking the distance between lincRNAs with their adjacent functional genes into consideration, we noted that the distance from lincRNA T13.19448 to TraesCS4D01G172200 is 98,620 bp (r = 0.555), from T16.1187 to TraesCS1B01G146800 is 28,814 bp, and from T10.12768 to TraesCS1B01G243100 is 6,625 bp. Thus, the strength of correlation of expression maybe increased as the distance decreased. This hints that two types of lncRNAs both play regulatory roles, in both positive and negative ways. However, the degree of regulation via lincRNA was influenced by the distance.

Co-Expression of Long Non-Coding RNAs and Allele-Specific Genes
Since lncRNAs, as miRNA targets or target mimics, competitively inhibited the interaction between miRNAs and target mRNAs to modulate gene expression, we identified five potential miRNA-targeted lncRNAs and seven mimic lncRNAs. Detailed information about the targets and target types of the lncRNAs are shown in Table 3. Intriguingly, all of the predicted target functional genes were identified as DE genes in the RNA-Seq profile. Among them, the linncRNA T16.13521 and phospholipid hydroperoxide glutathione peroxidase gene TraesCS6A02G246400 both are targeted and cleaved by tae-miR1137a. The linncRNAs T13.17661 and lincRNA T13.21716 are target mimics of miR339b and miR156d-3p, respectively, which were predicted to target TraesCS1B02G415800 (putative ubiquitin-conjugating enzyme E2) and TraesCS2D02G400500 (pentatricopeptide repeat-containing protein) for cleavage and translation inhibitor activity, respectively. The Q-PCR results showed that the gene expression level of lncRNA T13.17661 and TraesCS1B02G415800.1 (both targeted by miR399b) were low and stable after Bgt-inoculation, but the expression of miRNA399b was upregulated 4-7 fold compared to 0 hpi in N9134R resistant background (Figure 3). The expression of lncRNA T13.17661 and Ub-enzyme E2 were both induced at 12 hpi in N9134S susceptible background. However, their expression levels were lower than that of miR399b at subsequent time-points, especially at 36 and 48 hpi. As for miRNA1137a and its targets, the expression of lncRNA T16.13521 and functional gene TraesCS6A02G246400.1 were significantly repressed at all tested time-points, compared with 0 hpi, and accompanied by slight fluctuation, while the expression of miRNA1137a was slightly induced at 12 and 48 hpi. From Figure 3, the highest expression level of T16.13521 was observed at 24 hpi and was half that in the control. Interestingly, the trend of lncRNA T16.13521 was very similar to the functional gene TraesCS6A02G246400.1 (encoding glutathione peroxidase-like protein) at 0, 12, 24, 36, 48, and 72 hpi, while similar patterns could be seen comparing the expression of T13.17661 with TraesCS1B02G415800.1 in N9134R (Figure 3). Furthermore, the expression patterns were also similar between lncRNA T13.17661 and TraesCS1B02G415800.1 in N9134S, although the miR399b was induced at different time-points from the pattern in N9134R. The expression of TraesCS2D02G400500.1 (pentatricopeptide repeat-containing protein) steeply increased, with two peaks at 12 and 36 hpi in both N9134R and N9134S, while the lincRNA T13.21716 was upregulated at 24 hpi in N9134R. For miRNA156d as a translation inhibitor, the expression was generally stable after Bgt inoculation in both genotypes apart from slight down-regulation at 36 hpi in N9134R. Taken together, these results substantiated the view that lincRNA could competitively interact with functional genes via miRNA regulation.

Discussion
Plants are sessile and must continuously integrate both abiotic and biotic environmental signals for development and defense responses. Because plants lack circulating cells, they rely on systemic signals emanating from infection sites to trigger the innate immunity response [30]. In this process, thousands of genes have been involved, including noncoding RNAs (ncRNAs) encoded by specific genomic regions [31]. Large-scale sequencing analyses have revealed that most of the eukaryotic genome is transcribed to RNAs, including short and long ncRNAs [5,6,14], not just functional genes. In plants, analysis of over 200 Arabidopsis transcriptome data sets identified about 40,000 putative lncRNAs, including over 30,000 natural antisense transcripts (NATs) and over 6000 lincRNAs [10,13,15]. However, in-depth functional analysis of lncRNAs in wheat remains limited, although some preliminary reports have been given [7,14,32]. Here, we identified 461 neighboring genes with 249 DE-lncRNA in wheat after fungal infection, and then investigated the co-expression relationship

Discussion
Plants are sessile and must continuously integrate both abiotic and biotic environmental signals for development and defense responses. Because plants lack circulating cells, they rely on systemic signals emanating from infection sites to trigger the innate immunity response [30]. In this process, thousands of genes have been involved, including noncoding RNAs (ncRNAs) encoded by specific genomic regions [31]. Large-scale sequencing analyses have revealed that most of the eukaryotic genome is transcribed to RNAs, including short and long ncRNAs [5,6,14], not just functional genes. In plants, analysis of over 200 Arabidopsis transcriptome data sets identified about 40,000 putative lncRNAs, including over 30,000 natural antisense transcripts (NATs) and over 6000 lincRNAs [10,13,15]. However, in-depth functional analysis of lncRNAs in wheat remains limited, although some preliminary reports have been given [7,14,32]. Here, we identified 461 neighboring genes with 249 DE-lncRNA in wheat after fungal infection, and then investigated the co-expression relationship between lncRNA and the adjacent functional genes. Regretfully, because of missing the direction of cDNA in the second-generation RNA-Seq, here we were not able to detect possible co-regulation relationships between lncNATs and their target genes. There is no doubt that further lncRNA regulators will be detected using strand-specific RNA-Seq [33]. In any case, as far as we know, this is the first genome-wide study on the possible role of lncRNAs in regulating functional genes. Our results provide a powerful foundation for future functional and molecular research on wheat-fungus interactions.
LncRNAs, like functional genes, have key functions in transcriptional, post-transcriptional, and epigenetic gene regulation [2]. In transcriptional level, transcription factors (TFs) are important regulators of gene expression in plants when responding to abiotic and biotic stress [34,35]. The NAC, WRKY, AP2/ERF, and C2H2 TFs have all been reported as being involved in plant responses to pathogen attack [36]. In this study, we analyzed the co-expression of lncRNAs with NAC17L, NAC68L, WRKY55L, C2H2, and WRKY64/70. The results supported the idea that lncRNAs have a co-regulation relationship with neighboring functional genes, although some exhibited opposite and/or positive expression patterns in different genetic backgrounds, and others showed no relation at all. WRKY transcription factors play important and well substantiated roles in plant immunity responses at both transcriptional and post-translational regulation levels [37,38]. Several NAC and C2H2 have been identified with roles in antioxidant defense mechanisms in plants [39], in particular a natural allele of a C2H2 transcription factor in rice that was shown to confer broad-spectrum rice blast resistance [40]. Here, our results demonstrated that lncRNAs maybe play the role of regulating the expression of TFs, especially in incompatible interaction. This means that lncRNAs could indirectly regulate the downstream functional genes via TFs. Taking co-regulation of lncRNAs with allele-specific functional genes directly together, our results will be helpful for improved understanding of the regulatory mechanisms of TFs in plant immune responses to disease. Comparing linncRNA with lincRNA, the linncRNA seemed to have closer co-expression relationships with their nearby functional genes, while the correlation of lincRNA expression with their respective genes increased as the intervening distance decreased. These results suggest that lncRNAs might have regulatory functions for neighboring functional genes, including transcription factors, in both positive and negative ways.
Conversely, MicroRNAs are endogenous short ncRNAs (21-24 nucleotides) that play important regulatory roles by repressing gene translation or degrading target mRNAs at the post-transcriptional levels [6,41]. Some lncRNAs have been predicted to be targets of miRNA [9,14,32,42]. To further clarify the relationship between lncRNA and miRNA, here we analyzed the co-expression of miRNAs, lncRNAs and their targeted functional genes, and substantiated the view that lncRNAs, including lncRNA targets and target mimics, are able to interact competitively with functional genes via miRNA regulation at the post-transcriptional levels [20]. However, further work is required to determine whether lncRNAs could be immune resistance markers as functional genes. Taking the regulation of miRNA with lncRNA together, the regulatory network of functional gene expression will be probably complex and flexible. Although it is not yet clear how lncRNAs are involved in regulation and with which effectors they interact, the knowledge and resources gathered here will provide useful insights into the mechanisms that regulate defense pathways against fungi, and these results should facilitate future investigating into epigenetic resistance in wheat, as well as understanding the plant-pathogen interactions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/6/896/s1, Supplemental Table S1. PCR primers used for Q-PCR amplification of lincRNA and functional genes. Supplemental Table S2. List of the identified adjacent functional genes nearby DE-lncRNAs regulated by fungi. Supplemental  Table S3. List of the detailed relative expression values used to infer the correlation.