Genome-wide identification of long non-coding (lncRNA) in Nilaparvata lugens’s adaptability to resistant rice

Plants and insects

For this study, two virulent BPH populations were full sib-mated for at least 40 generations and fed with TN1 and YHY15 plants in a climate chamber under a 16-h light/8-h dark cycle at 25 ± 1 °C and 70% relative humidity at the Hubei Academy of Agricultural Sciences, Wuhan, China. Twenty insects of the two virulent BPH populations (hereinafter referred to as TN1 and YHY15) were obtained, respectively. Whole insects of the two samples were immediately frozen in liquid nitrogen and stored at −70 °C until RNA isolation.

Construction of lncRNA sequencing library and RNA sequencing

Total RNA was extracted from brown planthopper using a Trizol kit (Invitrogen, Carlsbad, CA, USA) and further purified with the RNeasy kit (Qiagen, Hilden, Germany). The purity and quality of RNA were determined using Nanodrop instrument. Total RNA (1-μg) from TN1 or YHY15 samples (three biological replicates for each treatment: a total of six samples) was taken for library construction. The ribosomal RNA was removed to retain all coding RNA and ncRNA. The first-strand of cDNA was synthesized using random primers. Then, the second-strand cDNA was synthesized using dNTPs (dUTP instead of dTTP), RNase H and DNA polymerase I and was purified by QIAquick PCR Purification Kit. After end repair, poly(A) and sequencing adaptor addition, the cDNA was digested by UNG (Uracil-N-Glycosylase) enzyme. The fragments with expected size were selected by agarose gel electrophoresis. Finally, the constructed libraries were sequenced with Illumina Hiseq 4000.

lncRNA prediction and new transcription analysis

The raw data were processed using filtering adaptors, containing more than 50% of low quality (Q-value ≤ 20) bases, and trimming the reads whose number of N bases accounted for more than 10% of the total by fastp (version 0.18.0) (Chen et al., 2018). The reference BPH genome and the annotation files were downloaded from the National Center for Biotechnology Information database with the accession NO. GCF_014356525.1 or InsectBase (http://insect-genome.com/planthoppers/) (Ma et al., 2021). The rRNA mapped reads were removed after short reads mapping to the ribosome RNA (rRNA) database by Bowtie2 (version 2.2.8) (Langmead & Salzberg, 2012). The remaining reads were further mapped to the reference genome using HISAT2 (version 2.1.0) with “-rna-strandness RF” and other parameters set as a default (Kim, Langmead & Salzberg, 2015). The transcript reconstruction was performed using software Stringtie (version 1.3.4) and HISAT2 (Kim, Langmead & Salzberg, 2015; Pertea et al., 2016). Two softwares, CNCI (version 2) and CPC (version 0.9-r2) (http://cpc.cbi.pku.edu.cn/) were used to assess the protein-coding potential of novel transcripts by default parameters (Kong et al., 2007; Sun et al., 2013). The intersection of both non-protein-coding possible results were chosen as long non-coding RNAs.

Differential expression genes screening and function enrichment analysis

Software StringTie quantified transcript abundances in a reference-based approach. A FPKM (fragment per kilobase of transcript per million mapped reads) for each transcription region was calculated to quantify its expression abundance and variations, using RSEM software (Li & Dewey, 2011). RNAs and lncRNAs differential expression analysis were performed using DESeq2 software between two groups (Anders & Huber, 2010; Love, Huber & Anders, 2014). The genes/transcripts with the parameter of FDR below 0.05 and absolute fold change ≥1.5 were considered differentially expressed genes/transcripts. The statistical power of this experimental design, calculated in PROPER is 0.72 (Wu, Wang & Wu, 2015). Differentially expressed coding RNAs were then subjected to enrichment analysis of GO functions and KEGG pathways (Boyle et al., 2004; Kanehisa et al., 2008). A Protein-Protein interaction network was identified using String v10 and visualized using Cytoscape software (v3.7.1) (Shannon et al., 2003; Szklarczyk et al., 2015).

lncRNA-mRNA association analysis

The interaction between antisense lncRNA and mRNA was predicted by the software RNAplex (version 0.2) (http://www.tbi.univie.ac.at/RNA/RNAplex.1.html) (Tafer & Hofacker, 2008). lncRNAs in less than 10-kb up/downstream of a gene were identified to cis-regulators, and the lncRNA co-expressed with protein-coding genes were considered to have trans-regulation function. The target genes were further used for GO and KEGG enrichment analysis.

Data validation by qRT-PCR

Total RNA was extracted using RNAzol® RT RNA Isolation Reagent (GeneCopoeia, Rockville, MD, USA). Total RNA from each sample was reverse transcribed in a 25-μL reaction using Surescript™ First-Strand cDNA Synthesis Kit (GeneCopoeia, Rockville, MD, USA). The sequences of the primers used are indicated in Table S1. The NlRPS11 gene of BPH was used as an internal control gene. qRT-PCR was performed using a BlazeTaq™ SYBR® Green qPCR mix2.0 protocol (GeneCopoeia, Rockville, MD, USA). The 20-μL reaction volume consisted of forward and reversed primers (1-μL), BlazeTaq™ SYBR® Green qPCR mix (10-μL), ddH2O (6-μL), cDNA (2-μL). The selected genes were verified using an ABI 7900HT Fast Real-Time PCR System with a cycling temperature of 60 °C and a single peak on the melting curve to ensure a single product. Each sample’s relative transcript levels were obtained using the 2^−ΔΔCt method (Rao et al., 2013). At least three replicates were tested per sample. The obtained data were subjected to unpaired a two-tailed Student’s t-tests using GraphPad Prism software (version 8). Different stars were used to denote the significant variations at the P < 0.05 level. All the numerical data in figures are presented as means ± standard deviation (SD) of three independent replications.

High-throughput sequencing of BPH lncRNAs

To investigate the dynamic variation of lncRNAs in BPH when fed with TN1 (susceptible rice variety) and YHY15 (moderately resistant rice variety), the whole-transcriptome strand-specific RNA sequencing for two BPHs was conducted using three biological replicates. In total, after removing low-quality reads, more than 506 million clean reads (above 99.85% raw reads) were generated with Q30 value in all libraries exceeding 92.95% though high-throughput sequencing (Table 1). After removing the reads mapped to rRNA, the mapping rate of clean reads against the BPH reference genome ranged from 42.87–80.30%. About 60% of the reads were mapped to the exon region, and about 20% were mapped to the intron and intergenic regions, respectively (Table S2).

Identification of lncRNAs in BPH

In this study, after the known coding mRNAs (transcripts and their splices) were filtered out from the mapped reads, largely expressed transcripts with length ≥200 bp and exon number ≥2 were selected. The results indicated that 10,937 potential lncRNA transcripts were identified. Then, these sequences were further processed using algorithms CPC2 and CNCI to assure the nonexistence of protein-coding domains. At last, 3,112 highly reliable lncRNAs were identified (Fig. S1).

According to the position of the new lncRNAs relative to protein-coding genes on the genome, new lncRNAs can be divided into five categories. In our results, the lncRNAs included 638 sense lncRNAs, 470 antisense lncRNAs, 77 intronic lncRNAs, 228 bidirectional lncRNAs, 3,594 intergenic lncRNAs, and 560 others (Fig. 1). Additionally, in terms of FPKM, the transcriptional abundance of lncRNAs was significantly lower than that of mRNAs, but there was no significant difference between TN1 and YHY15 (Fig. S2).

Analysis of differentially expressed lncRNAs in BPH

To analyze the difference in expression of lncRNAs between TN1 and YHY15, the normalized expression of lncRNAs in two BPH were compared. In this study, the false discovery rate (FDR) below 0.05 and absolute fold change ≥1.5 were used to identify expression lncRNAs and mRNA differentially (Figs. 2A, 2B). A total of 157 differentially expressed lncRNAs (including 84 upregulation and 73 downregulation), and 675 differentially expressed mRNAs (including 281 upregulation and 394 downregulation) were identified (Fig. 2C), respectively.

Function of differentially expressed lncRNAs

To explore the potential ncRNAs functions, lncRNA and mRNA association analysis was performed in three ways: base complementary pairing of lncRNA and mRNA (antisense), lncRNA’s location in 10 kb upstream and downstream of its adjacent protein-coding genes (cis-regulation), and lncRNA’s correlation with co-expressed protein-coding genes identified predictable targets (co-expression). As the results showed, three antisense targets (Table S3), nine pair cis-regulation targets (include eight lncRNAs and eight mRNAs) (Table S4), and 3,972 pair co-expressed targets (include 156 lncRNAs and 643 mRNAs) were identified (Table S5).

For the three antisense targets, the expression profile of the lncRNAs and corresponding target genes indicated a similar trend. The expression levels of lncRNA MSTRG.12888.1 and it predicted that target ncbi_111044590 (homology of Dus4l, dihydrouridine synthase 4) was higher in TN1 than in YHY15, MSTRG.26118.1, and MSTRG.26199.1 with their target genes ncbi_111053346 (pilin gene-inverting like protein) and ncbi_120355505 (homology of TrwC, DNA relaxase/conjugal transfer nickase-helicase) showed lower expression level in TN1 than in YHY15 (Fig. 3A).

For the nine pair cis-regulation lncRNAs and targets, except lncRNA XR_005571984.1 and target gene ncbi_111050469 (zinc finger CCCH domain-containing protein) exhibited opposite expression trend (lncRNA XR_005571984.1 were higher expressed in TN1 than in YHY15, but mRNA ncbi_111050469 were lower expressed in TN1 than in YHY15), the other lncRNAs and their targets exhibited similar expression trend (expression trend were simultaneous upregulation or downregulation) (Fig. 3B). Ncbi_111058532 coding a serine/threonine-protein kinase, ncbi_111053346 coding a pilin gene-inverting like protein, ncbi_111053362 coding a F-box DNA helicase 1 (Fbh1), ncbi_111061470 is a homology of jerky protein, which facilitates Wnt signalling.

The function annotation was further performed for the 3,972 pair co-expressed targets, including inc 156 lncRNAs and 643 mRNAs (Fig. 3C, Table S6). First, Gene Ontology (GO) analysis was conducted to categorize these protein-coding genes. There were 23 classes of biological processes and these protein-coding genes were mainly enriched in “single-organism process”, “cellular process”, and “metabolic process”. Moreover, some essential stress resistance genes were identified as lncRNAs targets, including “response to stimulus”, “immune system process”, and so on (Fig. 4). These stress resistance genes homology may mainly related to the MAPK signalling pathway, including upstream regulation genes such as mitochondrial uncoupling protein 2 (ncbi_111047800), peroxiredoxin 1 (ncbi_111050813), toll-like receptor (ncbi_111064393), superoxide dismutase (ncbi_120349329), or downstream regulation genes such as defensin-like (ncbi_111045303), glutaredoxin-3 (ncbi_111048799), transferrin (ncbi_111050388), nemo-like protein (ncbi_111044309), serine/threonine-protein kinase (ncbi_111045417, ncbi_111048296, ncbi_111055048), P450 (ncbi_120356101) and so on (Table S7). These findings proposed that these differentially expressed lncRNAs might be involved in different biological processes, including stress resistance by regulating the expression of related protein-coding genes—significantly enriched cellular components, including cell, cell part, organelle, and membrane. The differential expression targets were mainly enriched in catalytic activity and binding in molecular function categories.

Moreover, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, these co-expression target genes were mainly enriched in metabolism, environmental information processing, organismal systems, cellular processes related pathways (Fig. 5, Table S8). Among these pathways, a part of categories enriched many genes, including metabolic pathways, signal transduction, transport and catabolism, and endocrine system. In metabolic-related pathways, arginine and proline metabolism, glutathione metabolism, and carbon metabolism categories were enriched 10 genes (Table S8). These findings suggested that lncRNAs have a significant effect on the adaptability of BPH when exposed to susceptible and resistant rice varieties.

Validation of differentially expressed lncRNAs and mRNA in BPH

qRT-PCR was used to confirm the expression profiles of the RNA-seq data in two BPH. The results indicated that among 10 selected different expression mRNAs and lncRNAs, about 80% of mRNA expression trends were consistent with RNA-seq data, except ncbi_111060469 and ncbi_111061470 (Figs. 6A and 6B). The qRT-PCR results of lncRNA indicated that, except MSTRG.26199.1, MSTRG.26220.1, and XR_005573932.1, expression trends of the other seven lncRNAs were consistent with RNA-seq and qRT-PCR, the gene number ratios attained 70% (Figs. 7A and 7B).