Genome-wide identi cation and characterization of long non-coding RNAs conferring resistance to Colletotrichum gloeosporioides in walnut (Juglans regia)

Background Walnut anthracnose caused by Colletotrichum gloeosporioides (Penz.) Penz. and Sacc. is an important walnut production problem in China. Although the long non-coding RNAs (lncRNAs) are important for plant disease resistance, the molecular mechanisms underlying resistance to C. gloeosporioides in walnut remain poorly understood. Results The anthracnose-resistant F26 fruits from the B26 clone and the anthracnose-susceptible F423 fruits from the 4–23 clone of walnut were used as the test materials. Specifically, we performed a comparative transcriptome analysis of F26 and F423 fruit bracts to identify differentially expressed LncRNAs (DELs) at five time-points (tissues at 0 hpi, pathological tissues at 24 hpi, 48 hpi, 72 hpi, and distal uninoculated tissues at 120 hpi). Compared with F423, a total of 14,525 DELs were identified, including 10,645 upregulated lncRNAs and 3846 downregulated lncRNAs in F26. The number of upregulated lncRNAs in F26 compared to in F423 was significantly higher at the early stages of C. gloeosporioides infection. A total of 5 modules related to disease resistance were screened by WGCNA and the target genes of lncRNAs were obtained. Bioinformatic analysis showed that the target genes of upregulated lncRNAs were enriched in immune-related processes during the infection of C. gloeosporioides, such as activation of innate immune response, defense response to bacterium, incompatible interaction and immune system process, and enriched in plant hormone signal transduction, phenylpropanoid biosynthesis and other pathways. And 124 known target genes for 96 hub lncRNAs were predicted, including 10 known resistance genes. The expression of 5 lncRNAs and 5 target genes was confirmed by qPCR, which was consistent with the RNA-seq data. Conclusions The results of this study provide the basis for future functional characterizations of lncRNAs regarding the C. gloeosporioides resistance of walnut fruit bracts. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-020-07310-6.


Results
The anthracnose-resistant F26 fruits from the B26 clone and the anthracnose-susceptible F423 fruits from the 4-23 clone of walnut were used as the test materials. Speci cally, we performed a comparative transcriptome analysis of F26 fruit bracts and F423 to identify differentially expressed LncRNAs (DELs) at ve time-points (tissues at 0 hpi, pathological tissues at 24 hpi, 48 hpi, 72 hpi, and distal uninoculated tissues at 120 hpi). Compared with F423, a total of 14525 DELs were identi ed, including 10645 upregulated lncRNAs and 3846 downregulated lncRNAs in F26. The number of upregulated lncRNAs in F26 compared to in F423 was signi cantly higher at the early stages of C. gloeosporioides infection. A total of 5 modules related to disease resistance were screened by WGCNA and the target genes of lncRNAs were obtained. Bioinformatic analysis showed that the target genes of upregulated lncRNAs were enriched in immune-related processes during the infection of C. gloeosporioides, such as activation of innate immune response, defense response to bacterium, incompatible interaction and immune system process, and enriched in plant hormone signal transduction, phenylpropanoid biosynthesis and other pathways. And 124 known target genes for 96 hub lncRNAs were predicted, including 10 known resistance genes. The expression of 5 lncRNAs and 5 target genes was con rmed by qPCR, which was consistent with the RNA-seq data.

Conclusions
The results of this study provide the basis for future functional characterizations of lncRNAs regarding the C. gloeosporioides resistance of walnut fruit bracts. Background Walnut (Juglans regia L.) is a diploid tree species (2n = 32), with approximately 667 Mb per 1C genome [1]. It is an ecologically important 'woody oil' tree species worldwide [2], and its kernel is a rich source of nutrients with health bene ts for humans [3]. Regarding the commercial production of walnut in China, Colletotrichum gloeosporioides (Penz.) Penz. and Sacc infects the walnut fruit bract and causes anthracnose, leading to premature fruit drop and yield losses. The C. gloeosporioides lifestyle transitions associated with the infection of the host include the following three stages: attachment, biotrophy, and necrotrophy [4]. Speci cally, the formation of adherent cells is critical for fungal development during the C. gloeosporioides infection [5]. In a previous study, LAC2 was revealed to contribute to the formation of adherent cells to enhance the pathogenicity of C. gloeosporioides [6]. However, it is unclear how walnuts recognize and resist infections by C. gloeosporioides, and the regulatory network of hub and peripheral genes underlying the resistance of walnuts to C. gloeosporioides remains uncharacterized. Therefore, elucidating the molecular basis of this resistance mechanism is imperative for the breeding of walnut resistant to C. gloeosporioides [7][8][9].
Long non-coding RNA (lncRNA) is a type of RNA comprising 200-1,000,000nt and structural characteristics similar to those of mRNA, but it does not encode a protein [10]. The lncRNAs were initially considered to be the transcription 'noise' of protein-coding genes, and were often ignored in transcriptome analyses [11]. However, the continuous development of sequencing technologies and transcriptome analyses has revealed that many lncRNAs in Arabidopsis thaliana [12], Triticum aestivum [13], Zea mays [14], and other plant species are related to stress responses, morphological development, and fruit maturation. For example, a heat-responsive lncRNA (TCONS_00048391) is an eTM for bra-miR164a and may be a competing endogenous RNA (ceRNA) for the target gene NAC1 (Bra030820), with effects on bra-miR164a expression in Chinese cabbage (Brassica rapa ssp. chinensis) [15]. Qin et al con rmed that the DROUGHT INDUCED lncRNA regulates plant responses to abiotic stress by modulating the expression of a series of stress-responsive genes [16]. In A. thaliana, two lncRNAs, COOLAIR and COLDAIR, are associated with FLOWERING LOCUS C and are crucial role for vernalization [17,18].
Many recent studies have proven that lncRNAs are important for plant-pathogen interactions. A role for nine hub lncRNAs and 12 target genes in the resistance of Paulownia tomentosa to witches'broom was uncovered via a high-throughput sequencing experiment, and their functions were analyzed with an RNA-lncRNA co-expression network model [19]. In tomato (Solanum lycopersicum), the lncRNA16397-GRX21 regulatory network reportedly decreases the reactive oxygen species content and cell membrane damage to enhance the resistance to P. infestans [20]. Moreover, the involvement of the WRKY1-lncRNA 33732-RBOH module in regulating H 2 O 2 accumulation and resistance to P. infestans was determined based on a comparative transcriptome analysis [21]. In cotton (Gossypium spp.), a functional analysis demonstrated that a lack of two hub lncRNAs, GhlncNAT-ANX2 and GhlncNAT-RLP7, enhances seedling resistance to Verticillium dahliae and Botrytis cinerea, possibly because of the associated upregulated expression of LOX1 and LOX2 [22]. In wheat (Triticum aestivum L.), lncRNAs have a tissue-dependent expression pattern that can respond to powdery mildew infections and heat stress [23]. Additionally, four kinds of lncRNAs have important effects on Puccinia striiformis infections [24]. However, there are no reports regarding the role of lncRNAs in the walnut fruit resistance to anthracnose.
In this study, Illumina HiSeq 4000 sequencing was used to analyze the disease-resistant (F26) and susceptible (F423) fruit bracts at different C. gloeosporioides infection stages. The number and characteristics of lncRNAs were analyzed. Additionally, the hub lncRNAs related to disease resistance were screened and functionally analyzed to predict the role of lncRNAs in walnut fruit bract resistance to anthracnose. To the best of our knowledge, this is the rst report on walnut lncRNAs and their biological functions related to fruit bract resistance to C. gloeosporioides. Our data may be a useful resource for clarifying the regulatory functions of lncRNAs in uencing walnut fruit resistance to C. gloeosporioides.

Results
Whole genome identi cation of lncRNAs expressed in walnut fruit bracts To identify lncRNAs expressed in walnut fruits in response to C. gloeosporioides, we constructed 20 cDNA libraries from the anthracnose-resistant and the anthracnose-susceptible walnut fruits at the following ve infection stages: tissue at 0 hpi (hours post inoculation), infected tissue at 24, 48, and 72 hpi, and distal uninoculated tissue at 120 hpi (Additional le 1: Table S1). The libraries were sequenced with an Illumina HiSeq 4000 platform. A total of 265.4 Gb clean data were obtained, with an average of 13.27 Gb per library. Approximately 69.7% of the clean reads in all libraries were mapped to the walnut reference genome (Additional le 2: Table S2). The aligned transcripts were assembled, combined, and screened with the FEELnc software to obtain 22,336 lncRNAs (length ≥ 200 bp, ORF coverage < 50%, and potential coding score < 0.5), including 18,403 unknown lncRNAs (23.97%) and 3,933 known lncRNAs (5.12%) ( Fig   1).

Characterization of walnut fruit bract lncRNAs
A total of 58,369 mRNAs and 22,336 lncRNAs were obtained for the walnut fruit bracts (all samples combined) (Additional le 3: Table S3 S4). The lncRNAs were characterized according to their locations relative to the partner RNA. A total of 40,429 (67.57%) lncRNAs were located in intergenic regions (i.e., only 32.43% genic lncRNAs). Additionally, 19,767 (48.89%) and 7,302 (37.63%) of the intergenic lncRNAs and genic lncRNAs were located in the antisense strand, respectively (Fig 2a) (Additional le 5: Table S5). Most lncRNAs contained two or three exons, which differentiated them from mRNAs (Fig 2c). Moreover, there was considerable diversity in the distribution of mRNA and lncRNA lengths (Fig 2b). Furthermore, the expression level of most lncRNAs was signi cantly lower than that of mRNAs (Fig 2d).

Differentially expressed lncRNAs at various infection stages
The lncRNAs that were differentially expressed between the disease-susceptible F423 fruits and the disease-resistant F26 fruits at different C. gloeosporioides infection stages were analyzed. Compared with F423, a total of 14,525 DELs were identi ed, including 10,645 up-regulated lncRNAs and 3,846 downregulated lncRNAs in F26. The number of upregulated and downregulated lncRNAs in the various comparisons were respectively as follows: 7,668 and 1,386 in the F26_0hpi vs F423_0hpi comparison; 6,910 and 1,165 in the F26_24hpi vs F423_24hpi comparison; 1,721 and 1,593 in the F26_48hpi vs F423_48hpi comparison; 898 and 1,133 in the F26_72 hpi vs F423_72 hpi comparison; and 4,711 and 550 in the F26_120 hpi vs F423_120 hpi comparison (Fig 3a, b) Table S7). These results revealed the similarities in the expression of lncRNAs and mRNAs. And the number of upregulated lncRNAs and mRNAs in F26 compared to in F423 was signi cantly higher at the early stages of C. gloeosporioides infection.
Identi cation of co-expressed lncRNA modules To identify the hub lncRNAs and predict their potential target genes in trans-regulatory relationships, a weighted gene co-expression network analysis (WGCNA) was used to generate a correlation matrix of the expression levels of 10,645 upregulated lncRNAs and 15,247 upregulated mRNAs. A total of 19 expression modules were screened (Fig 4a) (Additional le 8: Table_S8). The relationships between modules and the resistance traits of the walnut fruit bracts were analyzed and four signi cantly correlated modules (r ≥ 0.8) were identi ed. The MEviolet module was correlated with F26_0hpi (r = 0.95, p = 9e−11), which contains 406 lncRNAs and 1350 mRNAs. The MElightyellow module was correlated with F26_24hpi ( r = 0.86, p = 1e−06 ), which contains 165 lncRNAs and 892 mRNAs. The MEbrown2 module was correlated with F26_48hpi (r = 0.82, p = 8e−0.86), which contains 128 lncRNAs and 224 mRNAs. The MEwhite module was correlated with F26_72hpi (r = 0.81, p = 1e−05), which contains 111 lncRNAs and 378 mRNAs (Fig 4c). Regarding F26_120 hpi, the rand p value for the MEorange module was 0.73 and 3e−0.4, respectively. The highest r value (0.77) for F423 was calculated for the MEdarkseagreen module and F423_48hpi (Fig 4b). And the MEorange module contains 76 lncRNAs and 227 mRNAs (Fig 4c). These results suggested that lncRNAs are closely related to the disease resistance of walnut fruit bracts.

Network analysis of hub lncRNAs
The hub lncRNAs are important for regulating the whole network. Therefore, we screened the 96 hub lncRNAs and 124 known target genes according to their weight value and connectivity in ve modules (Additional le 11: Table_S11). In the MEviolet module, the 25 known target genes for 15 hub lncRNAs were found to be involved in multiple functions (Fig 10a), such as probable galacturonosyl transferase 10 and ultraviolet-B receptor UVR8-like. In addition, target genes encoding receptor-like serine / threonineprotein kinase NCRK (XM_018958556.1) and eukaryotic translation initiation factor 5A-2-like (XM_018994862.1) are known resistance genes (Fig 10b). In the MElightyellow module,16 hub lncRNAs were generated and their 22 known target genes were involved in many functions (Fig 11a). And the target genes encoding G-type lectin S-receptor-like serine / threonine-protein kinase LECRK1 (XM_018950446.1), probably inactive leucine-rich repeat receptor-like protein kinase At2g25790 (XM_018989953.1) and TMV resistance protein N-like (XM_018961957.1) were known resistance genes (Fig 11b). In the MEbrown2 module, 24 hub lncRNAs and their 15 known target genes were generated (Fig  12a), the target gene encoding probable LRR receptor-like serine / threonine-protein kinase At3g47570 (XM_018962714.1) was konwn resistance gene (Fig 12b). In the MEwhite module, 23 hub lncRNAs were generated and their 38 known target genes were involved in many functions (Fig 13a). The target genes encoding putative disease resistance protein At1g50180 (XM_018965430.1) , probable LRR receptor-like serine/threonine-protein kinase At1g63430 (XM_018973294.1) and L-type lectin-domain containing receptor kinase IV.2-like (XM_018954279.1) were konwn resistance genes (Fig 13b). In the MEorange module, 18 hub lncRNAs were generated and their 24 known target genes were involved in many functions (Fig 14a). And the target gene encoding the inactive LRR receptor-like serine / threonine-protein kinase BIR2 (XM_018967526.1) was konwn resistance gene (Fig 14b). All disease resistance genes in walnut are listed in Additional le 12: Table_S12. These results suggested that lncRNAs may participate in the resistance of walnut bracts to C. gloeosporioides by acting on their target genes.

Discussion
In previous studies, lncRNAs were identi ed and analyzed in various biological processes important for seed development [35], photomorphogenesis [36], fruit development [37,38], and biotic and abiotic stress responses [16,39]. Additionally, there has been substantial research on the role of lncRNAs in plantpathogen interactions. In A. thaliana, lncRNAs reportedly enhance the resistance to Pseudomonas syringae pv. tomato DC3000 by promoting PR1 expression [40]. In tomato, lncRNA23468 functions as a ceRNA that modulates NBS-LRR gene expression by mimicking the target of miR482b, thereby increasing the resistance to P. infestans [41].Walnut anthracnose has been responsible for the premature fruit drop and yield losses that have adversely affected walnut production in China [7]. In this study, we investigated the role of lncRNAs in the resistance of walnut fruit bracts to anthracnose based on sequence analyses. Walnut anthracnose is caused by C. gloeosporioides, which completes its infection process as a hemibiotroph [4,42]. First, conidia germinate to generate appressoria, which produce invasion pegs that initiate the infection into susceptible plants. The primary mycelium produced in plant cells exists as a biotroph, after which the secondary mycelium produced in the infected site switches to necrotrophic growth [43,44]. We previously determined that the C. gloeosporioides life cycle in walnut tissue involves attachment at 24hpi, biotrophy at 48hpi, and necrotrophy at 72hpi (data unpublished). In this study, RNAseq was performed to build the lncRNA and mRNA pro les of the walnut fruit bract tissue at 0 hpi, infected tissue at 24, 48, and 72 hpi, and distal uninoculated tissue at 120 hpi. A total of 58,369 mRNAs and 22,336 lncRNAs were identi ed, including 3,933 known lncRNAs and 18,403 unknown lncRNAs. Consistent with the results of similar studies on other organisms, the identi ed putative lncRNA had fewer exons, shorter transcripts, and lower expression levels than protein-coding genes [45,46].
The development of bioinformatic analysis technology has enabled researchers to reveal that lncRNA functions and characteristics are far more complex than previously thought [10]. A recent comparative transcriptome analysis between wild-type and WRKY1-overexpressing tomato plants revealed 199 lncRNAs (DELs) and indicated that many of the lncRNA target genes that are likely affected by WRKY1 and associated with the resistance of tomato to P. infestans are involved in the response to biotic stimulus (GO:0009607) and plant-pathogen interaction (KO4626) [20]. In another recent study, 4,594 putative lncRNAs were identi ed in comprehensive dynamic lncRNA expression networks under heat stress conditions. Co-expression networks revealing the interactions among the differentially expressed lncRNAs, mRNAs, and microRNAs indicated that several phytohormone pathways are associated with heat tolerance, including salicylic acid and brassinosteroid pathways [15]. In the current study, we obtained 10,645 upregulated lncRNAs and 15,247 upregulated mRNAs among the ve comparisons (F26_0hpi vs F423_0hpi, F26_24 hpi vs F423_24 hpi, F26_48 hpi vs F423_48 hpi, F26_72 hpi vs F423_72 hpi, and F26_120 hpi vs F423_120 hpi). The number of up-regulated lncRNAs and mRNAs in the F26 vs F423 was signi cantly higher at the early stages of C. gloeosporioides infection.
The functions of lncRNAs cannot currently be inferred from their sequence or structure, but lncRNAs can function in trans mode to target gene loci distant from where the lncRNAs are transcribed [47]. In F26, a total of 5 modules related to disease resistance were obtained by WGCNA during the infection of C. gloeosporioides. Many target genes of lncRNAs in these modules are enriched in plant immune related items and pathways, such as "activation of innate immune response", "activation of immune response" in MElightyellow module, "defense response to bacterium, incompatible interaction" in MEbrown2 module, "defense response "and "immune system process" in MEwhite module. These results suggest that these genes may play important roles in the process of resistance to C. gloeosporioides of walnut fruit bracts. Phytohormones are known to be important in the regulation of defense responses in plants [48][49][50]. Plants can exhibit systemic acquired resistance through the salicylic acid (SA) / jasmonic acid (JA)-mediated signaling network [51,54]. In our study, a total of 32 genes were identi ed in the signi cantly enriched KEGG pathway "Plant hormone signal transduction". Meanwhile, there are 3 and 5 genes enriched in "jasmonic acid mediated signaling pathway"and "response to jasmonic acid" respectively. We also showed that some genes were enriched in "auxin-activated signaling pathway" and "cellular response to auxin stimulus" at 24 hpi . Therefore, auxin may play a role in the resistance of walnut bracts to C. gloeosporioides. In addition, our result showed that the phenylpropanoid biosynthesis was one of the most signi cantly enriched pathways in the process of resistance to C. gloeosporioides of walnut fruit bracts. In this pathway, phenylalanine ammonium lyase (PAL) is the key regulatory enzyme in altering the biosynthesis and accumulation of avonoids and lignin [55]. Lignin plays a structural role in the secondary cell walls formation [56], and avonoids mediate plants against UV radiation and act as a visual signal for attracting pollinators [57,58]. In Caragana korshinskii, C. korshinskii adjusts its phenylpropanoid biosynthesis process to water-de cit environments and activates PAL by drought stress [59].
During long-term evolutionary interactions with plants, several pathogens successfully cause effectortriggered susceptibility response (ETS) by producing a number of effectors. Simultaneously, plants have evolved R genes that recognize these effectors and function through highly speci c interactions between effectors and their corresponding nucleotide-bindingsite and leucine-rich repeat (NB-LRR) class receptors [60]. In tomato, lncRNA23468 reportedly increases the expression of the NBS-LRR target genes (encoding R proteins), resulting in enhanced resistance to P.infestans [41]. In the current study, we detected 10 R genes among the target genes of 96 hub lncRNAs. During the infection of C. gloeosporioides on the walnut fruit bracts, the results of RNA-seq showed that the expression of R genes (XM_018950446. ). These ndings imply that lncRNAs may help mediate the disease resistance of walnut fruit bracts through the target R genes. The speci c interaction between lncRNAs and R gene needs further veri cation. The expression levels of ve hub lncRNAs (MSTRG13585, MSTRG11713, MSTRG152205, MSTRG112028, and MSTRG62751) and their target genes were further con rmed by qPCR, the results of which were consistent with the RNA-seq data. The data presented here provides researchers with the biological basis for future investigations of the mechanism underlying the disease resistance of walnut fruit bracts.

Conclusions
In this study we generated the expression pro le of lncRNA in anthracnose-resistant F26 and anthracnose-susceptible F423 at ve times. Compared with F423, a total of 14525 DELs were identi ed, including 10645 upregulated lncRNAs and 3846 downregulated lncRNAs in F26. Bioinformatic analysis showed that the target genes of upregulated lncRNAs were enriched in immune-related processes , plant hormone signal transduction, phenylpropanoid biosynthesis and other pathways during the infection of C. gloeosporioides. And hub lncRNAs with high connectivity to disease resistant genes were predicted. These results contribute to our understanding of the potential mechanism by which lncRNAs involved in C. gloeosporioides resistance and will facilitate the functional veri cation of the lncRNA in the future.

Plant materials and fungal isolates
The scions of walnut seedling trees B26 was provided by walnut specialized farmers' cooperative of Dongliugang village, Baishi Town, Wenshang County, Shandong Province, China (35°46′56.2″N, 116°40′30.8″E). The 4-23 walnut tree was from F1 progeny of an intraspeci c cross between walnut cultivar 'Yuan Lin' (susceptible to anthracnose)× 'Qing Lin' (resistant to anthracnose)which was carried out by ourselves in 2002. The plant materials were conserved by patch budding onto walnut seedling rootstock at the Forestry Experimental Station of Shandong Agricultural University, Tai'an, Shandong Province, China (36°10′ 19.2″N, 117°09′ 1.3″E) in late May 2009. In 2015-2017, we evaluated the anthracnose resistance of each plant for three consecutive years followed by previously described [7,8], and it was found that B26 clone was high resistance to anthracnose in fruit bract, and the 4-23 clone was high susceptible to anthracnose in fruit bract. The fruits of B26 clone (i.e., F26) and 4-23 clone (i.e., F423) were used as experimental materials. The voucher specimen of F26 and F423 had been deposited to our lab but not to any publicly available herbarium.We didn't use wild plants in this study and according to national and local legislation, no speci c permission was required to collect these plants. C.gloeosporioidesm9 isolates (GenBank ID: GU597322) used in this study were maintained by our group.

Fungal pathogen inoculation of walnut fruits
Colletotrichum gloeosporioides was cultivated on potato dextrose agar medium for 5-7 days at 28℃.To prepare conidial suspensions, the colonies were washed with sterile distilled water containing 0.05% (v/v) Tween 80, passed through a lter (40-100 μm pores), quanti ed with a hemocytometer, and diluted with sterile distilled water to 105-106 conidia/ml [0.001% (v/v) Tween 80 nal concentration]. Healthy fruits from the east-, south-, and west-facing parts of each tree were collected in mid-June and disinfected with 0.6% sodium hypochlorite and rinsed with sterile water. The punch inoculation of the detached walnut fruits was completed as previously described [8]. Based on anatomical changes to the infected walnut fruit bract, samples of the inoculation site were collected at 0, 24, 48, and 72 hpi. Additionally, distal uninoculated tissue was collected at 120 hpi. Take two independent samples as biological replicates at each infection time (Additional le 1: Table S1). All samples were ash-frozen in liquid nitrogen and stored at −80°C until analyzed.
RNA extraction, library construction, and sequencing Total RNA was extracted from F423 and F26 samples with the Thermo Gene JET Plant RNA Puri cation Mini Kit (Thermo Fisher Scienti c Inc., USA). The purity and concentration of the extracted RNA were determined with the NanoDrop2000 spectrophotometer (Thermo Fisher Scienti c Inc.) (OD 260/280 ≥ 1.8, OD 260/230 ≥ 1.5, and concentration > 40 ng/µl). The RNA integrity was assessed by agarose gel electrophoresis. Ribosomal RNA was removed with the Ribo-Zero™ Magnetic Kit (Epicentre) and the remaining RNA (polyA+ and polyA−) was recovered. The RNA was randomly fragmented to approximately 200-bp sequences in Fragmentation Buffer (Thermo Fisher Scienti c Inc.) and then used as the template to synthesize rst-strand cDNA with random hexamers. The second cDNA strand was synthesized with dNTPs, RNaseH, and DNA polymerase I. The overhanging ends were lled in with T4 DNA polymerase and Klenow DNA polymerase to generate blunt ends, after which the A base was added to the 3′ end and the fragment was ligated to a linker. The AMPureXP beads were used for selecting fragments. The second cDNA strand containing U was degraded with the USER enzyme, after which a sequencing library was obtained by PCR ampli cation. A total of 20 sequencing libraries were constructed. The Qubit 2.0 DNA Broad Range Assay (Invitrogen, USA) was used for a preliminary quanti cation. The sequencing library inserts were detected with the Agilent 2100 Bioanalyzer. Finally, the effective library concentrations (> 2 nM) were accurately quanti ed by qPCR. Paired-end sequencing (2 × 150 bp) was performed with an Illumina HiSeq 4000 platform.

Identi cation of lncRNAs
The FEELncprogram (version 0.1.1) was used to identify lncRNAs in the transcriptome assemblies. First, the FEELnc lter was used to remove short transcripts (default 200nt) and assess single-exon transcripts [28]. The FEELnc codpot predictors were used to calculate a coding potential score. The assembled sequences were used for reconstructing the transcriptome. Finally, RNAs longer than 200bp and derived from ≥ 2 exons, with an ORF coverage < 50% and a potential coding score < 0.5 were designated as lncRNAs [29].

Classi cation of lncRNAs
The lncRNAs were analyzed regarding their corresponding positions in the reference genome and the positional relationships between lncRNAs and partner RNAs based on 10,000-100,000 fragments. The lncRNAs were then divided into genic lncRNAs (overlapping partner RNAs) and intergenic lncRNAs (lincRNAs). The genic lncRNAs were further divided as overlapping, containing, or nested subtypes.
Intergenic lncRNAs were divided as divergent, convergent, and same strand subtypes.

Analysis of differential expression patterns
Genes differentially expressed between the disease-resistant and susceptible fruits at ve infection stages were analyzed with DESeq2 (version 1.22.1) [30]. After assessing the signi cance of any differences, the genes with a p value ≤ 0.05 and a |log 2 foldchange| ≥ 1 were designated as differentially expressed genes. Prediction of lncRNA functions based on co-expression and genomic co-localization Co-expression modules were generated with the WGCNA package (version 1.67) as previously described [31] (http://lab.genetics.ucla.edu/horvath/Coexpression Network/). The lncRNAs and mRNAs that were not detected in at least one infection stage were not considered. In this analysis, the soft thresholding power was set to 12, after which the adjacency function was used to construct the adjacency matrix. A topological overlap measure map was constructed based on the adjacency matrix to calculate the similarity matrix of the lncRNA and mRNA expression between different nodes. The lncRNAs and mRNAs were hierarchically clustered based on the algorithm. To generate a number of clusters, modules were de ned after eliminating or combining branches. The co-expression module dynamic shear tree parameters were determined as described by Gerttula [32]. The minimum number of genes was set to 30, the split sensitivity (deep Split) was set to 2, and the other settings were software default parameters. The module was related to the trait, and the correlation matrix between the module and the trait was calculated. The module with the highest correlation coe cient and the smallest p value was designated as the module most relevant to the trait. In this study, a signi cantly correlated module was identi ed based on a correlation coe cient (r) ≥ 0.8 [29] and p < 0.05. The co-expression networks of lncRNAs and hub lncRNAs in highly correlated modules were generated with the Cytoscape software (version 3.7.1) [33].

Functional enrichment analysis
The genes targeted by lncRNAs were functionally annotated based on the GO and KEGG pathway (http://www.genome.jp/kegg/) databases. The KOBAS program (version 2.0) was used to determine the signi cantly enriched KEGG pathways among the target genes [34]. The data sets are included within the article and its Additional les. The raw sequencing data were deposited in NCBI Sequence Read Archive under the accession number GSE147083 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE147083 ).

Competing interests
The authors declare that they have no con icts of interest.

Funding
This study was supported by a Grant from the National Natural Science Foundation of China (No. 31870667 and No. 31470680). The funders had no role in the design of the study, collection, data analysis, data interpretation and manuscript writing.
Authors' contributions KQY conceived the idea and revised the manuscript. SF and HF conducted meta-analysis, drew gures and drafted the manuscript. XL collected the experimental materials and YD helped in drawing gures and drafting the manuscript. All authors listed have made direct and substantial efforts for improving the manuscript and approved the nal version. Validation of selected lncRNAs and mRNAs in a quantitative PCR assay. Blue and red represent the F423 and F26 samples, respectively. Expression data were normalized against the data for the18S rRNA housekeeping gene and are presented as themean ± standard error; *p < 0.05, **p <0.01.