We collected fresh leaf samples of four Hevea species (H. benthamiana, H. camargoana, H. pauciflora, and H. spruceana) from the rubber tree germplasm resource nursery of the Chinese Academy of Tropical Agriculture Science (N 19°34′31.53″and E 109°31′17.97″). The total genomic DNA from the leaves was extracted using the CTAB method (Li et al., 2013). The genomic DNA was used for library construction. For next-generation sequencing, libraries with insert sizes of 350 bp were constructed and sequenced using the Illumina Novaseq6000 platform (Illumina, USA). For third-generation sequencing, 15 Kb fragment libraries were constructed by the PacBio RS II platform (Pacific Biosciences, Menlo Park, CA, USA). A total of 10.2–16.6 GB of raw sequencing data, of which the third-generation sequencing raw data was 5.5–10.1 GB ( Supplementary Figure S1 ), and a mitochondrial genome depth of coverage of 138.5–489.8× was obtained from the four Hevea species.

The short reads were filtered using Trimmomatic v0.38 (Bolger et al., 2014), and their de novo assembly was performed using SPAdes v3.5.0 with different K-mer parameters (Bankevich et al., 2012). GSAT software (He et al., 2022) was used to calculate the mitochondrial genome graph-based framework, the genomic structure maps of H. camargoana and H. spruceana mitochondria shows complex relationships between genomic sequences ( Supplementary Figure S2 ). The obtained PacBio reads were de novo assembled using canu v2.2 and PAMT (Koren et al., 2017; Bi et al., 2024) in hifi mode, with a genome size set to 2M, and default parameters for the rest, and the draft assemblies were polished by Pilon (Walker et al., 2014). During the running of Pilon, only paired-mapped reads were used to correct, which can eliminate nuclear mitochondrial DNA segments (NUMT) as much as possible. A total of 14–110 contigs were obtained by canu assembler, the longest was 309,883 bp in length, and the average length was 53,687 bp. The circularization process consists of the following steps: A. The establishment of coherent linkages among segments, B. Extensive and uniform coverage at the sequence assembly junctions, C. The wholeness of gene segments along with consistent and meaningful coding sequences, and D. An appropriate correlation between the repetition rates in repetitive segments and the depth of sequencing, ultimately leading to the effective and comprehensive assembly of the entire genome. The assembly process was summarized in Supplementary Figure S3 . At the same time, we calculated base depth and corrected the mitochondrial sequence using genotypes that matched the sequencing depth of mitochondrial sequence. The IGV software to analyze the reasonableness of sequencing read coverage and to determine if there were complete third-generation sequencing reads spanning the low-depth regions. For regions with low sequencing coverage or complex structures, validation was performed with Sanger sequencing to ensure the circular genomic sequence was accurate and complete. Specifically, the Primer Premier software was used to design the primers for each low sequencing coverage or complex structures region. The genomic DNA was used as the template to amplify the target DNA fragment. After purification and identification, PCR amplification and sequencing was performed. Subsequently, the mitochondrial genome was annotated by GeSeq (Tillich et al., 2017) and ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/) using the Hevea brasiliensis (AP014526) mitochondrial genome as the reference genome. After GeSeq annotation, the BLASTN and BLASTP alignment tools were used to identify and compare the mitochondrial sequences (Zhang et al., 2000), and the annotation results were manually modified. In addition, the tRNAs were annotated using tRNAscan-SE 2.0 (Lowe and Chan, 2016) and ARWEN (Laslett and Canback, 2008). Finally, the circular mitochondrial genome of Hevea species was mapped using Proksee (https://proksee.ca/).

Relative synonymous codon usage values were analyzed using MEGA v7.0.26 (Kumar et al., 2016). Codon usage heatmap was plotted using the R package pheatmap v1.0.12. Next, the AT and GC skews were calculated using the formulas: AT skew = (A − T)/(A + T), and GC skew = (G − C)/(G + C), respectively. In addition, the long repeats in the mitochondrial genome were analyzed using REPuter (Kurtz et al., 2001). Finally, the simple sequence repeats (SSRs) were identified using the MISA v2.1 software (Beier et al., 2017), with 10, 6, 5, 5, 5, and 5 as the minimum numbers of repeats for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively. Tandem repeats were identified using the program Tandem Repeats Finder (TRF) (Benson, 1999).

We selected the H. benthamiana mitochondrial genome, which has the highest similarity to the rubber tree mitochondrial genome, for paired repetitive sequence detection using BLASTN, and to identify potential recombinations from these repetitive sequences. We selected two of them for PCR validation. The primers at the ends of repeats were designed by Primer Premier3 v4.1 (https://bioinfo.ut.ee/primer3-0.4.0/) ( Supplementary Table S12 ). The PCR amplification system was 40 µL in total, consisting of 1 µL of template gDNA, 2 µL of forward and reserve primer, 20 µL of 2× Rapid Taq Master Mix, and 15 µL of ddH₂O. The PCR amplification program was 94 °C denaturation for 3 min, 35 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 60 s, and a 72 °C final extension for 5 min.

The mitochondrial genomes of the four Hevea species were searched by BLASTN against the chloroplast genomes (Zhang et al., 2000), using parameters “-evalue 1e-6 -outfmt 6”, and length > 100 bp to identify transfer fragments between the two genomes. The Intergenomic sequence transfer fragments between chloroplast and mitochondrial genomes were then mapped using Circos v0.69 software (Krzywinski et al., 2009).

The Deepred-mt (Edera et al., 2021) was used for the prediction of RNA editing sites of the mitochondrial genomes. It predicts C-to-U editing events based on the 40 nucleotides flanking a given cytidine.

Phylogenetic analysis was performed for 21 Malpighiales species using Passiflora edulis as the outgroup. The mitochondrial genomes of the 21 Malpighiales species, including the four Hevea species were downloaded from the Genbank. Next, the mitochondrial protein-coding genes (atp1, atp4, atp6, atp8, atp9, ccmB, ccmC, ccmFc, ccmFn, cox1, cox2, cox3, cob, matR, nad7, nad9, nd1, nd2, nd3, nd4, nd4L, nd5, nd6, prl10, prl16, rps1, rps12, rps3, and rps4) common to all the species were extracted to construct a phylogenetic tree. The mitochondrial gene sequences of Malpighiales species were aligned using Muscle v.3.8.31 (Edgar, 2004). A phylogenetic tree was built with a maximum likelihood (ML) approach using RAxML v8.1.5 (Stamatakis, 2006) based on the best-fit substitution model (GTR+I+G) determined by jModelTest v2.1.7 (Posada, 2008). Branch support was assessed using bootstrap resampling with 1000 replicates.

Raw mitochondrial genome sequence data were obtained from H. benthamiana (OR663908), H. camargoana (OR663909), H. pauciflora (OR663910), and H. spruceana (OR663911) and were successfully assembled as a single circular genome. The sizes of the four newly sequenced Hevea mitochondrial genomes ranged from 935,732 to 1,402,206 bp ( Figure 1 ; Table 1 ), with GC contents of 44.17–44.24%. Based on the comparative analysis, the four Hevea mitochondrial genomes had 34–35 unique protein-coding genes, 35–38 tRNA genes, and 6–13 rRNA genes ( Supplementary Tables S1–S4 ). Notably, the number of H. spruceana mitochondrial genes was less than that of the other three species, mainly because the copy number of some genes of H. spruceana (ccmB, cob, and matR et al.) was only a single copy ( Supplementary Figure S4 ; Supplementary Table S5 ). The four newly sequenced Hevea mitochondrial genomes have been submitted to the GenBank with accession numbers OR663908–OR663911.

The length of the Hevea mitochondrial genomes was approximately 5.8–8.7 times that of chloroplast genomes (MT333859, MN781109, NC_059798, and NC_059799). There were 17,294–46,552 bp long common sequence fragments between the chloroplast and mitochondrial genomes in the four Hevea species, including intergenic and gene regions ( Figure 2 ; Supplementary Table S6 ), which accounted for 1.9–3.5% of the mitochondrial genome and 10.7–28.9% of the chloroplast genome. The transfer sequences from the chloroplast to the mitochondrial genomes in the four Hevea species contained eight intact genes (psaA, rrn16S, tRNA-Val, rrn5S, rrn4.5S, tRNA-Arg, tRNA-Asp, and tRNA-Asn), intergenic spacers regions and partial gene sequences (psaB, ycf3, tRNA-Ile, ycf1, rrn23S, ycf1, tRNA-Ala, ycf2, atpE, ndhF, and rps16). Most of these transfer sequences are located in the IR region (50/79) of the chloroplast genomes ( Supplementary Table S6 ).

The length of most Hevea mitochondrial protein-coding genes (PCGs) is conserved ( Figure 3A ; Supplementary Tables S1–S4 ). In addition, the majority strand PCGs had a negative GC and positive AT skews, while the minority strand PCGs had positive GC and negative AT skews. Notably, the PCGs GC contents and AT skews were consistent across the Hevea mitochondrial genomes ( Figure 3B ; Supplementary Table S7 ).

The relative synonymous codon usage analysis of the mitochondrial genome revealed that among the four Hevea species, there were 64 codons encoding 20 amino acids ( Figure 3C ; Supplementary Figure S5 ; Supplementary Table S8 ). The total number of codons ranged from 12,292 to 18,395, with TTT as the most frequent (456–661 instances) codon, and TGC as the least (55–87 instances) among the PCGs in the Hevea mitochondrial genomes. The codon usage pattern was largely consistent in the four Hevea mitochondrial genomes.

Analysis of the Hevea mitochondrial genome repeats revealed 453–677 long repeats with a total length of 322,471–555,902 bp, accounting for 23.0–59.4% of the mitochondrial genomes ( Supplementary Table S9 ). Most of the long repeats were 30 to 50 bp in length, and the longest was 116,525 bp in H. benthamiana.

Tandem repeats, consisting of core repeat units are repeated many times in a tandem pattern. We detected a total of 61–93 tandem repeats in the mitochondrial genome of Hevea mitochondrial genomes with lengths between 5 and 76 bp ( Supplementary Table S10 ).

In addition, 98–133 SSRs were identified in Hevea mitochondrial genomes, with a total size of 1,306 to 2,030 bp ( Supplementary Table S11 ). Mononucleotide repeats were the most common, followed by di-, tri- and hexa- repeats, respectively. Most SSRs were located in the intergenic or intronic regions, except for 1–2 SSRs that were in the nad1 protein-coding region.

There were two repeat sequences that could potentially mediate the homologous recombination in H. benthamiana mitochondrial genome. We performed PCR amplification and Sanger sequencing to verify the presence of the minor conformation. The scheme of the primer design was shown in Supplementary Figure S6A and Supplementary Table S12 . Both sets of specific primers were successful in amplifying fragments of the repeat sequence fragments. Upon interchanging the reverse primers, we observed that they retained the ability to amplify these repeat sequence fragments, as shown in Supplementary Figure S6B . Based on these validation outcomes, we can infer the potential types of homologous recombination in the H. benthamiana mitochondrial genome, illustrated in Supplementary Figure S6A . Both conformations A and B are loop structures, each spanning a total length of 1,325,823 base pairs, and are differentiated by the orientation of their intermediate sequences, which are either forward or reverse. Due to the presence of direct repeats, conformation C consists of loops that form a two-chromosomal structure.

We predicted a total of 436–443 potential C-to-U RNA editing sites in Hevea mitochondrial PCGs ( Supplementary Table S13 ; Supplementary Figure S7 ). The predicted RNA editing sites of the mitochondrial genes of the four Hevea species are shown in Supplementary Figure S6 . Among these mitochondrial PCGs, we identified 37 RNA editing sites in nad4 and 35 for ccmB gene, which were the most frequent among these PCGs. In contrast, the gene ycf3 had only 0–1 C-to-U editing sites in Hevea mitochondrial genes.

Phylogenetic analysis was performed based on 29 common mitochondrial genes, from 21 Malpighiales species. Phylogenetic analysis revealed that most branches of the phylogenetic tree had high ML bootstraps values. Four well-supported clades (C1, C2, C3, and C4 Clade) were recovered within Malpighiale. C1 Clade included family Euphorbiaceae. C2 Clade included families Rhizophoraceae and Erythroxylaceae. The Calophyllaceae and Clusiaceae families formed C3 Clade. C4 Clade included family Salicaceae. Salicaceae species formed a monophyletic clade and diverged earlier. The five Hevea species (Euphorbiaceae) were clustered together ( Figure 4 ).

The Hevea mitochondrial genomes have consistent GC contents, codon occurrence and AT skews. However, there are significant differences in their lengths and sequence repeats. In addition, there are 17,294–46,552 bp intergenomic transfer fragments between the Hevea mitochondrial and chloroplast genomes. The evolutionary position of Hevea species was verified by phylogenetic analysis based on PCGs of the mitochondrial genome of 21 Malpighiales species. The findings from this study provide valuable genetic resources for further studies on Hevea species.