Genome Sequencing of the Endangered Kingdonia uniflora (Circaeasteraceae, Ranunculales) Reveals Potential Mechanisms of Evolutionary Specialization.

Kingdonia uniflora, an alpine herb, has an extremely narrow distribution and represents a model for studying evolutionary mechanisms of species that have adapted to undisturbed environments for evolutionarily long periods of time. We assembled a 1,004.7-Mb draft genome (encoding 43,301 genes) of K. uniflora and found significant overrepresentation in gene families associated with DNA repair, underrepresentation in gene families associated with stress response, and loss of most plastid ndh genes. During the evolutionary process, the overrepresentation of gene families involved in DNA repair could help asexual K. uniflora reduce the accumulation of deleterious mutations, while reducing genetic diversity, which is important in responding to environment fluctuations. The underrepresentation of gene families related to stress response and functional loss of ndh genes could be due to lack or loss of ability to respond to environmental changes caused by long-term adaptation to a relatively stable ecological environment.

Introduction

Habitat destruction caused by changing climate and human activities has driven numerous plant species to endangered status (Yang et al., 2018). Current management of endangered species includes in situ and ex situ measurements (Nuijten et al., 2016). Owing to the rapidly growing number of threatened species, optimizing conservation practices based on evolutionary history and endangered mechanisms is imperative. In recent years, the development of genomic techniques has opened up new opportunities for species conservation. The main improvement is the enormous increase in loci allowing for a more detailed study of previous research questions (Nuijten et al., 2016). Specifically, genomic approaches provide the possibility to detect adaptive genetic variation across the entire genome (Kohn et al., 2006, Ouborg et al., 2010, Angeloni et al., 2012, Nuijten et al., 2016), which may help to address long-standing questions in conservation biology not fully resolved with traditional methods (Primmer, 2009, Avise, 2010, Frankham, 2010, Ouborg et al., 2010, Allendorf et al., 2010, Angeloni et al., 2012, Funk et al., 2012, Steiner et al., 2013, McMahon et al., 2014, Shafer et al., 2015; Garner et al., 2016, Benestan et al., 2016, Fuentes-Pardo and Ruzzante, 2017).

Plant lineages vary widely in their geographic distributions owing to numerous factors, e.g., ecological niche range, dispersal capacity, reproductive investments, and time since origin (Debussche and Thompson, 2003, Lavergne et al., 2004, Becker, 2010, Youssef et al., 2011), among which one crucial factor is adaptive capacity to respond to environmental changes. Lineages maintaining a small distribution are likely to possess low adaptive capacity to respond to geological and climatic changes at large scales, as well as habitat changes at small scales. For example, some asexual lineages, generating genetically and phenotypically identical individuals, are limited in their capacity to respond quickly to environmental fluctuation owing to low genetic diversity, which can lead to a status of endangerment or even extinction (Bell and Collins, 2008). In addition, species living in an equable environment for long periods might lack or lose the ability to defend against rapidly fluctuating environmental stress. Once the habitat is altered or destroyed, these species are incapable of colonizing new habitats. Shrinking habitats and low adaptive ability to new environments together lead some plants to endangered status. A recent study has shown that endangered species have more vulnerable genomes than that of congeneric non-endangered species (Hamabata et al., 2019); in fact, the endangered species exhibited significantly lower genetic diversity and proportion of duplicated genes and accumulated more deleterious variations than non-endangered species. Owing to limited genomic studies, genomic changes of endangered plants are largely unknown. To clarify which genetic factors contribute to vulnerability in endangered species, genome-wide genetic changes in functionality should be assessed (Hamabata et al., 2019).

Kingdonia uniflora Balf. f. and W.W. Sm. (Circaeasteraceae, Ranunculales), an alpine herb (diploid, 2n = 18), has a very narrow distribution (Figure S1). The habitat of K. uniflora represents an ecological environment of primeval forest with few disturbances. Specifically, K. uniflora is restricted to growing in high altitudes (∼2,800–4,000 m); cold, damp climates with deep humus; and usually under species of Abies. K. uniflora and Circaeaster agrestis Maxim. together constitute the early-diverging eudicot family Circaeasteraceae (Ranunculales) (Angiosperm Phylogeny Group, 2016). Previous estimates showed K. uniflora diverged from C. agrestis around 52 mya (Ruiz-Sanchez et al., 2012). Although no fossil record is known for K. uniflora, fossil fruits similar to those of C. agrestis have been reported from the mid-Albian of Virginia, USA (Crane et al., 1994, Drinnan et al., 1994, Sun et al., 2017). Remarkably, different from all other angiosperms, K. uniflora and C. agrestis possess an unusual dichotomous venation (Figure S2) similar to that found in ferns and Ginkgo (Sun et al., 2017). Old age for the lineage and dichotomous veins indicate an ancient relictual character for K. uniflora. Additionally, K. uniflora typically reproduces asexually, relying on rhizome systems to produce new individuals. These indications raise multiple questions, including: As an ancient relictual lineage, what history has K. uniflora experienced? As an endangered and habitat-dependent species, what genomic traits does K. uniflora possess? K. uniflora provides an ideal model to study the evolutionary mechanisms of ancient plant lineages that have an extremely narrow distribution and rely on a highly specialized habitat.

The ndh genes encode subunits of the thylakoid NADPH complex that mediate cyclic electron flow around Photosystem I and facilitates chlororespiration (Martín and Sabater, 2010). A series of studies suggests that the ndh genes can be dispensable under mild non-stressing environments (e.g., Casano et al., 2001, Martin et al., 2004, Rumeau et al., 2007, Martín and Sabater, 2010). A previous study (Sun et al., 2017) found rampant loss and pseudogenization of ndh genes in the K. uniflora plastome, which may indicate an association between the vulnerability of K. uniflora and genomic changes. In the present study, we provide a de novo genome sequence of K. uniflora using both Illumina and Pacbio sequencing technologies. We aim to investigate the evolutionary history of K. uniflora and reveal the potential mechanisms of its evolutionary specialization.

Results

Genome Assembly and Annotation

Genome size estimation using flow cytometry suggested a haploid genome size of 1,150 Mb for K. uniflora (Figure S3), whereas k-mer statistics indicated a similar genome size of 1,170 Mb, with very low heterozygosity (Figure S4 and Table S1). In the present study, we generated 236 Gb of Illumina reads and 106 Gb Pacbio reads with an N50 length of 12.876 kb (Table S2). A total assembly of 1,004.7 Mb (representing ∼86% of the estimated genome size), consisting of 2,932 scaffolds (scaffold N50 of 2.09 Mb; longest scaffold of 11.5 Mb) was achieved (Table 1). A total of 43,301 protein-coding genes were predicted (Table 1), among which 35,953 genes (83.03%) were functionally annotated (Table S3). In addition to protein-coding genes, various noncoding RNA sequences were identified and annotated (Table S4), including 1,124 transfer RNAs, 715 ribosomal RNAs, 125 microRNAs, and 1,751 small nuclear RNAs. The completeness of gene regions assessed by BUSCO (Benchmarking Universal Single Copy Orthologs) showed that 90.6% of the green plant single-copy orthologs were complete (Table S5).

Table 1

Genome Assembly of Kingdonia uniflora

Genome Features	Contigs/Scaffolds
Total length, bp	1,004,656,313
Total number of contigs	2,932
Longest length, bp	11,531,354
Length of N50, bp	2,099,369
Length of N90, bp	292,588
GC content, %	38.04%
No. of genes	43,301

We compared the draft genome of K. uniflora with the well-annotated genomes of the model plant Arabidopsis thaliana (Brassicaceae) and the Ranunculales species Aquilegia coerulea (Ranunculaceae). The genome size of K. uniflora is much larger than that of both references. The K. uniflora genome showed strong synteny with the genome of Aq. coerulea (Figure 1), but weak synteny with that of A. thaliana (Figure S5), which is not surprising given their placement in the angiosperm Tree of Life. The gene density in K. uniflora is lower than that in Aq. coerulea and A. thaliana, whereas the density of TEs (transposable elements) in K. uniflora was higher than that in Aq. coerulea and A. thaliana (Figures 1 and S5).

Figure 1

Comparative Analyses of Genomic Features between Kingdonia uniflora and Aquilegia coerulea

Tracks from inside to outside are collinearity between both genomes, number of chromosomes/contigs, gene density, GC content, and TE density.

We also compared our draft genome sequence with five other draft genomes of Ranunculales taxa, representing three of the seven Ranunculales families (Table S6); the quality of our assembly is comparable with that of all the five species, generating the longest N50 length and relatively fewer scaffolds. Comparatively, the genome of K. uniflora is larger than other Ranunculales species with sequenced genomes, such as Aq. coerulea (293.08 Mb), Eschscholzia californica (489.065 Mb), and Macleaya cordata (377.83 Mb) (Table S6).

Repeat Elements

Through a combination of approaches, we annotated 66.83% of the assembly as repetitive elements, among which long terminal repeats (LTRs) were the most abundant, occupying 40.62% of the genome assembly length; DNA elements and long interspersed nuclear elements occupied 5.0% and 3.0% of the genome, respectively (Table S7). The proliferation of LTR retrotransposons in K. uniflora was estimated to peak around 2.7 mya (Figure 2A). Analyses of age distributions built from synonymous substitutions per synonymous site (Ks) indicated that K. uniflora has undergone one recent whole-genome duplication (WGD) event, which occurred after its divergence from C. agrestis (Figure 2B). The inferred WGD event in the K. uniflora genome was further supported by dot-plot analysis of representative scaffolds, in which numerous paralogs derived from this event were identified (Figure 2C).

Figure 2

Evolutionary Analyses of the K. uniflora Genome

(A) Insertion time distribution of LTR retrotransposons. (B) Distribution of synonymous substitution rates (Ks) for pairs of paralogs/orthologs in/between K. uniflora and C. agrestis. (C) Dot plots of paralogs identified across contigs in the K. uniflora genome. (D) Venn diagram showing unique and shared gene families between genomes of K. uniflora and four other Ranunculales species.

Phylogenetic Tree Construction and Estimation of Divergence Times

Applying OrthoFinder (Emms and Kelly, 2015) to eight whole-genome and four transcriptome sequences including monocots, basal eudicots, and core eudicots, we identified a total of 18,742 orthogroups, among which 6,883 were shared by Kingdonia and four other Ranunculales species (Figure 2D). Among these orthogroups, 497 were identified as putative single-copy gene families. To further investigate the phylogenetic relationships within Ranunculales, we conducted both concatenated and coalescent analyses using the sequences of 497 single-copy genes in 12 species. The topologies from two analyses were identical, confirming the sister relationship between K. uniflora and C. agrestis and resolving Circaeasteraceae as sister to the clade formed by Ranunculaceae and Berberidaceae (Figure 3); Papaveraceae + Eupteleaceae was placed as the earliest-diverging clade (Figure 3). K. uniflora and C. agrestis were estimated to have diverged ∼51.8 mya in our analyses using MCMCtree with two calibration points (Figure 3). In addition, the phylogenetic analysis with an expanded group of taxa, which correspond to a larger taxonomic sampling but fewer loci indicated a similar phylogeny of Ranunculales, except the placement of Eupteleaceae (Figure S6).

Figure 3

Dated Phylogeny for 12 Plant Species with Oryza as an Outgroup

A timescale is shown at the bottom, and red points in some nodes indicate fossil calibration points.

Gene Family Overrepresentation and Underrepresentation

Comparisons of the genomes among 12 species identified a total of 111 gene families that are significantly (p < 0.01) overrepresented in K. uniflora and 22 gene families that are significantly underrepresented (Table S8). The results from Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) annotations showed that overrepresented gene families were considerably enriched in DNA repair pathways, such as homologous recombination, mismatch repair, DNA replication and nucleotide excision repair, whereas gene families showing significant underrepresentation in the K. uniflora genome were found to be involved in pathways related to stress or pest responses, such as the phenylpropanoid biosynthesis and secondary metabolites biosynthesis (Table 2).

Table 2

Functional Annotation of the Significantly Overrepresented and Underrepresented Gene Families in Kingdonia uniflora

Gene Families	KEGG Terms	Input No.	Background No.	p Value	Corrected p Value
Overrepresented gene families	Glycosphingolipid biosynthesis - globo series	9	9	5.28 × 10−7	1.59 × 10−5
	Homologous recombination	17	56	2.87 × 10−6	6.74 × 10−5
	Mismatch repair	12	39	7.31 × 10−5	0.000968079
	Sphingolipid metabolism	9	26	0.000282866	0.003082136
	DNA replication	12	50	0.000516123	0.005446878
	Nucleotide excision repair	14	69	0.000789672	0.007779277
	Peroxisome	16	87	0.000878589	0.008505414
	Galactose metabolism	12	55	0.001064483	0.010034842
	Plant hormone signal transduction	33	271	0.002344688	0.020456082
Underrepresented gene families	Cyanoamino acid metabolism	5	60	2.06 × 10−9	1.57 × 10−7
	Phenylpropanoid biosynthesis	5	157	2.06 × 10−7	6.24 × 10−6
	Starch and sucrose metabolism	5	202	6.95 × 10−7	1.96 × 10−5
	Biosynthesis of secondary metabolites	6	1,076	0.00020813	0.001073645
	Metabolic pathways	6	1,910	0.004094269	0.015173186

Dispensability of Plastid ndh Genes

To detect whether the ndh genes/segments lost from K. uniflora plastome were transferred to the nuclear genome, we conducted a BLASTN search using 11 intact ndh sequences extracted from C. agrestis as the query, using the assembled K. uniflora genome sequences as target. The result showed that no intact sequences for ndh plastid genes were discovered with the exception of ndhE and ndhJ, indicating functional copies of these genes likely have been lost (Table S9 and Figure 4).

Figure 4

Length Comparison of ndh Genes between K. uniflora and C. agrestis

Discussion

Species that live in stable habitats face less stress, which can cause lack or loss of ability to respond to environmental changes. In the present study, we de novo assembled the genome of K. uniflora, an ancient relictual species exclusively found in China. Given that K. uniflora has a larger genome than many of its close relatives, we hypothesize that the proliferation of LTR retrotransposons and the WGD event together are likely responsible for the increased genome size of K. uniflora (Michael, 2014). Several studies (Tenaillon et al., 2010, Michael, 2014) have suggested that the proliferation of TEs and specifically long terminal repeat retrotransposons (LTRs) in genomes is the primary driver of genome size differences in plants. This is because LTRs are expressed as RNA and reverse-transcribed into a new DNA element that can be inserted in every replication cycle (Wicker et al., 2007). A comparative study using high-quality genomes detected a correlation between intact LTR retrotransposons and genome size (El Baidouri et al., 2013). Abundant LTR retrotransposons were detected in the genome of K. uniflora, whereas considerably fewer LTR retrotransposons were identified in other members of the order with sequenced genomes (Figure 1). To investigate the evolutionary dynamics of the LTR retrotransposons, we estimated their insertion dates. The results indicate that the proliferation of LTR retrotransposons in K. uniflora was likely triggered around one rapid uplift of the Hengduan Mountain region, occurring between the late Miocene and late Pliocene (Kirby et al., 2002, Clark et al., 2005, Sun et al., 2011, Wang et al., 2012, Wang et al., 2014, Meng et al., 2016, Xing and Ree, 2017). WGDs have been shown to pervade the evolutionary history of angiosperms (Landis et al., 2018), and K. uniflora is no different (Figures 2B and 2C). Therefore, the relatively larger genome size of K. uniflora compared with close relatives might be promoted by both LTRs proliferation and WGD events.

Based on phylogenetic inference and estimated divergence times, we speculate that the speciation of K. uniflora was promoted by the Himalayan orogeny. Previous studies based on DNA commonly recognize K. uniflora being closely related to C. agrestis (Kim et al., 2004, Wang et al., 2009, Sun et al., 2017). Our phylogenetic analyses confirmed the sister relationship between K. uniflora and C. agrestis. These species were resolved as sister to the clade formed by Ranunculaceae and Berberidaceae, a result incongruent with previous placements of Circaeasteraceae and Lardizabalaceae as sister pairs (Kim et al., 2004, Wang et al., 2009, Sun et al., 2017). The divergence estimation between K. uniflora and C. agrestis (∼51.8 mya) is consistent with a previous estimate of ∼52 mya (Ruiz-Sanchez et al., 2012) and also coincides with the timing of the first stage of Himalayan orogeny (Rowley, 1996, Huang et al., 2015). Similar to K. uniflora, C. agrestis has a narrow distribution confined to China and the Himalayas. However, the distribution range of Circaeaster is relatively larger; in places where Kingdonia occur Circaeaster can always be found, whereas in most regions that Circaeaster is distributed Kingdonia is absent. Previous studies have shown that orogeny could create conditions favoring speciation of resident lineages (e.g., Hoorn et al., 2013, Wen et al., 2014, Favre et al., 2015). Hence, we hypothesize that the divergence between K. uniflora and C. agrestis was likely driven by the rapid uplift of contemporary Himalayan orogeny.

An asexual reproductive system and overrepresentation of DNA repair genes together reduce genetic diversity of K. uniflora. Levels of genetic diversity are often associated with reproductive strategies of species (Otálora et al., 2013). Colonization, population persistence, and extinction probabilities are all influenced by the reproductive systems of species (Stephenson et al., 2000, Barbará et al., 2009, Wornik and Grube, 2010, Beatty and Provan, 2011, Otálora et al., 2013). K. uniflora primarily reproduces asexually via rhizomes. This reproductive mode could produce identical individuals rapidly but lacks recombination and the possibility to create genetic variation in offspring, reducing the opportunities for adaptive evolution (Eckert, 2001, Castonguay and Angers, 2012). In addition, without segregation and recombination, obligate asexual multiplication may push a species into extinction owing to the steady accumulation of deleterious mutations (Thomas et al., 2016). During the long-term evolutionary history of K. uniflora, deleterious mutation accumulation cannot be ruled out; if something goes awry, such as the occurrence of a fatal mutation, whole clusters of clones can be wiped out. The integrated DNA-repair mechanism overrepresented in K. uniflora would allow for a reduction in the accumulation of deleterious mutations. This DNA-repair system might also reduce genetic diversity produced by mutations. Having high genetic diversity is important for plants to respond to environmental changes. We thus speculate that the extremely narrow distribution range of K. uniflora is associated with low genetic diversity, which restricted suitable environments to simplex, equable habitats rather than multiple, divergent habitats. Specifically, K. uniflora can only live in high elevations with minor human disturbances, while being characterized by perennial cold temperatures of below 0°C. The extreme temperature is likely to cause DNA, RNA, and protein damage. The overrepresentation of gene families for DNA repair might also be one kind of protection from extremely low temperature.

Underrepresentation of genes associated with stress response in K. uniflora leads to degeneration of adaptive ability to environmental changes. Phenylpropanoids are believed to contribute to all aspects of plant responses toward biotic and abiotic stimuli (Vogt, 2010). As concluded by La Camera et al. (2004), when plants suffer environmental stress, pests, or diseases, phenylpropanoids could evoke relevant response mechanisms to protect the plant from damage. Similarly, secondary metabolites also play a role in plant defense against environmental stresses, pests, and diseases (Bennett and Wallsgrove, 1994). Given the relatively undisturbed ecological environment of K. uniflora and lack of habitat stress during growth, adaptation to such conditions resulted in the functional degeneration of stress response systems. In addition, phenylpropanoids can also promote invasion of new habitats (Bais et al., 2003, Vogt, 2010). Our results indicate the underrepresentation of gene families involved in phenylpropanoid biosynthesis, which might be one reason causing low ability of K. uniflora to invade new habitats.

We conclude long-term living in highly equable habitats lead to the underrepresentation of stress response genes, which finally resulted in loss of ability to adapt to changing environments. The asexual reproductive strategy promoted overrepresentation of DNA repair genes, which reduced genetic diversity associated with adaptive capacity to environmental changes. Hence, both the underrepresentation of stress response genes and overrepresentation of DNA repair genes are responsible for the low adaptive ability of K. uniflora.

Equable habitats likely promoted the dispensability of most ndh genes in K. uniflora. Most angiosperms contain 11 plastid ndh genes, whereas all ndh genes, except for ndhE and ndhJ, were found to be either pseudogenized (ΨndhA, ΨndhB, ΨndhD, ΨndhH, and ΨndhK) or absent (ndhC, ndhF, ndhI, and ndhG) in the K. uniflora plastome (Sun et al., 2017). All 11 plastid ndh genes are intact in C. agrestis (Sun et al., 2017), indicating that the loss of ndh genes from K. uniflora occurred after the split between K. uniflora and C. agrestis, suggesting that within 52 million years most of the plastid ndh genes were lost from K. uniflora not only in the plastome but also in the nuclear genome. Among land plants, the plastid ndh loci have also been found absent in non-photosynthetic plants, epiphytes, Gnetales, conifers, and Erodium (Geraniaceae) (Kim et al., 2015, Lin et al., 2017, Nie et al., 2017). Evidence suggests that the thylakoid NADPH complex could optimize photosynthesis for plants under environmental stresses, while being found dispensable for plant growth under optimal growth conditions (Martín and Sabater, 2010). K. uniflora is extremely selective in habitat preference (Li et al., 2001) and is known as the indicator for natural ecological environment without disturbance. We hence speculate that the current habitats of K. uniflora might have promoted the dispensability of the plastid ndh genes. Additionally, within plants, NADPH supplies hydrogen for many anabolic processes (Antal et al., 2015). The underrepresentation of gene families related to metabolic pathways, as detected from our CAFÉ-based analyses, is likely related with the nonfunction state of plastid ndh genes.

Changing climate, shrinking habitats, and low adaptive ability to environmental changes together contributed to the extremely narrow distribution of K. uniflora. The overrepresentation of gene families involved in DNA repair could help reduce the accumulation of deleterious mutations during asexual reproduction, which is the dominant mode of reproduction in K. uniflora, while reducing genetic diversity, which is important in responding to environment fluctuations. The underrepresentation of gene families in charge of stress response and nonfunction of plastid ndh genes could be due to the adaptive degeneration caused by long-term adaptation to living in relatively stress-free environments. Considering the long evolutionary history of K. uniflora, and the fossil records from C. agrestis in the mid-Albian of Virginia, USA (Crane et al., 1994, Drinnan et al., 1994), we speculate it should have been widespread around the world. Changing climate, shrinking habitats, asexual reproduction, and adaptive degeneration caused by relying on easeful environment together led it to a current status of endangerment.

Limitations of the Study

We reported the draft genome of the ancient and currently endangered K. uniflora. Although we revealed the evolutionary history and potential endangered mechanism of K. uniflora through primary analyses, population genomic analyses are required to estimate genetic diversity, reveal genetic structure and its causes, and identify loci associated with environmental adaption.

Resource Availability

Lead Contact

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Yanxia Sun (sunyanxia@wbgcas.cn).

Materials Availability

This study did not generate new reagents.

Data and Code Availability

The accession number for the genome assembly and raw reads reported in this paper is GenBank: PRJNA587615 and PRJNA611722, respectively.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.