Comparative genomics and phylogenetic analysis of seven Ficus species based on chloroplast genomes

Plant material, DNA extraction, and sequencing

All seven Ficus species examined in this study (F. esquiroliana, F. pandurata, F. formosana, F. erecta, F. carica, F. hirta, and F. stenophylla) were sampled. Voucher specimen information and herbarium accession numbers listed in Table S1; all vouchers are deposited in the Herbarium of Guangxi Institute of Botany. Sampling and collection complied with the regulations of the Guangxi Zhuang Autonomous Region, and consisted of common, non-protected species collected from non-restricted locations.

Total genomic DNA was extracted from approximately 50 mg of silica-dried young leaf tissue using the modified cetyltrimethylammonium bromide (CTAB) method (Doyle & Doyle, 1987). DNA quality and concentration were assessed using agarose gel electrophoresis and a NanoDrop spectrophotometer. High-quality genomic DNA library construction and sequencing were performed by Beijing Gezhi Boya Biotechnology Co., Ltd. DNA was evaluated by Qubit, NanoDrop, and gel electrophoresis, then sonicated using a Covaris S220 to an average insert size of approximately 350 bp. Library construction was performed using the NEBNext^® Ultra™ II DNA Library Prep Kit. Sequencing was performed on an Illumina HiSeq system (Illumina, San Diego, CA) using a paired-end 2 × 150 bp configuration.

Chloroplast genome assembly and annotation

Raw sequencing reads were filtered to remove low-quality reads and adapters using Trimmomatic version 0.36 (Bolger, Lohse & Usadel, 2014). Clean reads were assembled into complete chloroplast genomes using NOVOPlasty v4.3 (Dierckxsens, Mardulyn & Smits, 2017), with the matK gene used as the seed sequence and the complete chloroplast genome of a related Ficus species (GenBank accession MN364706) supplied as a reference. The assembly accuracy was confirmed by mapping reads back to the assembled genomes. Annotation was performed with Plann version 1.1 (Huang & Cronk, 2015) and manually corrected in Geneious version 10.2.6 (Kearse et al., 2012). The circular genome map was visualized using OGDRAW version 1.3.1 (http://ogdraw.mpimp-golm.mpg.de/) (Greiner, Lehwark & Bock, 2019).

Repeat sequence analysis

Simple sequence repeats (SSRs) were detected using MIcroSAtellite version 2.1 (MISA, https://webblast.ipk-gatersleben.de/misa) (Beier et al., 2017) with default parameters (minimum repeat numbers: 10 for mononucleotide, 5 for dinucleotide, 4 for trinucleotide, and 3 for tetra-, penta-, and hexanucleotide repeats). Long repeat sequences (forward, reverse, palindrome, and complementary) were identified using REPuter version 1.1 (https://bibiserv.techfak.uni-bielefeld.de/reputer) (Kurtz et al., 2001) with a minimum repeat size of 30 bp and sequence identity greater than 90%. Tandem Repeats Finder version 4.04 (http://tandem.bu.edu/trf/trf.html) (Benson, 1999) was used for tandem repeat analysis, with parameters set to 2 for the alignment parameter match and 7 for mismatches and indels.

Codon usage analysis

Protein-coding genes were extracted, and codon usage patterns were analyzed using CodonW version1.4.4 (http://codonw.sourceforge.net/) (Sharp & Li, 1986) to compute the relative synonymous codon usage (RSCU) values. Heatmaps are drawn using TBtools version 2.2.2 (Chen et al., 2020), with RSCU = 1 defined as the neutral (no-bias) value and used as the central point of the color scale (white).

Comparative genome and sequence divergence analysis

Genomic structure visualization and global alignment among the seven Ficus cp genomes were performed using mVISTA (Frazer et al., 2004) in Shuffle-LAGAN mode, with F. hirta as the reference. The borders of the large single-copy (LSC), small single-copy (SSC), and inverted repeat (IR) regions were compared and visualized using IRscope version 3.1 (https://links.jianshu.com/go?to=https://irscope.shinyapps.io/irapp/) (Amiryousefi, Hyvönen & Poczai, 2018).

Nucleotide diversity (Pi) was calculated using DnaSP version 6.0 (Rozas et al., 2017) with a sliding window analysis (window size: 600 bp; step size: 200 bp) to identify highly variable regions among the cp genomes.

Phylogenetic analysis

Phylogenetic analyses were conducted on a concatenated nucleotide alignment of 79 complete chloroplast protein-coding genes (PCGs) from 26 Moraceae taxa (Table S2); Antiaris toxicaria (NC_042884) was used as the outgroup. Subgeneric assignments were adopted following established taxonomic and molecular treatments (Berg, 1989; Rønsted et al., 2008; Cruaud et al., 2012). Sequences were aligned with MAFFT version 7.450 (Katoh & Standley, 2013), and the best-fit model (General Time Reversible (GTR)+F) was determined using ModelFinder in IQ-TREE version 2.0.3 (Minh et al., 2020). Maximum likelihood (ML) analysis was conducted with IQ-TREE using 10,000 bootstrap replicates, and Bayesian inference (BI) was performed with MrBayes version 3.2.7 (Ronquist & Huelsenbeck, 2003) via four parallel Markov chain Monte Carlo (MCMC) runs of 1,000,000 generations each, sampling every 500 generations and discarding the first 25% of samples as burn-in. Trees were visualized using FigTree version 1.4.4 (Drummond et al., 2012).

General features

Complete chloroplast genomes for the seven Ficus species (F. esquiroliana, PQ526730; F. pandurata, PQ526731; F. formosana, PQ526732; F. erecta, PQ526733; F. carica, PQ526734; F. hirta, PQ526735; and F. stenophylla, PQ526736) were successfully assembled from Illumina HiSeq paired-end reads (150 bp) using NOVOPlasty software and deposited at National Center for Biotechnology Information (NCBI).

The Ficus plastomes consistently exhibited a typical quadripartite structure with conserved gene content, position, and orientation (Fig. 1), and their sizes, GC contents, numbers of genes, and other information are shown in Table S3. The cp genomes of Ficus species ranged in length from 160,340 bp (F. hirta) to 160,669 bp (F. formosana). The GC content was 35.9%. The total number of annotated genes in Ficus plastomes was 130, comprising 85 (79+6) protein-coding genes (PCGs), 37 (30+7) tRNA genes, and 8 (4+4) rRNA genes. Within these genomes, 18 intron-containing genes were identified (six tRNA and 12 PCGs). Among these intron-containing genes, 15 contain one intron (ndhA, ndhB, atpF, petB, petD, rpoC1, rps16, rpl2, rpl16, trnA-UGC, trnG-UCC, trnI-GAU, trnK-UUU, trnL-UAA, and trnV-UAC), while the other three genes (rps12, clpP and ycf3) contain two introns (Table S4).

Repeated sequences

A total of 506 complex repeat sequences, consisting of 210 interspersed repeats (64 forward, nine reverse, 134 palindrome, and three complementary) and 296 tandem repeats, were identified within plastome genomes using MISA, REPuter and Tandem Repeat Finder software, as described in the Material and Methods section (Fig. S1A). For interspersed repeats, the sequence length is mainly concentrated in 30–64 bp, regardless of the forward or palindrome repeats. As for the tandem repeats, most were in the range of 8–26 bp, and a four bp repeat was observed in F. formosana and F. carica within the ycf3 gene. These tandem repeats were mainly distributed in the non-coding LSC and SSC regions.

Microsatellites are small repeating units (one to six nucleotide) within a genome nucleotide sequence (Shukla et al., 2018). The high polymorphism rate of repeat sequences at the species level positions them as common molecular markers valuable for phylogenetic and population genetics studies (Zheng et al., 2020). An analysis of SSRs revealed a count ranging from 89 to 98 in the plastomes (Fig. S1B). The distribution of SSR types was predominantly mononucleotide repeats (59%), among which A/T repeats were the most numerous. Dinucleotides were the next most frequent at 22%, and tetranucleotides made up 9%, with other SSR types occurring at lower rates (Figs. S1C; S1D). This composition pattern suggests that mononucleotide repeats potentially play a greater role in genetic variability than other SSR types. Our finding that A/T repeats were most abundant is similar to reports from other studies (Munyao et al., 2020). Analysis of SSR locations further indicated that the majority were situated in non-coding, intergenic, and intron regions. These SSRs offer potential for developing specific markers crucial for systematic, evolutionary, and conservation studies in Ficus.

Codon usage

The patterns of codon usage and nucleotide composition contribute to establishing a theoretical foundation for the genetic modification of cp genomes (Mazumdar et al., 2017). A comprehensive analysis was conducted on the codon usage frequency of protein-coding genes within the chloroplast genomes of seven assembled species of Ficus. The study revealed 64 distinct RSCU values within the plastomes of these Ficus species. These genomes contain between 53,455 and 53,556 codons, with F. hirta containing the fewest and F. formosana containing the most. Among all the codons, leucine (Leu) emerged as the most prevalent amino acid, with a frequency ranging from 9.7 to 10.49%. This was followed by isoleucine (Ile), with a frequency between 8.64% and 9.35%. Conversely, tryptophan (Trp) was the least abundant, with a frequency of 1.21–1.35% (Fig. 2A). Methionine and tryptophan are each encoded by a single codon (ATG and TGG, respectively); therefore their RSCU values are 1.00 by definition and RSCU is not informative for assessing synonymous codon bias for these amino acids. Thirty codons were identified with an RSCU value greater than 1. Among them, except for UUG (Leu) and AGG (Arg), all codons ended with A or U(T) nucleotides (Fig. 2B). This observation suggested a preference for A and T as the terminal bases in codons.

Comparative analysis between seven Ficus species

The expansion and contraction of the IR region are major driving forces in the evolution of land plants (Kode et al., 2005; Raubeson et al., 2007; Yao et al., 2015; Abdullah Mehmood et al., 2020). In this study, we compared the positions of the LSC/IR junction (JL) and the IR/SSC junction (JS) across seven assembled Ficus plastomes (Fig. 3). The lengths of the IR regions were relatively uniform, ranging from 25,527 to 25,917 bp, indicating high conservation with minimal contraction/expansion across the seven Ficus plastomes. The JL (IR-LSC: rpl2 & rps19) boundary showed high similarity in seven Ficus plastomes, except for F. pandurata and F. stenophylla. At the IR-LSC boundary, the rps 19 gene crossed into the IR region by approximately 108 bp in F. pandurata and F. stenophylla, while it remained entirely in the LSC region in the other five species. The JS (IR-SSC: ycf1 & ndhF) boundaries were also highly similar in Ficus plastomes. The ycf1 gene crossed the IR-SSC border and extended into the IR region at approximately 1,116 bp. Concurrently, the ndh F gene consistently extended into the IRB region by 17 bp in all species, with the remainder of the gene located in the SSC. At the JSA (SSC-IRA) boundary, the ycf1 gene partially extended into the SSC region for all species, with its segment in SSC measuring 4,722 bp in six species and 4,728 bp in F. esquiroliana. At the JLA (IRA/LSC) boundary, the length of the intergenic spacer between the boundary and the trnH gene, which is located in the LSC region, exhibits significant variation among different species, ranging from 56 bp in F. pandurata to 437 bp in F. erecta.

We also analyzed the differences in the cp genomes among the seven species through global sequence alignment (Fig. 4). Using mVISTA for global sequence alignment, with F. hirta sequence and annotation file as references, we observed variations across different regions among the seven species. The results indicated that the seven chloroplast genomes are highly conserved. Furthermore, the alignment revealed higher sequence conservation in coding regions compared to non-coding regions, and in inverted repeat regions relative to single-copy regions (Fig. 4).

The nucleotide diversity (Pi, π) within Ficus plastomes ranged from 0 to 0.0141, with a mean value of 0.0022. Inverted repeat (IR) regions exhibited low nucleotide polymorphism, and the majority of variations were localized to the large single-copy (LSC) and small single-copy (SSC) regions (Fig. 5). While protein-coding regions demonstrated overall conservation across these plastomes, three gene regions (ccsA, ccsA - ndhD, and rpoB - trnC-GCA) displayed significantly elevated Pi values (>0.012). Notably, ccsA exhibited the highest divergence, reaching a Pi value of 0.0141 (Fig. 5). These polymorphic loci serve as promising candidate barcode sequences for phylogenetic inference and population genetic studies within Ficus.

Phylogenetic analysis

To ascertain the phylogenetic relationships among Ficus species, we employed 79 PCGs from 26 species to reconstruct phylogenetic trees. The results indicate that the topology generated by distance estimation using the best fit model (GTR+F) is consistent under two different reconstruction methods (ML and BI) (Fig. 6). All Ficus species formed a monophyletic clade, separate from the outgroup Antiaris toxicaria. Notably, all seven newly sequenced plastomes in this study are assigned to Subgen. Ficus (Table S2) and are distributed into two strongly supported subclades (Fig. 6). Within one clade, F. hirta, F. esquiroliana, F. erecta, F. stenophylla, F. pandurata, and F. langkokensis clustered together with 100% bootstrap support, indicating a close relationship among these species. Within the other clade, F. carica, F. polynervis, F. formosana F. heteromorpha, and F. ischnopoda clustered together, suggesting a close relationship. All branches had support rates above 70%.

Repeat Sequences and SSRs

We found a substantial number of complex repeat sequences and simple sequence repeats (SSRs) in the Ficus plastomes. The predominance of mononucleotide (A/T-rich) SSRs is a common feature in angiosperm chloroplast genomes (Munyao et al., 2020; Provan, Powell & Hollingsworth, 2001), likely reflecting underlying mutational biases and relaxed selection in non-coding regions. These identified SSRs and other repeat motifs represent valuable resources for developing high-resolution genetic markers, which can be applied in population genetics, phylogeography, and species identification within Ficus (Zhang et al., 2022a; Zhang et al., 2022b). The observed distribution of SSRs, mainly in non-coding regions, also supports their utility as informative hotspots for molecular marker development.

Codon usage bias

Analysis of codon usage revealed a strong bias toward codons ending in A or T, consistent with the compositional bias in angiosperm chloroplast genomes (Mazumdar et al., 2017). This bias may reflect selection for translational efficiency, mutational pressures, or a combination of both. The predominance of leucine and isoleucine codons, and the low frequency of tryptophan, is in line with findings from other plant cp genomes (Morton, 1998). Such codon usage patterns provide fundamental information for future studies of plastome gene expression and for potential transplastomic engineering in Ficus.

IR boundary dynamics

The expansion and contraction of inverted repeat (IR) regions are known to contribute to size variation and evolutionary novelty in plastid genomes (Wang et al., 2008; Xiong et al., 2009). Our comparative analysis of the seven Ficus chloroplast genomes revealed a generally conserved quadripartite structure, with IR region lengths ranging from 25,527 bp to 25,917 bp, indicating relative stability within this genus. However, detailed examination of the LSC/IR and IR/SSC junctions (Fig. 3) unveiled minor yet discernible variations, reflecting dynamic micro-evolutionary shifts. The presence of partial rps19 and rpl2 duplications at the respective boundaries suggests that both lineage-specific and potentially ongoing border shifts may occur within Ficus. However, the relatively limited IR boundary variation suggests a largely stable cp genome structure in this genus, compared to the more dramatic IR dynamics seen in some other angiosperms.

Potential molecular markers in Ficus

Global sequence alignment and nucleotide diversity (Pi) analyses demonstrated that, as expected, the IR regions were the most conserved while the single-copy regions especially non-coding and intergenic spacers exhibited higher sequence divergence. The identification of three highly variable regions (ccsA, ccsA - ndhD, and rpoB - trnC-GCA) with elevated Pi values highlights their potential as candidate DNA barcodes for Ficus systematics and population studies. These hypervariable loci may facilitate more precise species delimitation and phylogeographic analyses, potentially complementing or providing higher resolution than existing universal barcode regions such as matK and rbcL (CBOL Plant Working Group, 2009).

Phylogenetic analysis in Ficus

Phylogenetic reconstruction based on a concatenated alignment of 79 chloroplast protein-coding genes robustly resolved relationships among the seven newly sequenced Ficus species nested within a monophyletic Ficus clade (26 taxa sampled). The recovered topology, which is consistent between ML and BI analyses and shows high support across nodes, identifies two major subclades that correspond well with previously recognized divisions based on nuclear and plastid markers (Rønsted et al., 2008). In our sampling most newly sequenced taxa are assigned to Subgen. Ficus and are distributed across two well-supported subclades; a small set of Subgen. Urostigma taxa were included to provide taxonomic context.

While plastome-scale data markedly improve resolution relative to single- or few-locus approaches and can help resolve recent radiations, chloroplast genomes represent a single, typically uniparentally inherited locus and therefore may not fully capture complex reticulate histories (e.g., hybridization, introgression, incomplete lineage sorting) that are likely in diverse genera such as Ficus. Moreover, our sampling-seven new plastomes within a genus of >800 species-remains limited. To more rigorously assess subgeneric monophyly and the deeper relationships within Ficus will require denser taxon sampling across recognized subgenera and integration of nuclear (and mitochondrial) genomic datasets. Finally, organellar complementarity deserves mention. While chloroplast genomes provide a relatively straightforward single-locus phylogenomic perspective, mitochondrial genomes and nuclear datasets provide independent and complementary histories. The recent publication of the complete mitochondrial genome of F. hirta (Deng & Cai, 2025) highlights an opportunity for combined organellar analyses; for example, comparisons of shared repeats, transferred sequences, or conflicting topologies between organellar genomes may help detect events such as introgression or assembly artifacts. We therefore encourage future work that integrates chloroplast, mitochondrial and nuclear genomic data, coupled with denser taxon and geographic sampling, to provide a more complete picture of Ficus evolutionary history.