Next-generation sequencing based Cyperus niveus (Cyperaceae) complete chloroplast genome: A comparative analysis and phylogeny

Sample collection, DNA extraction, and sequencing

Fresh sample of C. niveus was collected on 12 September 2022 from hills of the National University of Sciences and Technology, Islamabad, Pakistan. The sample was brought to the Pakistan Museum of Natural History for confirmation and voucher number PMNH-041726 was assigned for future reference. Fresh leaves were silica dried and DNeasy Plant Mini Kit (Qiagen, Inc., Hilden, Germany) was used for whole genome DNA extraction. For quality and quantity checking 1% agarose gel electrophoresis and Multiskan GO (Thermo Scientific, Inc., Waltham, MA, USA) wereused. For sequencing, high-quality whole genome DNA was sent to Novogene, Hong Kong and sequenced from paired ends with 150 bp short reads and a 350 bp insert size using HiSeq2500 (Steemers and Gunderson, 2005). All other species used in this study are provided in Table 1 along with the accession number.

Table 1.

List of Cyperaceae species along with accession number.

Genome assembly and annotation

The chloroplast genome of C. niveus was de novo assembled according to methods of Abdulla et al. (2020) with Kmer values of 141, 131, 121, 111, 101, and 91 using Velvet 1.2.10 (Zerbino and Birney, 2008). The resulted contigs were then mapped to reference using Cyperusiria (accession No. MW123056.1) in Geneious R8.1 (Kearse et al., 2012) which generated whole chloroplast genome of C. niveus. The coverage depth analysis was performed using the Burrow-Wheeler aligner (Li and Durbin, 2010) and the boundaries of LSC, SSC, and IR (IRA and IRB) were determined using careful observation of overlapping contigs in Geneious R8.1 (Kearse et al., 2012).

Chloroplast genome of C. niveus was annotated using online server of GeSeq (Tillich et al., 2017) and presence of transfer RNA genes was further confirmed using tRNAscan-SE version 2.03 (Schattner et al., 2005) The circular map of chloroplast genome of C. niveus was drawn using OGDRAW (Organellar Genome Draw) (Greiner et al., 2019). For GenBank submission, the final annotated file was converted into GenBank format by using GB2sequin and after submission accession number OR418070 was assigned.

Genomic features of Cyperus niveus and comparative analyses

Genomic features of C. niveuswas analyzed in Geneious R8.1 (Kearse et al., 2012) and was compared with 10 other species of genus Cyperus. To access the degree of divergence and nucleotide variability (Pi) within the chloroplast genome, a group of 11 Cyperus species were aligned together using MAFFT v7.490 (Katoh and Standley, 2013) and manual adjustments were made using Geneious R8.1 (Kearse et al., 2012). An analysis using the sliding window approach was conducted in DnaSP v6.12.03 (Rozas et al., 2017) with the following settings: window length of 800 and step size of 100. The contraction and expansion of repeats regions in closely related four members of genus Cyperus (C.niveus, C. exaltatus, C. difformis, and C. aromaticus) was determined using IRscope (Amiryousefi et al., 2018).

Codon usage, amino acids frequency, and RNA editing sites

Amino acids frequency and codon usage of protein-coding genes of C. niveusweredetermined using Geneious R8.1 (Kearse et al., 2012). RNA editing sites weredetermined using Plant RNA Editing-Prediction & Analysis Computer Tool (PREPACT) version: 3.12.0 using NCBI BLAST 2.2.26+ (http://www.prepact.de/prepact-main.php) and Arabidopsisthaliana was selected as reference protein database (Lenz and Knoop, 2013).

Simple sequence repeats and oligonucleotide repeats

For simple sequence repeat (SSR) analysis Perl script MIcroSAtellite identification tools (MISA) v2.1, 2020-08-25 (Beier et al., 2017) was used at the following parameters; 10 for mono-, 5 for di-, 4 for tri-, 3 for tetra-, penta-, and hexanucleotide repeats. For oligonucleotide repeats REPuter (Kurtz et al., 2001) was used to detect forward (F), reverse (R), palindromic (P), and complementary (C) repeats at minimum repeat size set to 30 with 3 mismatches and minimum similarity set to 90%.

Substitution and InDel analysis

For substitution analysis, each region (LSC, SSC, IR) of C. niveus chloroplast genome was pairwise aligned with C. iria taking it as reference. using MAFFT v7.490 (Katoh and Standley, 2013) and then visualized in Geneious R8.1 (Kearse et al., 2012). For InDel events DnaSP (Rozas et al., 2017) was used.

Phylogenetic relationship family Cyperaceae

For phylogenetic relationship, chloroplast genome of 28 Cyperaceae species and 1 Poaceae species, Brachypodium-distachyon (accession No. NC_011032) which was used as outgroups. Complete chloroplast genome of 29 species was multiple aligned using MAFFT v7.490 (Katoh and Standley, 2013). All types of InDels were removed from alignment and only substitutions were used to construct Maximum Likelihood tree using online server of IQ-TREE at default parameters and substitution model GTR+FO was adapted automatically on basis of data type (https://www.hiv.lanl.gov/ content/sequence/IQTREE/iqtree.html) (Nguyen et al., 2014; Trifinopoulos et al., 2016). The ultrafast bootstrap value of 1,000 replicates wasapplied todataset (Hoang et al., 2018). The resulted Newick tree was brought into MEGA version 11.0.13 for visualization and beautification (Tamura et al., 2021). Different marker identifications were applied to differentiate species of each genus and outgroup.

Chloroplast genome features of C. niveus

Hiseq2500 paired end run generated 5.3 GB sequencing raw data that contained 17.86 million reads. The whole genome sequence data was used to de novo assemble the complete chloroplast genome of C. niveus with average coverage depth of 100×. The chloroplast genome of C.niveus had quadripartite structure with total length of 185,454 base pairs (bp). The chloroplast genome contained an LSC region of 99,339 bp, two IRs (IRA and IRB: each 38,021 bp) separated by SSC of 10,073 bp. The chloroplast genome of C. niveus was rich in adenine-thymine (AT) content constituting 66.7% and guanine-cytosine (GC) content was 33.3%. As a pivotal indicator of the evolutionary relatedness between distinct species, the GC content of C. niveus was similar to other species of genus Cyperus (Zhu, 2019; Ren et al., 2021; Yang et al., 2021). The GC content of three regions of chloroplast genome varies; the IR regions contain a high GC content (37.3%), compared to LSC (30.9%) and SSC (25.5%). One possible reason of higher GC content in IR regions can be the presence of rRNA and tRNA genes in these regions (Guo et al., 2018; Biju et al., 2019). Comparatively, the GC content of C. niveus was similar to C. rotundus while other species of Cyperus had lower GC content (Zhu, 2019; Ren et al., 2021; Yang et al., 2021). The detailed genomic features of C. niveus have been provided in Table 2. The complete annotated circular chloroplast genome of C. niveus is provided in Fig. 1.

Table 2.

Genomic features of Cyperusniveus.

Fig. 1.

Circular representation of the complete chloroplast genome of Cyperusniveus. Genes located within the circular structure undergo transcription in a clockwise direction, while those positioned outside the circle experience counterclockwise transcription. Genes possessing diverse functions are assigned distinct color codes. The inner region of the circle is shaded in a deeper gray to represent the GC content, while the lighter gray hue indicates the AT content. The circular illustration designates specific segments as LSC (large single copy), SSC (small single copy), and features inverted repeats known as IRA and IRB.

The chloroplast genome contained 140 unique genes, including 83 protein-coding genes, 39 tRNA genes, and 8 rRNA genes. Among these genes, 26 protein-coding genes, 8 rRNA genes, and 16 tRNA genes were present and duplicated in the IR regions, whereas 52 protein-coding genes and 21 tRNA genes were present in LSC region, and 5 protein-coding genes and 2 tRNA genes were present in SSC region. Furthermore, in chloroplast genes, 17 protein-coding genes and 8 tRNA genes containing introns were identified. Among these introns containing genes, clpP1 and pafl contained two introns while all other genes contain one intron. Seven introns containing genes (ndhA, ndhB, rpl2, rpl23, rps12, trnA-UGC, and trnl-GAU) were present in IR region and thus were duplicated. The other 11 introns containing genes were present in LSC and no intron containing gene were present in SSC. In these intron containing genes rps12 wasidentified as transsplicing gene as its 5′-end existed in LSC region as single copy while 3′-end existed in IR region as duplicate. The complete details of genes with their functions are given in Table 3.

Table 3.

Detail of genes, their types and function in chloroplast genome of Cyperusniveus.

Codon usage, amino acid frequency, and RNA editing sites

Deviation in codon usage holds significant importance in the modulation of gene expression and cellular processes. This phenomenon not only serves as a critical regulatory mechanism but also offers an extra avenue for investigating speciation and evolutionary dynamics at the molecular scale (Plotkin and Kudla, 2011). Codon usage analyses revealed that protein-coding sequence in C. niveus codes for a total of 20,997 codons. Of these codons, ATT showed maximum frequency (970) that codes for isoleucine while TGA showed minimum frequency which is a stop codon (Online Supplementary Material S1). The amino acid frequency analyses showed the abundance of leucine (10.90%) and isoleucine (9.20%) while cysteine frequency (1.10%) was in least. These observations were found consistent with other species of angiosperms (Biju et al., 2019). The PREPACT successfully anticipated 171 post-transcriptional RNA editing alterations within 83 protein-coding genes of C. niveus. Among these modifications, the matK and rpoB genes showed the highest RNA editing site predictions, with 15 sites each. Additionally, the ndhF gene displayed 13 predicted RNA editing sites, while ccsA, rpoC2, and ycf2 genes exhibited 9 RNA editing sites each. Moreover, the ndhD gene contained six predicted RNA editing sites. Notably, all these modifications occurred at the initial and secondary nucleotides of codons. The conversion rate at the second nucleotide was 1.5 times greater compared to the first nucleotide. The maximum number of amino acid shift was threonine to isoleucine, and leucine to phenylalanine. The detailed events of RNA editing sites are given in Online Supplementary Material S2.

SSR and oligonucleotide repeats

Dispersed repeats sequences are recognized for their pivotal role in driving chloroplast genome rearrangement and recombination processes (Weng et al., 2014; Wang et al., 2018). Furthermore, long repetitive sequences have proven to be valuable markers, extensively employed in the realms of plant evolution, comparative genomics, and phylogenetics research (Park et al., 2017). Perl script MIcroSAtellite identification tools (MISA) detected 146 total SSR and 28 compound SSR in chloroplast genome of C. niveus. The distribution of SSR was as follows; 93 in LSC, 21 in IR, and 9 in SSC. All compound SSR were distributed in LSC. These SSR include 90 mononucleotides (61.64%), 21 dinucleotides (14.38%), 11 trinucleotides (7.53%), 20 tetranucleotides (13.70%), 4 pentanucleotides (2.74%), and no hexanucleotides. Mononucleotides were made of 10–24 repeat units, dinucleotides were made of 5–13 repeat units, trinucleotides were made of 4–5 repeat units, tetranucleotides andpentanucleotides were made of three repeat units. In the case of complementary most of repeats belongs to A/T types (90), while AATCT/AGATT (1) were least (Online Supplementary Material S3). In closely related species C.esculentus 396 SSRs were identified and these SSRs encompassed 339 mononucleotides (constituting 85.6%), 22 dinucleotide (making up 5.56%), 10 trinucleotides (comprising 2.53%), 23 tetranucleotides (accounting for 5.81%), two pentanucleotide repeats (representing 0.51%), with no hexanucleotides detected (Ren et al., 2021). The SSRs pattern in C. niveus and C. esculentus were consistent with previous studies (Biju et al., 2019).

Two types of oligonucleotide repeats were found in chloroplast genome of C. niveus, which include 28 forward repeats and 21 palindromic repeats, whereas no inverse or complementary repeats were found. Most of repeats (91.83%) varied from 82–182 bp in length, and only four repeats were above 700 bp in length. In 49 oligonucleotide repeats, 13 repeats were only found in IR regions, five repeats were found only in LSC region, and one repeat wasfound only in SSC region. All other 24 oligonucleotide repeats were found in more than one region, either both in LSC and IR regions or both in SSC and IR regions. Furthermore, six forward oligonucleotide repeats were found both in IRA and IRB, and thus were duplicated. Among those 49 oligonucleotides repeats most (81.63%) were in intergenic spaces. In protein-coding genes, 6 oligonucleotide repeats were found (3 palindromic, 3 forward) in rpoB, while one forward repeat in each of following three genes; rpl32, rpoC2, and rps12 (Online Supplementary Material S4). In a preceding investigation, two categories of extended repeats were pinpointed within the chloroplast genome of C. esculentus. This encompassed 15 forward repeats and 34 palindromic repeats, with no detection of inverse repeats or complementary repeats. The majority of these repeats (93.8%) ranged from 64 to 385 base pairs in length, with only three exceeding 600 base pairs (Ren et al., 2021).

Substitution and InDel analysis

The different regions of chloroplast genome of C. niveus were pairwise aligned with C. iria taking it as reference. The 2,000 substitution was observed in complete chloroplast genome, 841 in LSC, 146 in SSC, and 1,013 in IR. A total of 290 InDel events were recorded; 183 in LSC, 25 in SSC, and 82 in IR (Table 4). Prior research findings indicate a greater prevalence of InDels and substitutions events within the LSC and SSC regions in comparison to the IR regions (Ahmed et al., 2012; Abdullah et al., 2019; Liu et al., 2019; Shahzadi et al., 2020). In the current study, InDels events were higher in LSC but substitutions events were more in IR as compared to LSC and SSC.

Table 4.

Type and number of substitution and InDel events in chloroplast genome of Cyperusniveus.

Out of these 2,000 substitutions 962 were transition and 1,038 were transversion type. The transition to transversion ratio (Ts/Tv) was 0.92 of complete chloroplast, while 0.95 in LSC, 0.92 in SSC, and 0.90 in IR. The current study’s Ts/Tv ratio wasconsistent with previous studiesin which authors reported Ts/Tv<1 (Ahmed et al., 2012; Liu et al., 2019; Abdullahet al., 2019).

Comparative analysis in genus

The chloroplast genome of C. niveus was compared with nine other Cyperus species whose chloroplast genome was already assembled. C. difformis has the largest genome size (187,808 bp), while C. rotundus has the smallest genome (182,986 bp). In the case of LSC individual size, C. glomeratus has the largest size (104,040 bp), and C. niveus has the smallest size (99,339 bp). The SSC and IRs have almost similar sizes in all species (Online Supplementary Material S5). The expansion and contraction of the IRs at their boundaries play crucial roles in the evolutionary processes, leading to variations in the length of plant chloroplast genomes (Kim and Lee, 2004). A thorough comparison was performed between the junctions of the IRs and two distinct single copy regions in the chloroplast genomes of four Cyperus species. This analysis encompassed regions such as LSC/IRA (JLA), LSC/IRB (JLB), SSC/IRA (JSA), and SSC/IRB (JSB), alongside the rearrangement of neighboring genes. Six genes (rps3 rps8, rpl14, rpl16, clpP1, and trnH) were present at junction of LSC and IRs. In C. niveus, clpP1 was entirely located in LSC region 124 bp closer to JLB, ndhA in IRB region 25 bp away from JSB, ndhl entirely in SSC, ndhF at JSA 926 bp in SSC and 1,279 bp in IRA, rps3 entirely in IRB, trnH in LSC 31 bp toward JLA. In C. exaltatus rps8, rpl14 was entirely in LSC region, rpl16 at the JLB 551 bp in LSC region and 905 bp in IRB region, ndhG in IRB region 97 bp toward JSB, ndhF entirely in SSC, ndhG in IRA 96 bp away from JSA, ndhl entirely in IRA, trnH in LSC 31 bp toward JLA. In C. difformis and C. aromaticus rps8, rpl14 both were in LSC, while rpl16 also in LSC but in C. difformis it was 124 toward JLB and 140 bp toward JLB in C. aromaticus. In both species ndhG were in IRB 87 bp toward JSB, ndhG in IRA 86 bp away from JSA, ndhl entirely in IRA. In C. difformis ndhF was entirely in SSC and trnH in LSC 31 bp toward JLA. In C. aromaticus psaC was entirely in SSC and while trnH was in LSC 29 bp toward JLA (Fig. 2). In genus Cyperus, variation was observed in genes presence and absence. In total 19 genes showed variation. In C. niveus, five genes were found unique as compared to other 10 Cyperus species. These include clpP1 gene, pafI gene, pafII gene, pbf1 gene, and rpl23 gene. The IR and SSC region were more divergent than LSC. The nucleotide variability of IR and SSC was higher (Pi > 0.3) than set threshold (Pi ≥ 0.2) The regions which showed greater divergence include tRNA (trnN-GUG, trnT-UGU, and trn-GUU) and genes (ccsA, ndhD, ndhE, ndhF, ndhG, and psac) (Fig. 3).

Fig. 2.

Comparative analysis of contraction and expansion of regions in Cyperus species. Colored boxes were employed above or below a central line to indicate adjacent border genes. The extent between these genes and the respective boundaries was quantified in terms of base lengths (bp). JSA (SSC/IRa), JSB (IRb/SSC), JLA (IRa/LSC), and JLB (IRb/LSC) serve as indicators of the specific sites where two corresponding regions of a genome meet and connect.LSC, large single copy; SSC, small single copy; IR, inverted repeats.

Fig. 3.

The sliding window analysis of complete chloroplast genome of Cyperusniveus. X-axis denotes the position of nucleotide and Y-axis denotes nucleotide variability (Pi) with threshold set to 0.2. The colored line below X-axis denotes different regions of chloroplast genome.

Phylogenetic relationship in family Cyperaceae

Complete chloroplast sequences hold significant value in unraveling phylogenetic relationships, particularly among closely related taxa. This is particularly relevant in cases involving recent divergence, rapid speciation, or slow genome evolution, where the extent of sequence variation may be limited (Daniell et al., 2016; Williams et al., 2016; Tonti-Filippini et al., 2017). Phylogenetic trees were constructed between C. niveus and 27 other species of Cyperaceae and 1 of Poaceae, Brachypodium distachyon as outgroup. The Maximum Likelihood tree showed that all species of the genus Cyperus are phylogenetically conserved as compared to other genera. The genus Carex showed maximum divergence which means that they are evolving rapidly and moving evolutionarily away from each other. The C. niveus in the family tree showed close resemblance with C. iria and C. exaltatus (Fig. 4). In this study, Cyperus emerged as single monophyletic group who has C4 photosynthetic cycle due to xeric conditions. A study was carried out to define phylogenetic relationship of giant genus Cyperuson the basis of nuclear ribosomal DNA (ETS1f) and plastid DNA (rpl32-trnL, trnH-psbA) reported that giant genus Cyperus consists of paraphyletic C3 genera and monophyletic C4 genera. The monophyletic C4 genera have nine nested segregate genera (Alinula, Ascolepis, Lipocarpha, Kyllinga, Pycreus, Queenslandiella, Remirea, Sphaerocyperus, and Voikiella) and they suggest these segregate genera should be incorporated as part of broad C4 genus Cyperus (Larridon et al., 2013). Our results wereconsistent with study as Cyperus emerged as single C4 monophyletic group. These findings were also congruent with the positions of other Cyperus species (Ren et al., 2021). Overall, these results are also consistent with classification of Cyperaceae as monophyletic, divided into two clades Mapanioideae and Cyperoideae on the basis of plastid rbcL and trnL-F (Muasya et al., 2009).

Fig. 4.

Maximum Likelihood (ML) tree of family Cyperaceae. The tree was constructed between 28 Cyperaceae and 1 Poaceae species Brachypodium distachyon (outgroup) using IQ-TREE at ultrafast bootstrap value of 1,000. Branch labels denote bootstrap confidence and colored identifier denotes different genus.