Plant materials and growth conditions

Ginkgo biloba L. (Ginkgoaceae) seedlings were harvested from Xi’an botanical garden (E, 108°93′, N, 34°17′, Shaanxi Province, Northwest China) and grown in a greenhouse under long-day conditions (16-h light/8-h dark cycle) at 28 ± 2 °C. Leaves were harvested from 8-week-old plants, and frozen in liquid nitrogen.

DNA isolation and PCR

The DNA was isolated using an improved CTAB protocol. Plant leaves (0.1 g) were ground into powder in liquid nitrogen. Then, 0.6 mL CTAB extraction buffer was added and the lysate was incubated at 65 °C for 30 min. The DNA was purified by adding an equal volume of a mixture of chloroform: isoamyl alcohol (24:1) followed by centrifugation at 8000 × g for 10 min at 4 °C. The supernatant was added to 2/3 volume of isopropanol and then subjected to centrifugation at 8000 × g. The precipitate was washed twice with 75% ethanol and then dissolved in 300 μL sterile water. NaAc (1/10 volume of 3 M, pH 5.2) and two volumes of ethanol were added to the tube followed by a 10-min incubation at −20 °C. The tube was centrifuged at 8000 × g for 5 min and the pellet was then washed twice with 75% ethanol and re-dissolved in 20 μL sterile water.

The primers of 82 G. biloba transcripts were designed based on the G. biloba chloroplast complete genome [AB684440], and the primer sequences are listed in (Additional file 1: Table S1). The PCRs were performed as follows: 95 °C for 3 min, 94 °C denaturing for 30 s, 53–60 °C annealing for 30 s, and an elongating time between 30 s and 1.5 min at 72 °C based on the DNA length (1 min per 1 kb). The PCR amplification products were electrophoresed on a 1% agarose gel and purified with E.Z.N.A^TM Gel Extraction Kit (OMEGA Bio-Tek, USA). The direct sequencing of cDNAs derived from these transcripts and of the corresponding genomic DNA (gDNA) was carried out by Sangon Biological Engineering Technology & Services (Shanghai, China).

RNA isolation and RT-PCR

The total RNA was extracted using E.Z.N.A^TM Plant RNA Kit according to the manufacturer’s protocol. The tissue was disrupted and homogenized as above, and the gDNA was preliminarily eliminated with a gDNA filter. The flow-through at the very last step was mixed with the membrane-binding solution and then loaded into the HiBind RNA Mini column. Finally, RNA was washed with RWC buffer and RNA wash buffer to remove protein, polysaccharide and salt contamination. The total RNA was treated with DNaseI to remove gDNA contamination. The cDNA was synthesized according to the PrimeScript RT Reagent Kit protocol (TaKaRa, Dalian, China).

RNA editing site identification

Direct sequencing was used in this paper. The PCR products were purified and sequenced at least three times. The editing sites were detected by aligning the DNA and cDNA sequences one by one using the EMBL-EBI ClustalW (http://www.ebi.ac.uk/Tools/msa/clustalo/). The sequences were analyzed using SeqMan of the Lasergene software package (https://www.dnastar.com/t-seqmanpro.aspx). According to Mower and Palmer [32], T and C appeared at the same site and clearly above the background, indicating partially edited sites.

Analysis of the protein structures, and their composition before and after editing

MegAlign of the Lasergene package was used to analyze protein similarities. The N-terminal signal peptide prediction was carried out by SignalP (http://www.cbs.dtu.dk/services/SignalP), and SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html) was employed to analyze the changes in the secondary structure. TMHMM (http://www.cbs.dtu.dk/services/TMHMM/) was used to predict alterations in the transmembrane region.

Evolution analysis of RNA editing genes

For the RNA editing evolutionary analysis, the ndh, pet, psa and psb gene families from 12 species were selected, and then a z-test was applied to detect selection constraints using Mega 5.1 software. The non-synonymous–synonymous (dN–dS) substitution rate analysis was also conducted for each gene according to the Goldman and Yang (GY-94) method in Hyphy, which estimates dS and dN substitution rates through a codon-based model [33–35]. Parameters were set as follows: Test hypothesis mode was set as Neutrality. Nei-Gojobori method was chosen in the substitution mode. In general, a dN value lower than dS (dN < dS) suggests negative selection, i.e. nonsilent substitutions have been purged by natural selection, whereas the inverse scenario (dN > dS) implies positive selection, i.e. advantageous mutations have accumulated during the course of evolution.

The homologues gene sequences and editing sites used in this paper

The 12 species used for the sequence alignments are listed as follows: A. belladonna [NC_004561.1]; S. lycopersicum [AM087200]; C. sativus [AJ970307]; A. formosae [NC_004543.1]; G. hirsutum [DQ345959.1]; A. thaliana [NC_000932.1]; C. taitungensis [NC_009618]; A. capillus-veneris [AY178864.1]; T. aestivum [AB042240.3]; N. tabacum [Z00044.2]; Z. mays [NC_001666.2]; G. biloba [AB684440.1]. Most of the editing site information was acquired from GenBank and RNA (http://dna.kdna.ucla.edu/rna/index.aspx) databases. Some editing sites were found in the literature.

G. biloba chloroplast transcripts undergo several editing events

We further analyzed the RNA editing frequencies of different gene groups in the chloroplast genome of G. biloba. The results showed that ndh genes exhibited the most editing cases, which were nearly 36% of the total editing sites, while the number of cases was not more than 10% in other genes (Fig. 1a).

To exclude interference by the gene length on the editing events, the number of corresponding editing sites was divided by the length of each gene group. ndh and clpP exhibited the same, and the highest, editing frequency, up to 8.5‰. Interestingly, rbcL had an almost undetectable editing frequency (Fig. 1b). These data suggested that ndh genes are more likely to be edited than other genes at the mRNA level.

The characteristics of the RNA editing sites in the G. biloba chloroplast genome

To gain further insights into the characteristics of the 255 RNA editing sites in the G. biloba chloroplast genome, we analyzed different types of editing codon positions. There were 63, 174 and 14 editing sites occurring at the first, second and third codon positions, respectively (Fig. 2). Editing sites occurred in second or third positions in one codon of the ycf3, psbB, rps14 and ndhD transcripts (Fig. 2, Additional file 2: Table S2). For the editing sites distributed in the first codon positions, there are 37 sites in front of purine (adenine or guanine at the second codon position), which makes up ~59% of the editing occurring in the first codon positions. In the second codon position, editing occur in a U_A context (50), followed by U_G (27), C_A (21), U_U (16), C_G (15) and U_C (14) context (the numbers in parentheses refer to the number of RNA editing sites in which editing occurred at the second position in a codon) (Fig. 3).

Most RNA editing sites exist in the protein-coding regions and often cause corresponding amino acid alterations. In addition to 16 silent editing sites, there were 239 sites that resulted in corresponding codon changes in G. biloba. Among them, 132 editing sites switched amino acids from hydrophilic to hydrophobic, and more than 60.5% of the editing events were serine to leucine, followed by serine to phenylalanine (24.2%) and threonine to isoleucine (8.3%). The amino acids maintained their hydrophobic properties at 80 editing sites, and the highest rate occurred in proline to leucine (60.0%), followed by histidine to tyrosine (20.0%) and leucine to phenylalanine (1.3%). Only 13 and 7 editing sites caused amino acids to change from hydrophobic to hydrophilic and to maintain their hydrophilicity, respectively (Fig. 4).

RNA editing events in G. biloba chloroplast genes may alter protein structures

In our attempt to understand whether RNA editing affects protein structure, we predicted the secondary structures of 82 proteins before and after editing using bioinformatics software. The results showed that many editing events might change the secondary structures of the corresponding protein. Most editing sites form a new α-helix structure in up- or down-stream regions around the editing codon (Additional file 3: Figure S1). A new cleavage site in the signal peptide within the 18th and 19th codon positions was created in ndhD-57 (Additional file 3: Figure S2). Five new transmembrane regions appeared in ndhD, ndhE, ndhF, psbB and psbN, respectively, after the corresponding codons were edited (Fig. 5a-e). In addition, a transmembrane region disappeared in petB when the amino acid at the 212 codon position changed from proline to serine due to editing (Fig. 5f).

Comparison of RNA editing sites in different species

Evolutionary pattern of RNA editing events in chloroplasts

To investigate the evolutionary tendency of RNA editing, the dN-dS values of the RNA editing sites in four photosynthesis-related gene families were calculated using the Z-test of selection in MEGA5.1 Beta software. The dN-dS values of most ndh and some psb genes were greater than zero (Fig. 6a and b), indicating that these editing sites may have undergone positive selection. The dN-dS values of most of the psa, a few psb and the pet genes studied, except for petB, were equal to zero (Fig. 6b, c and d), suggesting that editing sites in these genes may undergo neutral selection. However, we noticed the tendency of the dN-dS values to trend to zero in most psa, psb and pet genes was faster than in the ndh gene group (Additional file 3: Figure S3). This occurred because C-to-T point mutations at the genome sites in most of the psa, psb and pet gene families caused the editing sites to disappear. Moreover, the C-to-U editing at the mRNA level and the reverse mutations at the genome level can both increase codon conservation. For example, petA-329, psaA-725 and psbF-77 were edited in G. biloba, but they underwent a reverse mutation to T at the DNA level in A. thaliana, T. aestivum and Z. mays, causing an increase in the corresponding codon conservation in most of the species (Additional file 3: Figure S4). The results contradicted those of what is commonly referred to as neutral selection, in which mutations are neither beneficial nor detrimental to the ability of an organism to survive and reproduce [36]. In fact, the conservation of amino acids is restored in most of these gene classes due to C-to-T point mutations at the genome level. Thus, C-to-U edits at the mRNA level are unnecessary and even waste energy. As a result, editing sites in these essential genes gradually disappeared during evolution. The evolutionary tendencies of RNA editing in these gene classes acts more like a purifying selection, so, we termed this kind of evolution as ‘similar purifying selection’, in which dN–dS is equal to zero but purifying selection actually occurred to retain codon conservation.

Abundant RNA editing events are retained in the chloroplast genome of G. biloba

Except for the marchantiid subclass of liverworts, RNA editing has been observed in the chloroplasts of all of the investigated terrestrial plants. The number of C-to-U RNA editing sites in chloroplasts was variable among plants, ranging from 0 in Volvox globator to more than 900 in A. formosae. Over 300 chloroplast editing sites were known in early branching land plants, such as Anthoceros and Adiantum, but fewer sites were edited in angiosperm chloroplast RNAs (Fig. 7). In this paper, we reported that the chloroplast protein-coding transcripts of G. biloba contain 255 editing sites, which is by far the highest number of editing sites in a seed plant. A model for the evolution of editing in plant organelles proposed that RNA editing was of monophyletic origin, had a common ancestor with many editing sites during seed plant evolution, and that many of the original editing sites, particularly in seed plants, had been subsequently lost [37]. G. biloba is one of the oldest seed plants and appeared in the Early Jurassic period, in which the CO₂ concentration in the atmosphere may have reached high levels, accelerating climate warming [38]. All of these changes may cause G. biloba to acquire many mutations at the DNA level and RNA editing recovered the equivalent genetic information. In addition, comparisons of editing events among three gymnospermaes, G. biloba, Pinus and Cycas, showed that a large number of unique editing sites of G. biloba had been lost in Cycas and Pinus (Fig. 8). G. biloba may maintain a more ancestral version of the chloroplast genome than Cycas and Pinus. Moreover, G. biloba shares 11 and 3 editing sites with Cycas and Pinus, respectively, and three common editing sites, atpF-370, petB-634 and psbE-214, are shared among the three species (Fig. 8). This indicated that the evolutionary conservation of RNA editing is essential for only a few plastid editing sites, which is a common phenomenon among angiosperms and has been verified in many cases [39].

RNA editing might change the structures and functions of some proteins in G. biloba

RNA editing, especially at the second codon position, can alter the encoding amino acid and change the protein primary, secondary or tertiary structures, which might be necessary for the protein function. We analyzed the secondary structures of 82 transcripts in G. biloba before and after editing using bioinformatics methods. One editing site changed the signal peptide, eight editing sites could create five new transmembrane regions, and one RNA editing event occurred in petB, which caused an existing transmembrane region to disappear. All of the newly created signal peptides and transmembrane regions might play important roles in the localization or formation of the proper spatial structures of the proteins, especially for membrane proteins. Until now, a great deal of experimental evidence supported the view that most of the unedited proteins had lower functional levels than the edited proteins. In peas, unedited acetyl-coA carboxylase carboxyl transferase β is not functional and cannot catalyze the synthesis of fatty acids [40]. In maize chloroplast rpl2, the AUG initiation codon generated by a C-U editing of ACG is essential to seed development [41]. In Arabidopsis, an editing defect at atp1-C1178 has a strong impact on the assembly of the ATP synthase [42]. Of the 255 editing sites in G. biloba, two types mainly cause the conversion of amino acids from serine to leucine or phenylalanine and proline to leucine. The former might increase the hydrophobicity of the corresponding peptide and the latter did not change the peptide hydrophobicity, but it could recover the normal curl of the secondary structure or remove misfolding because proline is a helix-breaker.

In addition, the chloroplast genome of G. biloba has many partial editing sites. Among the 255 editing sites, 73 partial editing sites were detected. Tseng et al. found that partial editing may regulate plastid gene expression by using a different editing frequency in the non-photosynthetic tissues of Arabidopsis [43]. Karcher et al. have found that ndhB has different editing profiles in photosynthetic and non-photosynthetic organs in Z. mays [44]. Thus, many partial editing sites in G. biloba might be associated with different tissues and developmental periods. Further experiments are needed.

RNA editing may undergo diverse evolutionary patterns in different photosynthetic genes

The RNA editing phenomenon may be a relic of the ancient RNA world that is involved in primordial error correction, such as repairing UV damage or other uncertain factors at the transcriptional level. As a result, RNA editing appears in almost all land plants, except Marchantia polymorpha of the marchantiid subclass of liverworts [45]. With evolution, the number of RNA editing sites gradually decreases from the lower to higher plants (Table 2). To understand the evolution of plastid editing sites, we introduced the dN-dS method to predict the evolutionary mode of RNA editing sites. Comparisons of editing sites in the ndh, psa, psb and pet genes in 12 plant species revealed that the dN-dS values of psa, most of psb and the pet gene groups were nearly equal to zero (Fig. 6). Additionally, the tendency of the dN-dS values to trend to zero in most of the psa, psb and pet genes was faster than in the ndh gene group (Additional file 3: Figure S3). Thus, these genes may undergo similar purifying selection. Most of the genes had an important role in photosynthesis. For instance, the targeted inactivation of psaI affects the association of psaL with the photosystem I core. Namely, the absence of psaI indirectly leads to a defect in photosystem I function [46]. Varotto et al. disrupted the A. thaliana photosystem I gene, psaE, and observed several defective phenotypes, including a significantly increased light sensitivity and a decreased growth rate of ~50% under normal conditions [47]. Additionally, losing PsbJ in tobacco causes the photosynthetic performance to be drastically reduced, as well as an extreme hypersensitivity to light [48]. Salar Torabi et al. also reported mutants in psbN-F and psbN-R of N. tabacum were extremely light sensitive and failed to recover from photo inhibition [49]. Fiebig et al. proposed that essential genes cannot tolerate frequent T to C mutations at the DNA level [50]. For the essential genes, such as psa, psb and pet, most of them have abundant editing sites in ancient species, but many editing sites disappeared during plant evolution due to reverse mutations at the DNA level that restored codons to conserved amino acid residues. Those editing sites were probably essential for the structure and/or function of the encoded protein.

The plastid ndh genes encode a thylakoid Ndh complex that purportedly acts as an electron feeding valve to adjust the redox level of the cyclic photosynthetic electron transporters [51]. By far the highest number of plastid editing sites in flowering plants was found in the ndh group of genes [52]. In our research, ndh genes also possessed the most editing sites and had the highest editing frequency. The ndh gene groups might be unessential for plants growing under normal conditions. Burrows et al. hypothesized that the ndh complex was dispensable for N. tabacum growth under optimal growth conditions [53]. Ndh genes are absent in epiphytic plants [54] and are partially lost in Phalaenopsis, Aphrodite and Erodium [55]. In P. thunbergii, most of the ndh genes are pseudogenes. Thus, we speculated that the RNA editing sites of the ndh genes might be randomly lost and that the loss rate was slow. Therefore, ndh genes could keep more editing sites than other gene groups in modern plants. For the ndh gene group, we found that RNA editing in ndhD, ndhF and ndhG might create obvious structural changes, which created a new transmembrane region or caused an existing transmembrane region to disappear after editing (Fig. 5). To a certain extent, its occurrence implies that editing in those genes has biological significance. In Arabidopsis, the editing deficiency in ndhF was associated with a delayed greening phenotype [56]. The decline of the editing efficiency in ndhB and ndhD affected the flow of cyclic electrons and enhanced disease resistance [57]. Although the products of the majority of ndh genes were unnecessary under standard growth conditions, editing was probably most important for the proper function of the NDH protein complex under stress conditions [58, 59]. Due to the RNA editing, ndh genes might improve photosynthesis and stress tolerance under harmful conditions, and they may display positive selection during evolution. Some psb and pet genes, which have a dN–dS bias greater than zero, such as psbE, psbF, psbH, psbJ, psbL, psbT, petB and petL may have similar evolutionary mechanisms. Thus, RNA editing may be a post-transcriptional regulatory process of ancient genes, as well as part of an evolutionary model with diverse evolutionary directions [60]. We speculated that the editing sites in each gene may undergo diverse evolutionary paths depending on whether the edited codon was important or not for protein executive functions.