The transcriptomic sequencing data of WT (ZS6, Brassica napus variety “Zhongshuang 6”) and mutant transgenic rapeseeds under PEG, exogenous ABA treatment was obtained from a previous study in our lab (NCBI, BioProject accession number PRJNA1227215). Detailed information about the WT and mutant lines is described in detail in a previously published study (Luo et al., 2025). The specific treatments are as follows: WT and mutant lines were placed in dual-layer filter paper culture dishes (10 cm × 10 cm) containing water, 10% PEG, and 2 μM ABA solution, and incubated in a 25°C light-controlled incubator with three biological replicates. Samples were collected and sequenced at 0 h, 24 h, 48 h, and 72 h post-treatment.

The reference genome files and GFF file for B. napus were obtained from the Brassicaceae Database (http://brassicadb.cn/#/) (Chen et al., 2022), using version Brana_ZS_V2.0. Transposable element (TE) annotation of the genome was performed using the Extensive de novo TE Annotator (EDTA, v1.9.6) (Ou et al., 2019).

For transcriptome sequencing, raw data were processed with fastp (Chen, 2023) to remove adapter sequences and low-quality reads, resulting in clean data. The clean reads were aligned to the reference genome using Bowtie2 (Langmead and Salzberg, 2012), followed by BAM file sorting with Samtools (Danecek et al., 2021). Transcript assembly was performed with StringTie (Pertea et al., 2015), generating a GTF file for each BAM file. These GTF files were then merged using StringTie. The overlap regions between transcripts and TEs were extracted using Bedtools (Danecek et al., 2021) (parameters: -loj -wa -f 0.9). Transcripts overlapping more than 90% with TEs were considered TE-derived transcripts. The fastq files were re-mapped to the reference genome using STAR (Dobin et al., 2013), with both the reference genome GTF and the TE-derived transcript GTF provided separately as inputs for the –sjdbGTFfile option and –quantMode TranscriptomeSAM specified. The resulting transcriptome-based BAM files were used as input for RSEM to quantify the TE transcripts, while gene expression was quantified using featureCounts (Liao et al., 2014). The expression levels of transposable elements were quantified using TEspeX (https://github.com/fansalon/TEspeX) (Ansaloni et al., 2022). For the expression levels of TEs, TE transcripts, and genes, HTSeq (Putri et al., 2022) was used to convert read counts into FPKM (Fragments Per Kilobase of exon model per Million mapped fragments).

Genes, TEs, or TE transcripts with an average FPKM ≥ 1 across all samples were considered expressed in each abiotic stress experiment. Differential expression analysis for the comparisons PEG treatment vs. H₂O and ABA treatment vs. H₂O at three different time points (24 h, 48 h, and 72 h) was performed using DESeq2 (Love et al., 2014). Upregulated genes were defined as having a log2 fold change > 1 and adjusted P ≤ 0.05, while downregulated genes were defined as having a log₂|fold change| < -1 and adjusted P ≤ 0.05. The differential expression results for genes and TE transcripts across different groups were visualized using the UpSet R package (Conway et al., 2017).

The expression levels of TEs were log-transformed as log2(FPKM + 1), and principal component analysis (PCA) was performed using the prcomp function in R. The first two principal components were visualized with ggplot2 (Villanueva and Chen, 2019). The pheatmap package (Gu et al., 2016) was used to generate clustered heatmaps showing both sample correlations and TE expression levels, with sample correlations calculated using the cor function in R (the “method” parameter was set to “pearson”).

TE transcripts longer than 200 bp, with more than two exons and expression levels greater than 5 reads, and classified by gffcompare (Pertea and Pertea, 2020) with class codes “u,” “x,” and “i,” were selected for coding potential prediction. Three tools (PLEK (Li et al., 2014), CPAT (Wang et al., 2013), and CPC2 (Kang et al., 2017)) were used to assess the coding potential of these transcripts. Transcripts predicted as non-coding by both PLEK and CPC2, and with a predicted coding probability (Coding_prob) < 0.364 according to CPAT, were classified as lnc TE transcripts. Other transcripts that were predicted to have protein-coding potential by these tools were classified as protein-coding TE transcripts.

First, we downloaded the Pfam-A.hmm file from InterProScan (https://www.ebi.ac.uk/interpro/). Using HMMER (Johnson et al., 2010), we searched for Pfam-A.hmm within TE transcripts to identify the Pfam domains present in each TE transcript. We searched Pfam domains associated with transposases ( Supplementary Table S1 ). TE transcripts containing these transposase-related Pfam domains were classified as transposase.

Total RNA was extracted from different leaf samples of B. napus using the Plant RNA Purification Reagent (Invitrogen, Carlsbad, CA, USA). Genomic DNA contamination was removed with DNase I treatment (Invitrogen), and first-strand cDNA was synthesized from 2 μg of total RNA using the Omniscript RT Kit (Qiagen, Valencia, CA, USA). Gene-specific primers ( Supplementary Table S2 ) were designed using Primer Premier 5.0 (Lalitha, 2000) and verified for specificity through BLAST searches against the B. napus genome. Quantitative real-time PCR (qRT-PCR) was performed using the LightCycler 480 system (Roche), with amplification reactions conducted according to the instructions of the LightCycler 480 SYBR Green I Master Kit (Roche). The Actin7_141 gene was used as an internal reference (Lou et al., 2024). Each sample was analyzed in three biological replicates. Relative gene expression levels were calculated using the 2^–ΔΔCt method (Holzapfel and Wickert, 2007), based on the average threshold cycle (Ct) values for each sample.

In this study, we identified 212,800 TEs in the B. napus genome (Brana_ZS_V2.0), including 135,188 DNA transposons, 72,892 LTRs, 451 LINEs (Long Interspersed Nuclear Elements), and 4,269 MITEs (Miniature Inverted-repeat Transposable Elements). Our previous research revealed that BnaABI5 regulates ABA to mediate seed germination under drought conditions. Given that transcriptional activity in TE regions is associated with various biological processes, we quantified TE expression levels using RNA-seq data from 60 WT and mutant lines under different treatments. In the WT and mutant lines, 54,734 (25.72%) and 52,349 (24.6%) TEs were expressed, respectively ( Figure 1A ). Among them, DNA transposons had the highest number of expressed TEs (39,529 in the WT and 37,758 in the mutant lines), while LINEs had the fewest (146 in the WT and 162 in the mutant lines). There were 13,897 and 13,414 expressed LTRs, while there were 1,162 and 1,051 MITEs in the WT and mutant lines, respectively ( Figure 1B ). The chromosomal expression patterns of most TEs were similar across samples, but some genomic regions showed significant changes in TE expression between the control (treatment after 0h) and all treated groups. For example, the transcriptional abundance of DNA transposons was higher in the treated groups near the 22 Mb region (21406762- 22216053) on chromosome A10, the 18 Mb region (17630418- 18028290) on chromosome C08, and the 10 Mb region (10099783- 10969962) on chromosome C09 compared to the control ( Figure 1C ). The average transcript abundance of other transposons in all samples was LTR>MITE>LINE ( Supplementary Figure S1 ). Changes in epigenetic modifications due to abiotic stress may contribute to these differences in TE expression (Liang et al., 2021), and these altered TE transcripts may play a role in the seed germination of B. napus under drought stress.

Additionally, this study predicted the protein-coding potential of transcripts derived from TE regions. A total of 5,555 and 5,192 protein-coding TE transcripts were identified in WT and the mutant lines, respectively. For TE transcripts lacking protein-coding potential, we classified them as lnc RNAs, with 936 identified in WT and 874 in the mutant lines. We found that protein-coding TE transcripts were primarily located in intergenic regions. In WT and the mutant lines, 3,607 and 3,390 TE transcripts with protein-coding potential were located in intergenic regions, respectively ( Figure 1D ). Based on the annotated protein domains of these TE transcripts with protein-coding potential in intergenic regions, we identified 884 TE transcripts with Pfam protein domain annotations ( Supplementary Table S3 ), among which only a minority (65) were associated with transposases, likely due to the inherent nature of TEs as mobile genetic elements. Additionally, we identified 4754 TE transcripts without annotated transposase-related domains, suggesting the potential emergence of novel functions in these TE transcripts.

To explore the expression characteristics of TEs under different treatment conditions and time gradients, we systematically analyzed the TE expression patterns in the WT and mutant lines across various treatments and time points. PCA analysis of the TE expression profiles revealed similarities and differences among the samples under different treatments and time points. The PCA analysis showed that at the same time point, even under different stress treatments, the TE expression profiles of the samples were still highly similar. However, under the same stress treatment, the differences between samples at various time points were more significant, indicating that the influence of time factors on TE expression changes is greater than that of stress treatments ( Figure 2A ). This finding suggests that during seed germination, the expression levels of TEs are mainly regulated by developmental stages, while stress treatments have a relatively minor impact, a trend that is more pronounced in the mutant lines.

Further clustering analysis of TE expression revealed differences in the response of TEs to abiotic stress and developmental stages ( Figure 2B ). Most TEs exhibited low transcription levels. However, in both the WT and mutant lines, some TEs showed high expression only after stress treatments, while maintaining low expression in control samples (0 h after treatment). Clustering analysis based on TE transcription levels indicated that samples under different stress treatments exhibited similar expression patterns, and at the same time point, the samples from both genotypes also showed similar expression characteristics. This further confirms that during rapeseed seed germination, TE transcription patterns are more consistent across the same time point than across different stress treatments. To further confirm the features of TE expression patterns, we calculated and clustered the Pearson Correlation Coefficient (PCC) of TE expression levels (FPKM values) across multiple samples ( Figure 2C ). The results verified that the changes in TE transcription are primarily influenced by time factors, followed by stress treatment. These findings provide new insights into the biological functions of TE transcription during plant development and stress response, offering a foundation for further research.

TE transcripts were classified into 15,808 protein-coding TE transcripts and 1,739 non-coding TE transcripts. Functional enrichment analysis revealed that 4,832 protein-coding transcripts located in intergenic regions are primarily enriched in biological processes such as monocarboxylic acid metabolic process (34.5%), water-soluble vitamin biosynthetic process (29.9%), cellular nitrogen compound biosynthetic process (12%), and cellular modified amino acid metabolic process (6.0%) ( Figure 3A ). These findings indicate that TE regions transcribe transcripts not only encoding transposase-related enzymes but also produce transcripts associated with proteins involved in other biological functions. Furthermore, these protein-coding TE transcripts play an important role in the synthesis and metabolism of various compounds. This broadens our understanding of TE transcription biology and offers clues for uncovering TEs with novel biological functions.

This study comprehensively revealed the expression patterns of protein-coding TE transcripts under various stress treatments and dynamic time courses. Compared to water treatment, WT plants exhibited 38, 38, and 25 differentially expressed protein-coding TE transcripts after 24, 48, and 72 hours of ABA treatment, respectively. In contrast, the mutant lines showed 23, 8, and 18 differentially expressed protein-coding TE transcripts at the corresponding time points ( Figure 3B , Supplementary Table S4 ). Combined with our findings that mutant lines display higher seed germination rates under drought treatment, this suggests that the absence of BnaABI5 promotes seed germination under stress conditions, potentially by disrupting the regulatory networks between genes and between genes and non-coding RNAs (ncRNAs), which in WT plants may inhibit germination as a survival strategy during stress. This further indicates that these TE transcripts may be embedded within these intricate interaction networks and are regulated by the BnaABI5 mutation. Additionally, under PEG treatment, 39, 7, and 28 differentially expressed protein-coding TE transcripts were identified in WT plants at 24, 48, and 72 hours, respectively, compared to water treatment, while the mutant lines exhibited 32, 17, and 34 differentially expressed protein-coding TE transcripts at the same time points ( Figure 3C ). Notably, there was minimal overlap in the differentially expressed protein-coding TE transcripts identified under drought stress conditions between the two genotypes (32 transcripts), with the majority of differentially expressed protein-coding TE transcripts (65 transcripts) being specific to one genotype under a particular time point after the different stress conditions. This further demonstrates that the BnaABI5 CRISPR editing altered the drought response patterns in B. napus.

Among the 128 differentially expressed protein-coding TE transcripts ( Figure 3D ), five transcripts (MSTRG.43480.1, MSTRG.45876.3, MSTRG.59783.1, MSTRG.71359.1, and MSTRG.82357.1) encode proteins with transposase-related functional domains. Additionally, 27 TE transcripts, including MSTRG.47289.1, MSTRG.131575.1, and MSTRG.98597.2, exhibited high expression levels at 0 hours after ABA and PEG treatments but showed relatively lower expression at 24, 48, and 72 hours of treatment. In contrast, nine TE transcripts, such as MSTRG.10280.1, MSTRG.40701.1, and MSTRG.123670.1, displayed higher expression levels in the mutant lines at 24 hours of PEG treatment compared to other treatment groups, while their expression was significantly lower in WT plants after 72 hours of PEG treatment ( Figure 3E ). To further validate the differential expression results, we randomly selected six TE transcripts (MSTRG.40701.1, MSTRG.43480.1, MSTRG.65455.1, MSTRG.90684.1, MSTRG.129909.1, and MSTRG.10280.1) for qRT-PCR analysis ( Supplementary Figure S2 ). The qRT-PCR results showed that the relative expression levels and expression patterns of these TE transcripts across different samples were largely consistent with the RNA-seq data, supporting the reliability of our transcript quantification. These results indicates that the proteins translated from these nine differentially expressed protein-coding TE transcripts are primarily involved in drought response at 72 hours.

In addition to encoding proteins, TEs can regulate gene expression by transcribing lncRNA. In this study, 1,739 lnc TE transcripts were identified among 17,547 TE transcripts ( Figure 4A ). The length distribution of lnc TE transcripts is similar to that of protein-coding TE transcripts, with the majority falling within the 100–700 bp range. Among the protein-coding TE transcripts, 37 were longer than 3,000 bp, while 22 lnc TE transcripts exceeded this length ( Figure 4B ). These results suggest that lnc TE transcripts exhibit diverse length distributions and may share certain structural characteristics with protein-coding TE transcripts. Compared to water treatment, 9, 9, and 3 differentially expressed lnc TE transcripts were identified in WT plants after 24 h, 48 h, and 72 h of ABA treatment, respectively. 13, 3, and 5 differentially expressed lnc TE transcripts were identified in mutant lines under the same treatment and time points ( Figure 4C ). Similarly, compared to water treatment, WT plants exhibited 8, 1, and 1 differentially expressed lnc TE transcripts after 24 h, 48 h, and 72 h of PEG treatment, respectively. Meanwhile, 12, 3, and 12 differentially expressed lnc TE transcripts were identified in mutant lines under the same conditions. These findings suggest that lnc TE transcripts play distinct roles in response to abiotic stresses and that the mutation of BnaABI5 alters the expression patterns of these transcripts under stress conditions. These TE-derived lnc transcripts provide a new resource for identifying functional non-coding RNAs.

Given that lncRNAs regulate gene expression, we analyzed the regulation of protein-coding genes with differential expression under various treatments by lnc TE transcripts through co-expression analysis. Four co-expression networks were constructed for different genotypes and treatments ( Figure 5 , Supplementary Table S5 ). There were 16 differentially expressed TE transcripts ( Supplementary Table S6 ). Under PEG treatment, the co-expression network in WT contained 3 lnc TE transcripts and 44 DEGs, while that in the mutant lines included 9 lnc TE transcripts and 43 DEGs. Notably, the lnc TE transcript MSTRG.101546.5 was co-expressed with multiple DEGs in both WT (38 DEGs) and mutant lines (28 DEGs), suggesting that this transcript may act as a key regulatory factor involved in gene expression under PEG stress, with conserved functions across different genotypes. Under ABA treatment, the co-expression network in WT included 7 lnc TE transcripts and 27 DEGs. Among them, MSTRG.108925.4 was co-expressed with 16 DEGs, while MSTRG.109003.7 was co-expressed with 6 DEGs. Further analysis revealed that both lnc transcripts regulated BnA10g0418090.1 (encoding PHD finger protein ALFIN-LIKE 1) and BnC07g0801050.1 (encoding maltose excess protein 1, chloroplastic-like), indicating that lnc TEs may participate in ABA signaling or related biological processes by regulating these key genes. In the mutant lines under ABA treatment, the co-expression network contained 3 lnc TE transcripts and 8 DEGs. Among these, MSTRG.109003.7, MSTRG.61653.4, and MSTRG.68277.10 were all found to regulate BnC06g0773950.1 (encoding metal tolerance protein 9), suggesting that these lnc TE transcripts may be associated with metal ion homeostasis or resistance mechanisms.