Transposable Element Landscape in a Monotypic Species Barthea barthei (Hance) Krass (Melastomataceae) and Its Role in Ecological Adaptation

1. Introduction

Transposable elements (TEs) make up a substantial proportion of most plant genomes. The hallmark of active TEs is their ability to copy/cut themselves and integrate elsewhere in the host genome. According to their transposition mechanisms [1], TEs are grouped into two classes: Class I (Retrotranspons, REs) and Class II (DNA transposons). REs proliferates via an RNA intermediate called ‘copy and paste’ strategy, whereas DNA transposons propagate via a DNA intermediate using a ‘cut and paste’ strategy. Alongside this classification system, TEs can also be divided into autonomous elements or non-autonomous elements based on their ability to transpose on their own or dependence on others. Furthermore, each class is hierarchically organized into orders, superfamilies, families, and subfamilies based on DNA sequence structure features and phylogenetic relationships [2]⁠. In plant genomes, the order LTR (long terminal repeats) REs are usually the most prevalent type, and are subdivided into two superfamilies, Copia (RLC) and Gypsy (RLG), with size ranging from a few hundred base pairs to 25 kb [3]⁠. Other orders in the non-LTR REs are mainly composed of LINEs (long interspersed sequences), SINEs (short interspersed sequences), DIRS (DIRS-like elements), PLEs (Penelope-like elements). Class II elements in the plant genomes are grouped into the order TIRs (terminal inverted repeats) or Helitrons. According the variations of target site duplications (TSDs) and terminal inverted sequences, the order TIRs are categorized into 6 superfamilies including Tc1–Mariner (DTT), hAT (DTA), Mutator (DTM), P (DTP), PIF– Harbinger (DTH), CACTA (DTC), respectively. With two terminal inverted repeated sequences but truncated or deleted internal transposase (Tase), MITEs (miniature inverted-repeat transposable elements) are another non-autonomous group of common DNA transposons in plant genomes. MITEs are usually short in length (less than 500 bp) but high in copy number (tens of thousands), such as the Stowaway family of MITEs and the Tourist family of MITEs in rice. However, according to the unified classification system [3]⁠, the designation of MITEs has no descriptive power in a taxonomic sense due to shared structural features with different origins. With a typical 5’TC and 3’CTRR motif and a short hairpin structure lying 15 ~ 20 bp before the 3′ end, the other order class II elements, Helirtrons, replicate via a unique rolling-circle mechanism and frequently duplicate passenger genes.

Once considered “junk DNA” or “parasites” in the host genome, TEs are now increasingly recognized as major determinants of genome architecture and organism evolution [4,5]. For example, genome size variations among different plant species are mainly attributed to their TE content, ranging from less 3% of the minute genome of the carnivorous bladderwort plant Utricularia gibba (77 Mb) [6]⁠, to over 70% of the giant genome of Norway spruce (Picea abies, 19.6 Gb) [7]. Counter-intuitively, the expectation that there is a positive correlation between genome size and TE diversity only held true when the genome size was less than 500 Mb. Beyond this threshold, this correlation no longer exists [8]⁠. Despite having low transposable element diversity at the family level, the Norway spruce genome was dominated by a relatively small family number of long terminal repeat (LTR) elements with high copy numbers, which were estimated to range in age from 5 to 60 million years ago (MYA). In contrast, medium-sized genomes with high transposable element diversity usually exhibited high turnover rates, as evidenced by the vast majority of TEs being less than 5 million years old [7,9]⁠. Why TE diversity was so variable among different species? Why have some transposable element (TE) superfamilies/families been so successful in certain lineages? So far, there are still few clues as to how environmental and genetic factors have modulated the accumulation and diversification of TEs across species [10]⁠. With the increasing availability of plant whole-genome sequences, it is promising that we can characterize TE landscapes across a broader phylogenetic scale and untangle the causal mechanisms and processes underlying the observed variations.

Transposable elements are an important source of genetic variation within and across species, contributing to phenotypic innovations and adaptations through mechanisms like domestication, exaptation, host gene regulation, retrogene formation, and enhanced genomic plasticity [11,12,13]⁠. Typically, TEs are rendered inactive in the genome through various silencing mechanisms like DNA methylation, histone modification, small RNA-based silencing, chromatin modifications [14]⁠. However, these silenced TEs can become activated in response to stresses or developmental cues [15,16]⁠. Diverse TE elements have been found to regulate a wide range of agronomically or ecologically relevant traits, such as flower color variation in morning glories [17]⁠, fruit coloration in apple [18],⁠ sex determination in melon [19]⁠, drought tolerance in maize seedlings [20]⁠, and photoperiod sensitivity in maize [21]⁠. While these TE-driven trait variations have been primarily reported in crops and model species, the impact of TEs may be underestimated in natural populations, where a wider range of abiotic and biotic factors can influence their expression. With the increasing availability of non-model plant genomes, the hidden diversity of TEs is likely to be further revealed.

The evergreen shrub species Barthea barthei is the sole representative in the monotypic genus Barthea (Melastomataceae). It is exclusively found in southern mainland China, specifically in Hunan Province, Guangdong Province, Guangxi Province, and Fujian Province, as well as in Taiwan Province [22]⁠. Two varieties were previously recognized based on the width of the capsule wings, but recent population genetic analysis did not provide support for this distinction [23]⁠. In contrast to the uniform red flowers found in the Melastoma genus, the flower colors of B. barthei exhibit polymorphism among populations, with two distinct colors observed in natural populations spanning altitudes between 400 and 2500 meters (Figure 1A). The red-flowered populations tend to occupy open areas within the forest or mountain-tops, while the white-flowered populations are typically found in the understory forest, where light resources are limited. Based on their adaptation to varying light intensities, we designated these two population types as the 'red ecotype' and 'white ecotype,' respectively. In this study, we utilized de novo genome assembly of B. barthei and performed RNA-seq sequencing on the two ecotypes with the aim of characterizing the transposable element (TE) landscape of this unique species. Additionally, we sought to determine the role of transposable elements in the adaptation to contrasting light environments.

2. Materials and Methods

2.1. Plant Material and Sequencing

Flower buds of the two ecotypes of B. barthei were collected in the field. One ecotype with red petals were sampled from Erhuangzhang Nature Reserve, Yangcun county, Guangdong, China (EHZ population, E111°25′22.59″, N21°53′14.37″, alt.611m) and red ecotype from Lianhua Mt. Huidong county, Guangdong, China (LHS population, E115°13'56.10", N23°3'42.43", alt. 1280m). Flower buds for each ecotype were frozen using the Liquid-Nitrogen in the field, and taken back for total RNA isolation using the RNAprep Pure Plant Kit (TIANGEN Biotech Co. Ltd., Beijing, China). The quality of total RNA was determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The qualified RNAs were subject to paired-end library construction using Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA), and sequenced on the Illumina Hiseq2500 platform using 150bp paired-end reads.

2.2. Transposable Elements Annotation of the Genome B. barthei and Phylogenetic Analysis

The published chromosome-level genome assembly of B. barthei were obtained using Pacific Biosciences single-molecule real-time sequencing on PacBio Sequel II platform and high-throughput chromatin conformation capture mapping [24]. Using the Extensive de novo TE Annotator (EDTA) package [25], we curated a non-redundant TE library for the B. barthei and annotated the genome. For the predicted LTR-RTs or non-LTRs in the categories of unknown superfamilies, a de novo classification tool with convolutional neural network implemented in the package DeepTE was used to classify them with a default probability threshold of 0.6 [26]⁠. Using RPSBLASTN v.2.11.0+, all retroelements were searched against NCBI’s conserved domain database (CDDs, http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) [27]⁠, and only those retroelements hit to Pfam00078 (RVT_1, Ty3/Gypsy) or Pfam07727 (RVT_2, Ty1/Copia) with a significance E-value of 0.001 were extracted. To classify these filtered retroelements at the family level, one RT domain protein sequence from each family in the REXdb (Viridiplante v3.0) [28]⁠ were randomly selected. For those extracted nucleotide sequences, the putative amino acid sequences were determined using an ORF finder program (RRID:SCR_016643). Only those amino acid sequences with size over 50 were aligned using MAFFT v7.407 with default settings [29]. Prior to maximum likelihood phylogeny inferences, amino acid substitution models were determined using ModelTest-NG v0.1.7 [30]⁠. Using RAxML v8.2.12 [31]⁠, the phylogeny inferences based on RT domains were implemented with 100 nonparametric bootstrap replicates. Customer scripts used in the data analysis were deposited at Github project (https://github.com/altingia/Barthea_TE_manuscript).

2.3. Insertion Time Calculation of Intact LTRs and History of TEs Proliferation Inference

The two long terminal repeats for LTR-RT were identical at the initial insertion moment, thus the divergence between them can be used to approximate the insertion time for an intact LTR-RT. Using the LTRpred v 1.1.3 tool [32]⁠, which is optimized for detecting intact LTR-RTs, the insertion ages were calculated based on a mutation rate of 1.3E-8 substitutions per site per year [33]. Additionally, the Kimura genetic distances between individual TE family copies and their consensus sequences were retrieved from the RepeatMasker output using the ‘parseRM.pl’ script. This information was then used to summarize and visualize the accumulated coverage of TE classes and subclass across different age ranges, providing insights into the proliferation history of transposable elements, as detailed in custom R scripts (https://github.com/altingia/Barthea_TE_manuscript /03.TE_islands/).

2.4. Relationship between Transposable Elements and Adjacent Genes and Transposable Element Islands in the Host Genome

Using TSSFinder [34], the Barthea specific model for transcription start site (TSS) were trained based on the initial model of Arabidopsis, then TSSs with TATA-box motif in the core promoter region in the upstream of B. barthei genes were predicted, only those TE elements located less upstream 1000 bp to the corresponding TSS were considered to be promoter insertion. We developed a versatile pipeline to identify the relative position and distance between each transposable element and its neighboring genes. We classified the relationship between a TE and its nearest gene into following categories: ‘upstream (promoter)’, ‘5-overlap’, ‘5-UTR’, ‘CDS’, ‘intron’, ‘3-UTR’, ’3-overlap’, ‘downstream’. For nested TE elements, we further determined the spanned regions by integrating the coordinates of the TE into the coordinates array of the corresponding gene features, which were recorded in a hash table. This approach allowed us to precisely map the positional relationships between TEs and gene features (https://github.com/altingia/Barthea_TE_manuscript/04.TE_gene_relationships).

We summarized the frequency of different categories, along with their associated distances or overlap lengths across various classes of TE elements. To assess TE abundance, we analyzed the genome using 500 kb sliding windows with a step size of 100kb along each linkage group. Windows with over 50% TE content coverage were designated as TE islands, and consecutive windows meeting this criterion were merged. Gene ontology (GO) terms and KEGG pathways for genes within TE islands, as well as for the entire genome, were assigned by aligning with the KOBAS v3.03 databases [35] with an E-value threshold of 1e-3. These annotations served as the foreground and background in enrichment analysis. The significance of the results was evaluated using the Benjamini-Hochberg (BH) test and False Discovery Rate (FDR) correction methods [36].

2.5. Transcriptome Profiles for Transposable Element and Genes between White and Red-ecotypes

Raw RNA-seq reads of flower buds, with three biological replicates for each population, were filtered and trimmed using Trimmomatic v0.39 [37]⁠ with default parameters (ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 TOPHRED33). The trimmed clean reads were then mapped to the reference genome assembly using the aligner STAR 2.7.9 [38]⁠ allowing multiple mapped reads with parameters settings as “ --outFilterMultimapNmax 100 -- winAnchorMultimapNmax 200 --outSAMtype BAM SortedByCoordinate”, in accordance with recommendations for TE analyses. For both TEs and genes were normalized to Counts Per Million (CPM, DESeq_default). Differential expression analysis between ecotypes was performed using the TEtranscripts v2.0.3 package [39]⁠ with the following parameters: “-stranded no -mode multi -p 0.05 -i 10”. In this analysis, only genes or TEs with reads mapping along the entire length were considered to be transcribed.

2.6. Transposable element insertion polymorphisms among B. barthei populations

According to the recommendations for specificity and sensitivity [40], we used SPLITEREADER pipeline beta-1.2 [41] and TEPID pipeline v0.10 [42] to identify non-reference insertion sites and reference absence variants respectively. For SPLITREADER pipeline, it identifies potential non-reference insertion sites individually. These sites are then merged and filtered by TE family across the population, keeping sites where at least one individual show 3 or more supporting reads. The process examines negative coverage patterns in nearby regions (100bp upstream and downstream) for both presence and absence variants. For presence variants, true non-reference sites show reduced coverage compared to surrounding areas. Similarly, true absence variants show coverage drops at their edges. Low-coverage genomes that lack sufficient evidence for either insertion presence (positive coverage) or absence (negative coverage) are marked as NA, as neither pipeline can make reliable calls for these cases. Single nucleotide polymorphisms (SNPs) were identified using the GATK pipeline following best practices [43,44]. Sequencing reads were aligned with BWA-MEM, duplicates were marked using Picard, and base quality scores recalibrated with GATK's BaseRecalibrator. Variants were called using GATK HaplotypeCaller, which accounts for local haplotype structures. Post-calling, SNPs were filtered using default GATK parameters: QD < 2.0, FS > 60.0, MQRankSum < -12.5, ReadPosRankSum < -8.0, and QUAL < 30.0 to ensure high-confidence variants. A principal component analysis (PCA) using identified TIPs was performed using the ‘prcomp’ function from the stats package v3.2.3 in R v3.4.4 [45]. The first two eigenvectors were retained to create a two-dimensional plot. We analyzed linkage disequilibrium (LD) patterns for each transposable element insertion polymorphism (TIP) using PLINK v1.9 [46]. For each TIP, we calculated the LD coefficient (r²) between the TIP and 300 SNPs on each side (upstream and downstream), as well as pairwise r² values among all 600 SNPs surrounding the TIP. We then compared the proportion of TIP-SNP pairs showing high LD (r² > 0.4) with the proportion of SNP-SNP pairs showing high LD. Similar patterns were observed using a less stringent threshold of r² > 0.2. To account for regional variation in LD patterns, we classified TE variants into 'low', 'medium', or 'high' LD categories by comparing their ranked r² values (TIP-SNP) against the median ranked r² value of SNP-SNP pairs in the same region.

3. Results

3.1. Composition and Proliferation History of Transposable Elements in the B. barthei Genome

Using the EDTA package, we curated a lineage-specific, non-redundant library of 1081 consensus TE sequences for B. barthei. Subsequent RepeatMasker analysis of the genome assembly revealed a total of 64,866 TE copies, accounting for 16.76% of the genome based on a 40% divergence threshold from the consensus sequences (Table 1). Additionally, we identified 19,316 copies of other repeat elements, including low-complexity regions and microsatellites, comprising 2.1% of the genome. The initial EDTA output classified 5923 LTR and 358 non-LTR copies as unclassified at the superfamily level. Further analysis with DeepTE successfully assigned 5665 of these LTR/unknown copies to Ty3/Gypsy (3768 copies) and the Ty1/Copia (1897 copies) superfamilies (Table 1). The unclassified non-LTRs were further categorized into 224 PLEs, 95 LINEs, 38 SINE/tRNAs, and 1 DIR. In total, we identified 24,296 Ty3/Gypsy (8.82%) and 9602 Ty1/Copia (2.47%) copies, together comprising 11.29% of the entire genome assembly. The non-LTR fraction included 588 LINEs/L1, 8 LINEs/I, 278 Penelope copies, and 220 tRNA copies (Table 1).

Class II DNA transposons accounted for 5.28% of the genome, with TIRs contributing 3.44% and Helitrons 1.63%. MITEs made up an additional 0.21% (Table 1). Among the TIR DNA transposons, the Mutator superfamily was the most abundant, with 13,868 copies (2.7%), followed by CACTA (2,820 copies, 0.49%), PIF-Harbinger (747 copies, 0.12%), hAT (516 copies, 0.09%), and Tc1-Mariner (127 copies, 0.04%). Of the 64,864 TEs identified, 1108 copies, representing 2.1% of the genome, were found to be intact. These including 120 LTR/Copia, 220 LTR/Gypsy, 1 nonLTR/PLE, 426 TIR/DTM, 118 TIR/DTC, 81 TIR/DTA, 38 TIR/DTH, 15 TIR/DTT, and 89 Helitron. The remaining 14.66% of the genome consisted of truncated or fragmented TEs (Table 1).

The age distribution of all TEs showed a sharp peak at zero, indicating recent proliferation of TEs, followed by a steep decline and a broad distribution spanning from 3 to 10 million years ago (MYA) (Figure 1B). Among the superfamilies, the order of decreasing proliferation was as follows: LTR/Gypsy, DNA/Helitron, DNA/DTM, LTR/Copia, DNA/DTC, and DNA/DTH, with most elements being younger than 15 MYA. We identified 782 intact LTR elements, consisting of 188 Ty1/Copia, 381 Ty3/Gypsy, and 213 unclassified elements. These insertions occurred between 0.001 and 1.064 MYA (Figure 1C), with a notable concentration within the last 1.00 MYA. Insertion times for Copia and Gypsy largely overlap, showing minor peaks under 0.5 MYA, whereas unclassified LTRs had an older peak (Figure 1C). Overall, this subset of intact LTR elements, being relatively recent, suggests they might still be mobile and functional.

Figure 1. Transposable element landscape of Barthea barthei. (A) Photographs for two ecotypes of B. barthei; (B) The age distributions of different superfamiles measured by divergence with corresponding consensus sequences; (C) Insertion time distribution for intact LTR retrotransposons (LTR-RTs); (D) Classifications for LTR retrotransposons (LTR-RTs) based on phylogenies of reverse transcriptase (RT) domains for superfamily Ty1/Copia (above), superfamily Ty3/Gypsy (below); (E) Relative content of exonic and TE- derived sequences along the eight largest scaffolds of the B. barthei genome. Shown are DNA transposons (DNA), long interspersed nuclear element (LINE) and long terminal repeat (LTR) retrotransposons, as well as other TEs (other). The genome is well structured into TE- poor regions (“low- density regions,” LDRs) and TE- rich regions (“TE islands,” orange highlights).

Figure 1. Transposable element landscape of Barthea barthei. (A) Photographs for two ecotypes of B. barthei; (B) The age distributions of different superfamiles measured by divergence with corresponding consensus sequences; (C) Insertion time distribution for intact LTR retrotransposons (LTR-RTs); (D) Classifications for LTR retrotransposons (LTR-RTs) based on phylogenies of reverse transcriptase (RT) domains for superfamily Ty1/Copia (above), superfamily Ty3/Gypsy (below); (E) Relative content of exonic and TE- derived sequences along the eight largest scaffolds of the B. barthei genome. Shown are DNA transposons (DNA), long interspersed nuclear element (LINE) and long terminal repeat (LTR) retrotransposons, as well as other TEs (other). The genome is well structured into TE- poor regions (“low- density regions,” LDRs) and TE- rich regions (“TE islands,” orange highlights).

3.2. Retrotransposon Classification Based on RT Domains

By analyzing RT domains, we identified 342 copies of Ty1/Copia elements and 214 copies of Ty3/Gypsy elements, which matched the RVT_2 and RVT_1 domain, respectively. This analysis included 17 Ty1/Copia and 12 Ty3/Gypsy family representatives, which were used to constructed a maximum likelihood phylogenetic tree for the LTR-RT sequences using PROTGAMMAJTT amino acid substitution model and a rapid bootstrap method. Within the Ty1/Copia, we identified nine distinct lineages. The most common were Osser (74 copies), Bianca (73 copies), Angela (47 copies), followed by less prevalent families such as Ale (20 copies), Sire (10 copies), Tork (12 copies), and Tar (12 copies) (Figure 1D). For the Ty3/Gypsy elements, the majority were unclassified (108 copies), while the remaining copies were categorized into Athila (36 copies), Renia (28 copies), CRM (18 copies), Tekay (16 copies), Galadriel (3 copies), and Chlamyvir (1 copy) (Figure 1D).

3.3. Identification and Characterization of Transposable Element Islands in Gene Spaces

Our investigation revealed an uneven distribution of transposable elements across the 20 linkage groups. We identified 13 TE islands, ranging in size from 20.9 kb to 1.4 Mb, which account for 11.7% of the genome’s total TE content (Figure 1E, Table S1, Figure S1). The largest of these TE islands, located on chromosome 6, extends over 1.4Mb. These TE islands contained 1037 protein-encoding genes, which exhibited significantly lower gene density compared to regions with fewer TEs presence (Fisher’s exact test, p-value < 2.2e-16). However, these regions were hotspots for genes involved in the photosynthesis pathway (ath00195, FDR = 2.2E-09, Benjamini-Hochberg (BH) test) and the tryptophan metabolism pathway (ath00380, FDR = 0.04, BH test). Additionally, 29 GO terms, particularly those related to abiotic stress (e.g., 'response to hydroxyurea', 'response to aluminum ion', 'response to UV-B'), were significantly enriched in these TE islands (Table S2). This suggests that the genetic evolution of B. barthei has been shaped by its adaptation to the acidic, aluminum-saturated soils and variable light conditions in South China.

3.4. Genome-wide Survey of Transposable Element Insertion Preferences in the B. barthei Genome

A genome-wide survey was conducted to map transposable element insertion sites in the B. barthei genome using a flexible pipeline (http://github.com/alitniga/04.TE_gene_relationships) (Figure S2, Table S3). Of the 64,864 identified TE copies, 63,665 were mapped relative to neighboring genes, while 1199 were found unannotated scaffolds (Table 2). Analysis revealed that transposable elements (TEs) showed preferential insertion in upstream regions of genes (45.8%; 29,144 copies) compared to downstream regions (41.7%; 26,562 copies). For LTR-RTs, more higher proportions were located more than 5 kb away from the nearest genes in both upstream and downstream regions (Figure S3A, S3B). In contrast, DNA transposons (both TIRs and Helitrons) were predominantly found within 5 kb of genes (Figure S3C, S3D). In total, 17,991 upstream and 14,845 downstream TEs were located within 5 kb of neighboring genes. Additionally, 577 TE copies were identified within gene promoters, affecting 431 genes (Table 2). Approximately 26.1% of B. barthei genes had at least one TE copy within 1 kb - a lower frequency compared to Arabidopsis (36%) or maize (86%) [47,48]⁠. The most common TE types near genes were LTR/Gypsy, DNA/TIRs, LTR/Copia, DNA/Helitrons, and nLTR (Table 2). The study identified. Only a small proportion of genes with adjacent TEs showed significant GO term enrichment, including 'ADP binding '(GO0043531) and 'response to ethylene ' (GO0009723) (Table S4, Table S5). Additionally, 760 TE copies overlapped with gene 5' ends and 630 with the 3' ends.

Intriguingly, 6569 copies (10.3%) were found nested within host genes, classified as intronic, coding sequence (CDS), or untranslated region (UTR) TEs (Table 2). Single-intron nested TEs were most prevalent, with 5,454 copies distributed across 2560 genes. The gene Barthea36295, homologous to LUTEIN DEFICIENT 5 (a cytochrome P450 member in Vitis vinifera) [49]⁠, contained the highest number of copies (37), with 23 located in its third intron (Figure S4a). Genes containing single-intron nested TEs showed significant enrichment for ‘cytosol’ (GO:0005829, FDR = 3.3E-6) and ‘cytoplasm’ (GO:0005737, FDR = 3.2E-4) (Table S6).

For exon-nested transposons, 635 copies were identified within 227 genes, predominantly in UTR regions (190 in 5' UTR, 370 in 3' UTR) rather than CDS regions (75 copies in 43 genes). The gene Barthea44534, homologous to beta-glucosidase 12-like [50], contained the highest number of UTR-nested TEs (22 copies) (Figure S4b). Particularly, Barthea28554, homologous to the beta' subunit of RNA polymerase, contained the highest number of DNA/DTM copies in its sequences (Figure S4c). Notably, while eight genes lacked significant hits, 35 others resembled conserved proteins such as Ribonuclease, Reverse transcriptase, plant disease resistance polyprotein, heat shock protein, and others (Table S7). Genes containing exon-nested TEs were associated with cellular and stress responses, including protein folding chaperone activity and response to temperature stimulus (Table S8). Additionally, 480 TE copies spanning multiple introns/exons were found in 618 genes, linked to developmental processes and epigenetic modification (Table S9).

3.5. Expression Profiles of Genes and Transposable elements During Flower Bud Development between Contrasting Ecotypes

Initial validation of RNA-sequencing samples through hierarchical clustering demonstrated high reproducibility within ecotypes (Figure S5). Transcriptome analysis identified 40,538 expressed genes and 4495 transcribed TE copies (Figure 2, Table S10). Differentiation expression analysis, using stringent criteria (baseMean > 10, |log2FoldChange| ≥ 2, and Padj < 0.05), identified 2544 differentially expressed genes (DEGs), comprising 1357 up-regulated and 1187 down-regulated genes between the two ecotypes (Table S11).

KEGG pathway enrichment analysis of DEGs revealed three significantly pathways (Table S12). The photosynthesis-antenna proteins (KO00196, FDR = 3.23E-8, BH test) included 20 DEGs corresponding to 12 Arabidopsis proteins, including LHCA1~LHCA6, LHCB1~LHCB 6 (Table S13), which are essential for energy equilibrium under variable light [51]. The phenylpropanoid biosynthesis (KO00940, FDR = 1.50E-7, BH test) comprised 51 DEGs, homologous to 30 Arabidopsis proteins (Table S14, S15), playing crucial roles in lignin and flavonoid synthesis and light-induced stress response [52]. The cutin, suberine, and wax biosynthesis (KO00073, FDR = 2.09E-4, BH test) contained 15 DEGs annotated as homologs of various Arabidopsis proteins, including HXXXD-type acyl-transferase family protein, FAR2, CYP86A1, CER1, CYP704B1, CYP86B1 in Arabidopsis (Table S16), which are involved in producing protective cuticular waxes against environmental stressors [53,54].

Figure 2. The distribution of differentially expressed genes and transposable elements between the white petal ecotype and purple ecotype during flower bud development of Barthea barthei. The significant levels were determined using a log2FoldChange > 2 and an adjusted p-value of 0.05.

Figure 2. The distribution of differentially expressed genes and transposable elements between the white petal ecotype and purple ecotype during flower bud development of Barthea barthei. The significant levels were determined using a log2FoldChange > 2 and an adjusted p-value of 0.05.

Among the transcribed TEs, retrotransposons were predominant, with 1457 transcripts identified, including 1059 LTR/Gypsy, 206 LTR/Copia, 148 LTR/unknown elements, and 44 LINEs/unknown elements. These elements were frequently observed near genes or within single introns, with LTR-REs maintained an average distance of 5738.5 bp from adjacent genes (Table 3, Table S17). Additionally, 1938 DNA transposon transcripts (43.1% of expressed TE transcripts) were detected, comprising various types: 403 DNA/DTC, 626 DNA/DTM, 48 DNA/DTH, 45 DNA/DTA, 27 DNA/DTT, 140 MITEs, and 649 Helitrons. Of these, 855 transcripts (19.0%) were located within intron/UTR/CDS regions, while the remainder were positioned at an average distance of 1881 bp from their respective genes (Table 3, Table S17). The consistent detection of DNA transposon-derived transcripts across samples suggests these are not mere DNA contaminants, as previously suggested in other RNA-seq studies [55,56]. Rather, their distribution pattern-42.8% nested within genes and 57.2% located within 2000 bp of host genes – indicates probable passive co-transcription. Furthermore, 1083 unclassified TE transcripts were identified, with 391 nested within introns/UTR/CDS and 605 positioned near host genes at an average distance of 7995.8 bp (Table 3).

Differential express analysis of TEs revealed 343 significantly differentially expressed elements between ecotypes, comprising 90 upregulated and 253 downregulated TEs (Figure 2, Table S18). These included 76 LTR/Gypsy elements, 16 LTR/Copia elements, 9 LTR/unknown elements, 125 DNA transposons, and 117 elements of unknown type. Notably, 30 of these differentially expressed TEs were located near or nested within genes (Table S19). Statistical analysis demonstrated that overall, TE expression levels were significantly lower than genes expression levels (t-test, p= 1.91e-07).

3.6. Transposable Element Insertion Polymorphisms among Two Ecotypes

After stringent filtering, we identified 3859 high-confidence transposable element insertion polymorphisms (TIPs), including 3615 non-reference TE presence variants and 244 reference absence variants (Table S20). The majority of these TIPs were attributed to Gypsy (2,764) and Copia (231), along with DNA transposon superfamilies DNA/DTM (488) and DNA/DTA (128) (Figure 3A). Site frequency spectrum (SFS) analysis of TIPs, based on informative genomes at each site, showed that 1.62% of TIPs had a minor allele frequency below 0.05, though this proportion may be underestimated due to limited sample size (Figure 3B).

ecotype from YC population and purple ecotype from HD population of Barthea barthei. A. Superfamily components and proportions of non-reference transposon element insertion variants (Left) and reference absence variants (Right); B. The minor allele frequency (MAF) distribution of transposable element insertion polymorphism for B. barthei; C. Counts of TE variants with different minor allele frequencies within each genomic feature classified as coding regions (CDS), intergenic regions, intron and untranslated regions (UTR) regions; D. Principal component analysis for the samples of red ecotype and white ecotype based on non-reference TE insertion variants and reference TE absence variants; E. Pattern of transposable element variants and single nucleotide polymorphisms (TE – SNP) linkage, Boxplots showing distribution of minor allele frequencies for Boxplots showing distribution of minor allele frequencies for each LD category (high, mid, low); Proportion of TE insertions, TE deletions in each LD category (below).

Distribution analysis of TIP variants with estimated minor allele frequencies (MAF) revealed 1,699 in intergenic regions, 415 in introns, 264 in UTRs, and 95 in CDSs. Within each genomic category, variant frequency decreased consistently across increasing MAF bins from 0-0.1 to 0.4-0.5 (Figure 3C). Principal component analysis of TIPs demonstrated significant genetic differentiation between HD and YC populations, underscore the role of TIPs in distinguishing the genetic makeup of these populations (Figure 3D). Linkage disequilibrium (LD) analysis with nearby single nucleotide polymorphisms (SNPs) classified 51.4%, 39.0% and 9.6% of TE variants as having low, intermediate, and high LD, respectively. A positive correlation was observed between MAF and LD state, with higher MAF variants more frequently exhibiting high LD. This pattern was consistent across both TE insertions and deletions (Figure 3E). Notably, the classification of TE variants into high, intermediate, or low LD categories was largely consistent between TE insertions and TE deletions (Figure 3E), supporting the observation that common alleles tend to exist in high-LD states.

Twenty-nine fixed TIPs associated with 28 genes were identified between HD and YC populations (Table S21). These TIPs were distributed across various genomic regions: 22 in intergenic regions, 3 in introns, 2 each in CDS and UTR regions. The affected genes included homologs of various transcription factors (AP2/ERF, bHLH, MYB-like, C2H2) and key enzymes such as Homocysteine S-methyltransferase (HMT), alpha/beta hydrolase fold, and 2OG-Fe (II) oxygenase superfamily (Table S22). Three genes containing fixed TIPs showed differential expression between white and red color ecotypes. Barthea13541 (log|fold change| = -2.70, P-adj = 0.0003), containing a non-reference DNA/DTM insertion in YC population, is homologous to Pmr5/Cas1p GDSL/SGNH-like acyl-esterase family protein, which influences cell wall modification, plant-pathogen interactions, stress responses, development [57]. Barthea35147 (log|Fold Change| = 3.90, P-adj = 0.002), with a fixed DNA/DTM present in HD but absent in YC, is homologous to transcription factor AP2/ERF, which regulates plant morphogenesis, stress responses, hormone signaling, and metabolism [58,59]. Barthea36291 (log|Fold Change| = 4.05, P-adj = 0.003), containing a fixed DNA/DTA present in YC but absent in HD, is homologous to Cytochrome P450 (CYPs), which synthesizes various compounds crucial for membrane structure, hormones, UV protection, pigments, signaling, and volatile compounds that mediate biotic and abiotic interactions [60].

4. Discussion

4.1. Impact and Distribution of Transposable Elements in the Compact Genome of B. barthei

The B. barthei genome, spanning 235,034,099 bp, exhibits a relatively low proportion of TE components, consistent with its compact nature. This proportion aligns with the wide range of TE content observed across plant genomes, from 2.5% in the bladderwort U.gibba [6]⁠ to 85% in maize [61]⁠. Among plant genomes, the Ty1/Copia and Ty3/Gypsy superfamilies predominate, distinguished by their distinct arrangements of RT (Reverse Transcriptase) and INT (Integrase) sequences within the POL (Polymerase-like protein) open reading frames. The higher abundance of Ty3/Gypsy copies in Viridiplantae compared to Ty1/Copia [62] may reflect their distinct genomic distributions: older Ty3/Gypsy elements typically occupy heterochromatic regions, while younger Ty1/Copia elements are more prevalent in euchromatic areas where recombination suppression is reduced [63,64]⁠.

The B. barthei genome contains five well-recognized superfamilies of Class II DNA transposons, though notably lacks the P superfamily. The relative abundance of these elements varies considerably among plant species [65]⁠, with TIR and Helitron transposons showing particularly marked interspecific variation compared to the more stable prevalence of LTR retrotransposons. For instance, Arabidopsis harbors more Helitron copies than TIRs, and DNA transposon abundance can differ significantly even between closely related species [66]⁠. Although DNA transposons typically constitute a small fraction of plant genomes, their frequent proximity to genes and role as mutagens significantly influences genome structure [65]⁠. Mutator-like elements (MULEs), particularly active in plants, can increase mutation rates up to 50-fold [67]⁠. Pack-MULEs, can generate biased gene modifications through selective insertion and DNA capture [68]⁠. In B. barthei, Mutator elements emerge as the predominant DNA transposons, likely serving as important contributors to genomic diversity through their mutagenic activity.

This analysis reveals the complex dual nature of TEs in the B. barthei genome: they simultaneously serve as sources of genetic variation and potential threats to genomic integrity. Their ultimate impact on genome evolution and adaptation depends on both their regulatory mechanisms within the genome and the broader ecological and evolutionary context of the organism. The balance between promoting adaptability and maintaining genome stability appears to be particularly refined in B. barthei’s compact genome.

4.2. TE Island Facilitates the Adaptation to Tropical Forest for B. barthei

The TE islands in B. barthei, despite their sparse gene distribution, exhibit non-random gene enrichment patterns, particularly in pathways related to photosynthesis, tryptophan metabolism, abiotic stress tolerance. Photosynthesis, the fundamental process of converting carbon dioxide (CO2) into sugars, is notably sensitive to various abiotic stressors including ultraviolet (UV) radiation, light intensity fluctuations, and hypoxic or anoxic conditions [69]. Additionally, the tryptophan metabolism pathway plays a crucial defensive role against pathogenic infections, as demonstrated in both rice or Arabidopsis [70,71]⁠. This pathway’s enrichment in B. barthei’s TE islands may confer evolutionary advantages against pathogens and herbivores, particularly relevant given the Janzen-Connell hypothesis, which emphasizes their role in maintaining tropical forest diversity [72].

B. barthei stands out among Melastomaceae species for its extension into southern subtropical areas of South China. This distribution pattern aligns with MacArthur’s (1972) ecological paradigm, which posits that species' geographical ranges along environmental gradients are constrained by abiotic factors at one extreme and biotic pressures at the other [73]. To overcome these constrains, B. barthei has developed various adaptive traits, including enhanced photosynthetic capacity, polymorphic flower colors, and diverse growth forms. These adaptations are particularly evident in high-altitude populations, such as the purple-flowered ecotype from Huidong county, Guangdong province. This ecotype, adapted to an environment characterized by mountain mists, variable rainfall, strong winds, and intense sunlight, exhibits distinctive features including deep purple flowers, thick coriaceous leaves, and compact growth structure.

While the adaptive significance of TE islands has been well-documented in species such as the invasive ant C. obscurior [12], where they facilitate relaxed selection and enable genetic modifications that enhance chemical perception, learning, and insecticide resistance, their role in plant adaptation remains less understood. In B. barthei, the diverse phenotypes observed in flower color and leaf texture suggest that TE islands may play a crucial role in facilitating expansion into various ecological niches. Further research is needed to elucidate how transposable elements adjacent to candidate genes regulate these adaptive traits, potentially revealing the mechanisms by which TE islands contribute to B. barthei ‘s successful adaptation to tropical forest environments.

4.3. Transposable Element Insertion Preferences and Polymorphisms Associated with Ecological Divergence between Ecotypes

The distribution of transposable elements in the B. barthei genome reveals distinct insertion patterns and their potential role in ecological adaptation. While only 26.1% of B. barthei genes are associated with TEs – markedly lower than the 78% observed in maize [74] – both species show similar patterns of TE distribution, with most copies residing within 5 kb of genes rather than within their coding regions. Our genome-wide survey identified 3859 high-confidence TE insertion polymorphisms, predominantly contributed by Gypsy (2764) and Copia (231) elements, along with DNA transposon superfamilies DNA/DTM (488) and DNA/DTA (128). The distribution of these polymorphic insertions mirrors the general pattern of TE insertion preferences, with the majority occurring in intergenic regions (1699) followed by introns (415), UTRs (264), and CDS regions (95).

The non-random distribution of TEs reflects the interplay between natural selection and genetic drift [75]⁠. The decreasing frequency of TIP variants across increasing minor allele frequency (MAF) bins suggests ongoing selection pressures, particularly against potentially deleterious insertions within gene bodies. This selective pressure is further evidenced by the linkage disequilibrium (LD) patterns, where 51.4% of TE variants showed low LD with nearby SNPs, while only 9.6% exhibited high LD. The positive correlation between MAF and LD state suggests that beneficial or neutral TE insertions are more likely to persist and become common variants.

The impact of TE insertions varies considerably depending on their genomic location [76]. Insertions in regulatory regions (promoter, enhancers) can lead to three distinct outcomes: gene expression nullification, enhanced expression through new cis-regulatory sites, or silencing via repressive chromatin marks [77]. This regulatory potential is exemplified by cases such as the Copia-like element insertion into upstream of the ruby gene, a key MYB transcriptional activator, affecting anthocyanin production in blood oranges [78]⁠; and two independent insertions by a MITE and an LTR element in the teosinte branched1 (tb1) locus led to increased apical dominance and altered ear morphology during maize domestication [79].

The significant genetic differentiation between HD and YC populations revealed by TIP-based PCA analysis underscores the role of TE polymorphisms in ecological divergence. Particularly noteworthy are the 29 fixed TIPs associated with 28 genes between these populations, including three differentially expressed genes that may contribute to ecotype-specific adaptations. The fixed DNA/DTM insertion affecting Barthea13541 in the YC population potentially influences cell wall modification and stress responses through its homology to Pmr5/Cas1p GDSL/SGNH-like acyl-esterase. Similarly, the ecotype-specific insertions near Barthea35147 (AP2/ERF homolog) and Barthea36291 (Cytochrome P450 homolog) likely contribute to divergent morphogenesis, stress responses, and secondary metabolism between white and red ecotypes.

While intron-nested insertions constitute 8.6% of total insertions and have been linked to natural phenotypic variations in other plants - such as the double-flower variants in Japanese morning glory [80] and yellowhorn [81] – CDS insertions remain rare due to strong negative selection. The few surviving CDS-nested TEs, like the Copia insertion affecting salt sensitivity in soybean [82], typically represent ancient, conserved elements. This pattern of insertion distribution, combined with the fixed polymorphisms between ecotypes, suggesting that TEs have played a significant role in the adaptive divergence of B. barthei populations across different ecological niches.

5. Conclusions

Our comprehensive analysis of the transposable element landscape in B.barthei reveals several significant findings regarding genome organization and ecological adaptation. The compact genome of B. barthei (235 Mb) contains a relatively low proportion of TEs, with distinct patterns of distribution dominated by Ty3/Gypsy elements and Mutator-like DNA transposons. Despite their limited abundance, these TEs demonstrate a clear non-random distribution pattern, with most elements residing within 5 kb of genes, suggesting their potential regulatory roles in genome evolution and adaptation.

The identification of TE islands enriched in genes related to photosynthesis, tryptophan metabolism, and stress response provides compelling evidence for their role in ecological adaptation. These genomic features likely facilitated B. barthei's successful colonization of diverse habitats across tropical and subtropical regions of South China, as evidenced by distinct ecotypes such as the purple-flowered variant in Huidong county.

Analysis of 3,859 high-confidence TE insertion polymorphisms revealed significant genetic differentiation between populations, particularly between HD and YC ecotypes. The discovery of 29 fixed TIPs associated with 28 genes, including those involved in cell wall modification, stress responses, and secondary metabolism, suggests that TE-mediated genetic variations have contributed substantially to ecological divergence in B. barthei. The observed patterns of TE distribution and polymorphism, characterized by preferential insertion in intergenic regions and strong selection against CDS insertions, reflect the delicate balance between generating adaptive variation and maintaining genome stability.

These findings enhance our understanding of how TEs contribute to plant genome evolution and adaptation, particularly in the context of tropical forest environments. Future research should focus on elucidating the specific regulatory mechanisms by which TE insertions influence adaptive traits, potentially offering insights into the role of mobile genetic elements in plant speciation and ecological divergence.