1 INTRODUCTION

The chloroplast NADH dehydrogenase-like (NDH) complex is one of the thylakoid membrane protein complexes that involved in photosynthetic electron transport chains. It catalyzes the transfer of electrons from ferredoxin (Fd) to plastoquinone to generate trans-thylakoid protons to drive the production of adenosine triphosphate (ATP) (Burrows et al., 1998; Shikanai, 2007). Cyanobacterial NDH-1 can be divided into two major parts that contain the membrane segments (NdhA-NdhG) that participate in proton translocation and the peripheral segments (NdhH-NdhK) that carry redox centers, Fe-S clusters, and flavin mononucleotide. It forms an L-shaped structure, which is conserved in both photosynthetic and respiratory NDH complexes (Matsubayashi et al., 1987; Suorsa et al., 2009; Efremov et al., 2010; Pan et al., 2020). Additionally, some new subunits that are required to stabilize the NDH complex were found in higher plants. Thus, the complete structure of the chloroplast NDH complex is composed of five subcomplexes, including subcomplex M (NdhA-NdhG), subcomplex A (NdhH-NdhO), subcomplex B (PnsB1-PnsB5), subcomplex L (PnsL1-PnsL5), and subcomplex EDB (Electron donor binding) (NdhS-NdhV) (Ifuku et al., 2011).

Since the defective ndhB was verified to influence photosystem I cyclic electron transport (PSI-CEF) in Synechocystis sp. PCC 6803 (Ogawa, 1991), chloroplast NDH was proven to be involved in PSI-CEF (Burrows et al., 1998; Shikanai et al., 1998; Shikanai, 2007). Although NDH-dependent PSI-CEF subtly contributes to the total PSI-CEF under optimal growth conditions, the chloroplast NDH complex plays essential roles in the regulation of plant tolerance to abiotic stress. For example, the chloroplast NDH complex functions to alleviate the oxidative stress caused by heat and high light stress (Li et al., 2004; Wang et al., 2006; Essemine et al., 2016; Yang et al., 2017; Bertelli and Unsworth, 2018; Tan et al., 2020a). Moreover, NDH-dependent PSI-CEF is necessary for the normal growth of many plants (Shikanai, 2016). There are 11 plastid-encoded Ndh subunits (NdhA-NdhK) that are homologous to their counterparts in the mitochondrial NADH dehydrogenase (Ohyama et al., 1986; Shinozaki et al., 1986). The chloroplast NDH complex is regarded to be coupled with the plastid terminal plastoquinone oxidase to mediate respiratory electron transport, which is referred to as chlororespiration (Peltier and Cournac, 2002).

To enhance the turnover of electron flux via Fd and effectively resist stress, the NDH complex has attained dozens of novel subunits during the evolution of land plants (Shikanai, 2016). Subcomplex M and several subunits in Subcomplex A (NdhH-NdhK), which forms the skeleton of NDH complex, are highly conserved in Escherichia coli, cyanobacteria, and land plants (Efremov et al., 2010). Although plastid-encoded ndh (ndhA-ndhK) are found in some ancestral algae, e.g., Nephroselmis olivacea or Mesostigma viride (Turmel et al., 1999; Lemieux et al., 2000), they are absent in most algal species, including the red and green algae (Peltier and Cournac, 2002). Subcomplex B, which participates in the maintenance of the complex stability, emerged during the early stage of evolution of land plants (Ruhlman et al., 2015). Marchantia polymorpha, which is considered to be the earliest case of divergence from the land plant lineage, merely contains a portion of Subcomplex L (Ueda et al., 2012), indicating that this subcomplex is unique to terrestrial higher plants. NdhS and NdhV in Subcomplex EDB form the Fd-binding site and are conserved in cyanobacteria and terrestrial higher plants, while other subunits in Subcomplex EDB are specific to terrestrial higher plants (Battchikova et al., 2011; He et al., 2015; Gao et al., 2016). The integral chloroplast NDH complex is primarily found in terrestrial higher plants, although some species in the Orchidaceae, Pinaceae, Gnetaceae, and Geraniaceae have lost their chloroplast NDH complex during evolution (Braukmann et al., 2009; Weng et al., 2014; Kim and Chase, 2017; Lin et al., 2017). The reason for the drastic variation in evolution of chloroplast NDH complex still remains unclear.

Zostera marina (Alismatales: Zosteraceae) is one of the most productive types of seagrass and widely distributed in the northern Pacific Ocean and northern Atlantic Ocean (Olsen et al., 2016). As the main component of seagrass meadows, it has crucial ecological values in nutrient cycling, sediment stabilization, and the provision of habitats and food for other organisms (Waycott et al., 2009). Approximately 200 million years ago, it migrated to terrestrial conditions and returned to the sea environments approximately 140 million years ago (Les et al., 1997). During this process, it occurred numerous gene losses and gains to adapt to its marine life style (Olsen et al., 2016).

In our study, all 31 ndh genes of the entire NDH complex were identified in Z. marina, while they have rarely been detected in marine macrophytes even among the seagrasses in the order Alismatales (Iles et al., 2013). Correspondingly, our previous research suggested that Z. marina possesses a highly efficient NDH-dependent PSI-CEF and chlororespiration (Tan et al., 2020a, b). Although the structure, physiological function, and assembly of chloroplast NDH has been reported in other species (Peng et al., 2011; Shikanai, 2016), the specifically-evolved and highly efficient Z. marina NDH has yet to be investigated to our knowledge. The whole genome sequences that have been reported enabled us to conduct a systematic research of ndh genes in Z. marina (Olsen et al., 2016; Xing and Guo, 2018). We comprehensively analyzed the phylogeny, gene structures, patterns of expression, and regulation of ndh genes in Z. marina. This study provided clues to understand how the ndh coordinate with each other to drive the highly efficient NDH-dependent PSI-CEF in response to light exposure and chlororespiration.

2 MATERIAL AND METHOD 2.1 Plant material and treatment

Zostera marina with intact rhizome systems were collected from subtidal seagrass beds in Yandunjiao, Rongcheng (37°91′N, 120°73′E), Shandong Province, China. Sampling was performed in the late afternoon during consecutive days in May 2020. Samples were identified by experienced taxonomists on the basis of their morphology. The identifications were also confirmed by sequence similarity with the whole genome sequencing of Z. marina (Olsen et al., 2016). The collection of plant materials complied with institutional, national, and international guidelines. Samples were cultured in an aquarium with seawater that was continuously aerated and renewed daily. Before experimentation, the plants were pre-cultivated for 3 days under 15 ℃ with a photoperiod of 10 h꞉ 14 h (light꞉dark) in minimum saturation light intensity (100 μmol photons/(m²·s)).

The leaves obtained from dark-adapted overnight plants were exposed to 300 μmol photons/(m²·s), followed by recovery under darkness. For reverse transcription-PCR (RT-PCR) assays, the samples were collected after light exposure and recovery for 3 h, respectively. For quantitative real time-PCR (qRT-PCR) assays, the samples were collected after light exposure and recovery for 10 min, 30 min, 1 h, and 3 h, respectively. The leaves under dark-adapted overnight conditions were used as the control. Moreover, the roots, leaves, flowers, stems, and rhizomes were collected at the flowering stage for tissue-specific expression analyses using RT-PCR and qRT-PCR. There were three biological replicates for each sample.

2.2 Characterization of Ndh subunits and sequence analysis

BLASTp analyses were performed against the Z. marina database from the NCBI (https://www.ncbi.nlm.nih.gov/) and Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) databases using the Arabidopsis thaliana Ndh protein sequences. The reannotated and misannotated genes were submitted to GenBank and can be retrieved with ID 2387379 and 2387367. The number of amino acids, molecular weights, and theoretical isoelectric points of Ndh subunits were attained by submitting protein sequences to the ExPasy website (http://web.expasy.org/protparam/). The conserved domains of each ndh gene were confirmed using Pfam and SMART (http://smart.embl-heidelberg.de/). The exon-intron compositions of the ndh genes were constructed using DNAMAN software by comparing the cDNA with genomic sequences. The duplication modes of the gene pairs were identified at the Plant Duplicate Gene Database (PlantDGD, http://pdgd.njau.edu.cn:8080) (Qiao et al., 2019)

2.3 Phylogenetic analysis

Phylogenetic trees were generated from the amino acid alignments of both the entire concatenated complex (31 subunits) and the individual Ndh subunits using MEGA5.0 with the following settings: Poisson model, pairwise deletion, and 1 000 bootstrap replications (Tamura et al., 2011). Trees containing evolutionary model species in Monocots, Dicots, Amborellales, Lycophytes, Bryophytes, Charophyceae, Chlorophyta, and Cyanobacteria were analyzed, using E. coli as the outgroup. The protein sequences of each species were acquired from the Phytozome and NCBI databases. All of the sequences were aligned using ClustalW2 (Larkin et al., 2007). The Interactive Tree of Life (https://itol.embl.de/) was used to visualize the trees.

2.4 Structural characterization

The motifs of Ndh subunits were identified using the MEME (Multiple Em for Motif Elicitation) online program (http://meme.nbcr.net/meme/intro.html) with the following parameter sets: any number of repetitions, maximum number of 20 motifs, and the optimum motif widths of 6 to 200 amino acid residues (Bailey et al., 2009). Trans-membrane helices in the Ndh subunits were predicted using TMHMM v2.0 (Krogh et al., 2001). The multiple sequence alignments and structure analysis were conducted with ESPript 3 in six species, including Thermosynechococcus elongatus BP-1, Zea mays, Sorghum bicolor, A. thaliana, Populus trichocarpa, and Z. marina (Robert and Gouet, 2014).

2.5 Nucleic acid isolation and cDNA synthesis

Total DNA was extracted from 100-mg leaf tissue using the cetyltrimethylammonium bromide (CTAB) method (Chen et al., 2011). Total RNA was isolated using a FastPure Plant Total RNA Isolation Kit (Vazyme, Nanjing, China). The quality of RNA was examined by electrophoresis on a 1.0% agarose gel and quantified using NanoQuant (TECAN Group Ltd., Männedorf, Switzerland). After wiping off the residual DNA, the cDNA was synthesized with a HiScript^® Ⅱ 1st Strand cDNA Synthesis Kit (Vazyme) using 1 μg of total RNA.

2.6 PCR amplification and sequencing

Primers were designed for polymerase chain reaction (PCR) and RT-PCR amplification using the DNA sequences obtained from GenBank (Supplementary Table S1). DNA and cDNA templates were amplified using PrimeSTAR^® Max DNA Polymerase (TaKaRa, Kyoto, Japan) with the following procedure: 98 ℃ for 10 s, 35 cycles of 55 ℃ for 10 s, and 72 ℃ for 30 s. The amplification products were separated on 1.0% agarose gels and purified using a FastPure Gel DNA Extraction Mini Kit (Vazyme) before Sanger sequencing.

2.7 Analysis of splicing patterns

The RT-PCR products of ndh genes from five tissues and leaves exposed to different light periods were visualized using a Gel Doc XR+ system (Bio-Rad, Hercules, CA, USA). The specific electrophoretic bands corresponded to differently spliced products. The possible protein isoforms were predicted by DNAMAN.

2.8 RNA editing analysis of chloroplast-encoded ndh genes

To confirm the RNA editing sites in the Z. marina ndh genes, the sequenced cDNA and DNA from different light periods were aligned using Vector NTI. The editing sites for Spirodela polyrhiza (Wang et al., 2015), A. thaliana (Tillich et al., 2005), Z. mays (Maier et al., 1995), and Amborella trichopoda (Hein et al., 2016) were obtained from respective publications.

2.9 Analysis of gene expression

Quantitative real time-PCR assays of the ndh genes obtained from five tissues and different light periods were conducted on a Bio-Rad CFX96 Real Time PCR System with AceQ Universal SYBR qPCR Master Mix (Vazyme). The housekeeping gene gapdh from Z. marina was used as an internal control. The qRT-PCR was programmed as follows: 95 ℃ for 10 s, followed by 40 cycles of 56 ℃ for 10 s and 72 ℃ for 30 s. The qRT-PCR data were calculated using the 2^-ΔΔCt method. Sequences of the primers used in qRT-PCR are described in Supplementary Table S1. The tissue-specific expression of A. thaliana was downloaded as the publicly available RNA-Seq data in the CoNekT-Plants database (www.conekt.plant. tools). The RNA Seq of Z. marina under light stress was analyzed in our previous study (Tan et al., 2020b). The differentially expressed genes were selected according to P < 0.05 and fold change > 1.3 (for up-regulation) or fold change < 0.75 (for down-regulation). Heatmaps were constructed using the transformed log2 (TPM+1) or log2 (FPKM +1) values by TBtools software (Chen et al., 2020).

2.10 Analysis of cis-acting elements in the ndh gene promoters

To obtain the possible cis-acting elements in the promoter regions of the ndh genes, 1 500-bp genomic sequences upstream of the start code ATG were analyzed online using Plantcare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

2.11 Statistical analysis

The statistical analyses were performed using SPSS 22.0 (IBM, Inc, Armonk, NY, USA). All data were analyzed by a one-way analysis of variance (ANOVA) and Tukey's tests with P < 0.05 set to be statistically significant.

3 RESULT 3.1 Identification of the Ndh subunits

A total of 31 Ndh subunits homologous to their A. thaliana counterparts were obtained from Z. marina. Gene characteristics, including the length of the protein sequence, protein molecular weight (MW), isoelectric point (pI), subcellular localization, and domain type were analyzed (Fig. 1, Table 1). Duplication events NdhB and PnsB3 occurred. The annotations of these subunits were validated using the available transcriptome data of Z. marina. Gene Zosma266g00030.1 encodes the NdhV subunit and Zosma266g00040.1 encodes a hypothetical protein that are part of the same transcript (designated Zosma266g0003040.1). Following the confirmation by RT-PCR and sequencing, in silico translation revealed that Zosma266g0003040.1 translated a 191 amino acid protein. Furthermore, the misannotated coding sequence of ndhB was observed during the analysis of sequenced cDNA fragments in the section of RNA editing. The duplicated ndhB had a uniform coding sequence with a 509 amino acid protein. To explore the evolutionary clues of duplicated ndh gene pairs, five widely accepted duplication modes (whole genome, tandem, proximal, transposed, and dispersed duplication) were systematically scanned. The gene pairs ndhB-1/ndhB-2 and pnsb3-1/pnsb3-2 originated from the expansion of large inverted repeat regions (IRs) and transposed duplication, respectively.

3.2 Phylogenetic analysis of ndh genes

To investigate the evolutionary relationship among Viridiplantae NDH complexes, the amino acid sequences of 31 Ndh subunits were concatenated to construct a phylogenetic tree. As shown in Fig. 2, the phylogeny of NDH complex was almost consistent with the phylogeny of species. The NDH complex appeared to be increasing in complexity from algae to bryophytes to land higher plants with the exception of gymnosperms. Therefore, the most integrated NDH complex was primarily found in the terrestrial angiosperms with an obvious diversification between monocots and dicots. Unlike other marine species, Z. marina had a complete NDH complex, which clustered with terrestrial angiosperms, indicating a special evolutionary status.

The phylogenetic tree of the entire NDH complex was also compared with those constructed by individual Ndh subunits. The phylogenies of most Ndh subunits, including NdhA, NdhE, NdhF, NdhG, NdhI, NdhJ, NdhN, NdhT, NdhU, PnsL1, PnsL2, PnsL3, and PnsL4, were consistent with the concatenated phylogeny (Supplementary Fig.S1a), indicating that these individual subunits co-evolved with the NDH complex. NdhB, PnsB3, NdhM, NdhO, and NdhS of Z. marina were at the basal branches of angiosperms (Supplementary Fig.S1b), indicating their earlier evolution status.

3.3 Motif features of the Ndh subunits

The motifs of each Ndh subunit in the conserved domain are shown in Supplementary Fig.S1 and Supplementary Table S2, respectively. It revealed that closely related species usually shared a similar motif composition. The motifs in conserved domain were widely distributed, while the motifs in the N- and C-terminal regions varied, particularly for nuclear-encoded Ndh subunits. Some are clearly unique to the Z. marina Ndh subunits, i.e., motif 20 in NdhG, motif 13 in NdhH, and motif 14 in NdhI, whereas some conserved motifs were absent, i.e., motif 11 in NdhD and motifs 7, 14, and 17 in NdhF (Fig. 3). It is worth noting that the antiporter-like subunits NdhD and NdhF, which are involved in the transmission of the protons, lost the motifs that are either located in the functional domain of the subunits or predicted to be a transmembrane helix.

Multiple sequences of NdhD and NdhF were aligned in six representative species of angiosperms. As shown in Fig. 4, the sequences in Proton_antipo_M domain were highly conserved. NdhD harbors 14 conserved TM (transmembrane) helices, and NdhF harbors 16 conserved TM helices. Unlike other species, Z. marina lost TM1 and TM13 in NdhD, which was ascribed to the substitution from motif 11 to 5. Moreover, it lacked the long amphipathic helix in NdhF that was attributed to the absence of motifs 7, 14, and 17.

3.4 The splicing patterns of ndh in different tissues and under different light stress

The alternative splicing (AS) types of the ndh genes that contained introns, including the chloroplast-encoded genes (ndhA, ndhB) and most of the nuclear-encoded genes (with the exception of ndhS and pnsl2), were analyzed using RT-PCR (Fig. 5 and Supplementary Fig.S2). Most types of AS events among these five tissues and leaves after exposure to light stress was intron retention. Ten ndh genes (ndhB, ndhL, ndhM, ndhU, ndhV, pnsl1, pnsl4, pnsb2, pnsb3-1, and pnsb3-2) showed tissue-specific AS patterns (Fig. 5a), indicating that the AS of ndh genes possibly participated in the regulation of functional specification in different tissues. Remarkably, the AS transcript of ndhB (ndhB-AS-2) that consisted of exons 1 and 2 and intron 1 exhibited the opposite pattern of expression with the fundamental transcript (ndhB-1). The ndhB-AS-2 was the primary isoform in stems and rhizomes, while ndhB-1 was preferentially expressed in leaves. Similarly, pnsb2 produced two transcripts: the functional transcript (pnsb2-1) that was the primary isoform in leaves and flowers and the unspliced transcript (pnsb2-2) that was the main isoform in rhizomes, stems, and roots. Furthermore, most AS transcripts of the genes described above did not respond to the light stress with the exception of ndhM whose AS transcripts (ndhM-2) turned into the dominant isoform after high light treatment (Fig. 5b). Subsequently, the gene structures of different transcript isoforms and their influence on protein translation were analyzed. As shown in Fig. 6, the AS transcription either led to functional impairment (ndhB-AS-2, ndhL-2, ndhM-2, ndhU-2, ndhU-3, ndhU-4, ndhV-2, pnsb2-2, pnsb3-2-2, pnsb3-2-3, pnsb5-2, pnsl1-2, pnsl1-3, and pnsl1-4), or untranslated mRNAs (ndhL-3 and psnb3-1-2).

3.5 RNA editing analysis of the chloroplastencoded ndh genes

A total of 35 editing sites were detected in 11 plastid-encoded ndh genes (Table 2), with partial editing of ndhD and ndhF and a lack of editing for ndhC. Compared with other species, Z. marina possessed more editing sites than those in Spirodela polyrhiza (the close relative of the seagrass), A. thaliana (the model dicot) and Z. mays (the model monocot) but fewer than those in A. trichopoda (the basal angiosperm). Z. marina shared more editing sites with A. trichopoda (12) than with S. polyrhiza (10), A. thaliana (8) and Z. mays (6). Furthermore, more types of editing, including C to U (48.6%), U to C (20%), A to G (14.3%), G to A (8.6%), C to A (5.7%), and G to U (2.8%), existed in Z. marina, while only the C to U editing existed in the other plants. Many of the primitive editing sites and types suggested that there was an ancestral editing pattern in Z. marina. Additionally, most of the editing events tended to occur in the third bases of codons that led to silent editing. The second bases of codons were also edited that resulted in the alteration of identity of amino acids. The majority of non-synonymous editing events in Z. marina either restored evolutionarily conserved amino acid sequences, such as the methionine to isoleucine conversion of the ndhB 72^th codon, the isoleucine to methionine conversion of the ndhD 435^th codon and the asparagine to aspartic acid conversion of the ndhH 50^th codon, or generated lineage-specific residues in Z. marina compared with S. polyrhiza, A. thaliana, and Z. mays, i.e., the ndhD 220^th, ndhI 45^th, and ndhK 51^th codons (Fig. 7a & b; Supplementary Figs.S3–S4).

The RNA editing events in different illumination conditions were also detected in the ndh genes of Z. marina. Interestingly, partial editing of the 361^th codon in ndhF that resulted in the incomplete serine to tyrosine conversion was observed in the dark adaptation period but subsequently disappeared in the light exposure and recovery periods (Fig. 7c; Supplementary Fig.S4), which reflected the response of RNA editing to the light exposure.

3.6 Expression patterns of ndh in different tissues

The expression of all 31 ndh genes derived from the five types of tissues of Z. marina were also investigated, with the level of expression of the ndhB representing the sum of duplicated ndhB owing to their consensus coding sequences. As revealed in Fig. 8, plastid-encoded genes (ndhA-ndhK) were preferentially expressed in leaves but barely expressed in roots and rhizomes. Moreover, nuclear-encoded genes exhibited various patterns of expression across the five tissues. Five genes (ndhL, ndhU, ndhS, pnsl2, and pnsl4) in the roots, six genes (pnsb1, pnsb2, pnsb3-1, pnsb4, pnsb5, and pnsl1) in the flowers and seven genes (ndhM, ndhN, ndhO, ndhT, ndhV, pnsb3-2, and pnsl3) in the leaves exhibited relatively high transcript abundances, while pnsl5 was constitutively expressed. It was notable that pnsb3-1 was highly expressed in flowers. In contrast, its duplicated sister (pnsb3-2) was particularly abundant in leaves, implying a functional divergence in the pnsb3 gene pairs. Unlike the other ndh genes, pnsl4 and ndhL showed different patterns of expression with their counterparts in A. thaliana (Supplementary Fig.S5), suggesting their functional divergence in different species.

3.7 Expression patterns of the ndh under the high light radiation and dark recovery

To analyze the dynamic expression of ndh genes, the 11 undetectable (ndhA-K) and eight differentially expressed ndh genes (ndhL, ndhM, ndhS, ndhT, ndhV, pnsb3-1, pnsb3-2, and pnsl5) in the available transcriptome were selected (Supplementary Fig.S6). As shown in Fig. 9, four ndh genes in conserved Subcomplex M (ndhA), Subcomplex A (ndhJ and ndhK) and Subcomplex EDB (ndhV) were rapidly induced and maintained a relatively high level of expression throughout their exposure to light. Moreover, three other genes in Subcomplex M (ndhC, ndhD, and ndhE), as well as the two genes in Subcomplex B (pnsb3-2) and Subcomplex L (pnsl5) that exist exclusively in terrestrial plants, were up-regulated during the middle and late stages of light exposure. Furthermore, these genes were also up-regulated during the subsequent dark recovery period. In contrast, ndhB, ndhH, ndhI, ndhL, ndhM, and pnsb3-1 were repressed during light exposure and then significantly induced during the recovery period. The dynamic patterns of expression of the ndh genes indicated that these conserved subcomplexes during evolution could be the most sensitive parts involved in the activation of the NDH complex and there were two distinct responsive and regulatory mechanisms of the NDH complex at transcriptional level during the periods of light exposure and the subsequent dark recovery.

3.8 Cis-elements analysis for the ndh promoters

To explore the possible regulatory mechanism of the Z. marina ndh genes, related cis-elements were scanned. There were three types of cis-elements, as shown in Fig. 10. Different amounts of light-responsive elements, including Box4, G-box, the GATA, GT1 and TCT motifs, were present in each ndh gene, with the G-box comprising the most abundant element. Among the 31 ndh genes, pnsl5 and ndhM had more light responsive elements. The second type of cis-elements was related to hormones, including those responsive to salicylic acid-, abscisic acid- and auxin, which were detected among 11 genes (ndhA, ndhG, ndhI, ndhM, ndhO, ndhS, pnsb1, pnsb2, pnsb3-2, pnsb5, and pnsl4), indicating that the hormones also played roles in regulating ndh genes expression. In addition, defense- and stress-responsive and WUN-motif and TC-rich repeats associated with wound responses were identified in nine genes (ndhF, ndhH, ndhI, ndhM, ndhU, pnsb3-2, pnsl1, pnsl2, and pnsl5), which contributed to plant tolerance to various stresses. These results indicated that the ndh genes can be regulated by numerous factors related to plant resistance with light serving as the primary factor.

The upstream length to the translational start site refers to the scale at the bottom.

4 DISCUSSION

The loss of ndh genes was universal in marine macrophytes (Peredo et al., 2013; Lee et al., 2018), which was consistent with that the complexity of NDH gradually increased during the transition from marine plants to land plants. Considering that most of the ndh mutants resulted in impaired NDH activity (Kofer et al., 1998; Horvath et al., 2000; Rumeau et al., 2005; Ishikawa et al., 2008; Fan et al., 2015), it was believed that each Ndh subunit was indispensable, functioning either in the stability or assembly of the NDH complex. In many cases, the contribution of the NDH dependent pathway is dispensable (Munekage et al., 2004). However, it could alleviate stromal over-reduction under stress conditions (Munekage et al., 2004; Okegawa et al., 2008). The NDH complex in most marine photosynthetic organisms was absent, which could be associated with their relatively stable habitats compared with terrestrial environments. Z. marina, as a marine seagrass, possesses a complete NDH that could possibly have been retained during the migration from land to sea. It could be associated with its special physiology and provides the foundation for highly efficient NDH dependent PSI-CEF and chlororespiration.

The comparison of the phylogeny between individual Ndh subunits and the entire NDH complex revealed that some genes were vertically or horizontally transferred among Viridiplantae together with the entire NDH complex, although some ndh genes evolved independently. The earlier diversification of NdhM, NdhO, and NdhS suggests a possible more primitive evolutionary status of these subunits in Z. marina.

Among the 31 ndh genes, two gene pairs (ndhB-1 and ndhB-2, chloroplast genome encoded; pnsb3-1 and pnsb3-2, nuclear genome encoded) probably originated from gene duplication events. It is commonly accepted that the majority of gene duplications in chloroplast genomes are caused by the expansion of large IRs (Xiong et al., 2009). In this study, there is only one ndhB in the outgroup species, while the duplicated ndhB paralogs in Cycas (the most ancient extant plants) and the angiosperms were all located on the chloroplast IRs (Wu et al., 2007; Xing and Guo, 2018). This implied that the duplication of ndhB might emerge in the common ancestors of seed plants through the expansion of IRs. However, the phenomenon that some gymnosperms (Ginkgoaceae and Pinaceae) also possessed a single copy of ndhB (Zhang et al., 2014) is attributed to the loss or contraction of IRs. Nuclear-encoded pnsb3 paralogs occurred as "single gene transposition-duplication". This type of duplication could generate two gene copies that are neither neighboring nor colinear. The non-duplicated pnsb3 in algae, bryophytes, lycophytes, and gymnosperms, together with the independent orthologous groups of the duplicated pnsb3 (pnsb3-1 and pnsb3-2) in the monocots and dicots, indicated that the duplication events happened after the divergence of gymnosperms and angiosperms but before the divergence of monocots and dicots. In Z. marina, both paralogs possessed the core angiosperm motifs 1, 2, 3 and 4, while the additional motifs 5, 10, and 15 were specific to pnsb3-1. The divergent motif composition combined with the different tissue-specific pattern of expression suggested that the neo-functionalization and/or sub-functionalization of pnsb3 paralogs could enable plants to adapt to coastal environments (He and Zhang, 2005).

Motifs in the functional domains were conserved in most of the Z. marina Ndh subunits with the exception of NdhD and NdhF. Previous studies (Efremov et al., 2010) showed that there is one near-continuous long amphipathic α-helix between the 15^th and 16^th TM helices of NdhF connected to other antiporter-like subunits. This could transmit conformational changes in the hydrophilic domain of NDH for proton translocation. The absence of the TM helices in NdhD and NdhF imply a possible modified mechanism for proton transmission in Z. marina. Moreover, the non-photochemical quenching value, which depends on the pH gradient generated in the thylakoid, is at a low degree in many seagrasses (Schubert et al., 2015). Thus, there could possibly be a special proton transmission mechanism that is related to the low non-photochemical quenching capacity.

Many AS transcripts were predominantly observed in leaves that primarily contained chloroplasts. Clearly, ndh are involved with photosynthetic-related regulation. Many splice variants generate small interfering peptides that participate in the formation of multi-protein complex but lack functional domains to compete with the functional complex (Staudt and Wenkel, 2011; Reddy et al., 2013). Intriguingly, the AS transcripts of ndhB, ndhM, and pnsb2, yielded truncated proteins without functional domains and exhibited opposite patterns of expression with their fundamental transcriptions indicating one negative regulation in these ndh genes.

RNA editing, another post-transcriptional process, caused numerous base transitions (Chateigner-Boutin and Small, 2010). In addition, the canonical C to U editing and the abundant uncanonical U to C, A to G, and G to A conversions that had been reported to exist exclusively in ancestral land plants (Chateigner-Boutin and Small, 2010; Uthaipaisanwong et al., 2012) were unexpectedly present in Z. marina. Moreover, the minor G to T and C to A editing events that had never been reported in other plant species were detected. It is possible that Z. marina retained numerous ancestral editing patterns that had a monophyletic origin but were followed by lineage-specific losses and gains. As suggested by Jobson and Qiu (2008), the RNA editing events could change the protein structure or interaction through the non-synonymous replacement of conserved amino acids (Jobson and Qiu, 2008; Chen et al., 2011). Therefore, the editing events adjacent to the charged lysine residue in the TM7 helix of NdhD, which is crucial for energy transduction (Efremov et al., 2010), as well as in the second helix of NdhK and the N-terminal amphiphilic helix of NdhI, which had an impact on quinone-binding (Pan et al., 2020), could influence the functional NDH in Z. marina.

Similar to A. thaliana, most of the ndh genes in Z. marina were primarily expressed in leaves, which implied their conserved function in photosynthesis of the NDH complex during evolution. However, the prolyl cis/trans isomerases, pnsl4 and pnsl5, could be involved in various physiological processes in addition to the assembly of NDH complex (Sirpiö et al., 2009), such as hormone-mediated plant development and brassinosteroid-mediated flowering (Romano et al., 2005; Zhang et al., 2013), thereby leading to their unbiased expression in leaves. Additionally, the ndhL in Z. marina was predominantly expressed in roots and flowers, which differed from that in A. thaliana. Enriched bicarbonate in seawater can be converted to CO₂ in the surface of leaves to support the concentration of limited inorganic carbon (Ci) in fully marine conditions (Larkum et al., 2017). Therefore, the high level of expression of ndhL, which functioned in the concentration and transport of Ci (Zhang et al., 2005; Shimizu et al., 2008), in roots and flowers might contribute to the transport of Ci to the leaves in Z. marina.

Most of the ndh-encoded components in Subcomplex M, Subcomplex A, and Subcomplex EDB were significantly induced or repressed during different light periods. Thus, the three subcomplexes could be the most sensitive parts involved in the activation of NDH complex. More light responsive elements in the promoter regions of pnsl5 and ndhM could be related to the high and continuous level of expression of these genes that responded to light stress. In addition, the up-regulated ndh genes differed during the light exposure and dark recovery period. This suggested that there were two diverse responsive and regulatory mechanisms at the transcriptional level in NDH-dependent PSI-CEF and chlororespiration. The complete structure, upregulated gene expression level and multiple post-transcriptional regulations could provide a molecular basis for the highly efficient NDH-dependent PSI-CEF, which contributed to the generation of ΔpH and repair of the photo-inactivated oxygen-evolving complex (Tan et al., 2020a, b), to maintain effective photosynthetic performance.

5 CONCLUSION

Thirty-one ndh genes were identified, of which ndhB and pnsb3 occurred as duplication events during evolution. The long amphipathic helix in NdhF was lost, which reflects a possible alternative mode in the generation of trans-thylakoid proton gradient. The AS events and RNA editing events showed the transcriptional regulatory mechanism of the NDH complex during the light stress. The RNA editing in Z. marina exhibited the ancestral pattern with many of the primitive editing sites and types. The dynamic profiles of expression in response to light stress suggested that there were two diverse responsive mechanisms of the NDH complex in PSI-CEF and chlororespiration.

6 DATA AVAILABILITY STATEMENT

Individual sequences were submitted to GenBank (https://www.ncbi.nlm.nih.gov/genbank/) and can be retrieved with accession numbers MW051562 and MW051563. The voucher specimen (specimen number: HY202005) is available in the Herbarium of Ocean School of Yantai University, Shandong Province, China.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S2 and Figs.S1–S6) is available in the online version of this article at https://doi.org/10.1007/s00343-021-0027-z.