A genome-wide analysis of the chloroplast N ADH dehydrogenase- like genes in Zostera marina*

Mingyu MA, Mingyu ZHONG, Quansheng ZHANG, Wei ZHAO, Mengxin WANG,Chengying LUO, Bin XU

Ocean School, Yantai University, Yantai 264005, China

Abstract The chloroplast NADH dehydrogenase-like (NDH) complex, homologous to respiratory complex I, participates in photosystem I cyclic electron f low (PSI-CEF) and chlororespiration in photosynthesis.Phylogenetic analyses indicated that Zostera marina, a widely distributed seagrass, has a complete NDH complex, which is rarely observed in marine macrophytes. We identif ied all 31 ndh genes necessary for the functional NDH complex, of which ndhB and pnsb3 occurred as duplication events. Secondary structural analyses of antiporter-like subunits showed that the long amphipathic helix of NdhF was lost in Z. marina,which could exhibit an alternative mode in the generation of trans-thylakoid proton gradient. The splicing pattern of ndh exhibited tissue-specif ic patterns and responded to light stress. RNA editing in Z. marina presented the ancestral pattern with many of the primitive editing sites and types. The partial editing in ndhF ref lected the link between light stress and RNA editing. Moreover, the predominant expression in leaves of most ndh genes suggested that their major function is in photosynthesis. The quantitative real time-PCR results show that the expression of ndh was signif icantly upregulated in response to light stress.Nevertheless, there were two diverse responsive mechanisms of the NDH complex in PSI-CEF and chlororespiration. Overall, the presence of a complete structure, upregulated gene expression level, and multiple post-transcriptional regulations could provide a molecular basis for the powerful NDH complex and enable Z. marina to maintain eff ective photosynthetic performance.

Keyword: Zostera marina; NADH dehydrogenase-like complex; genome-wide analysis; prof iles of expression

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 f lavin 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 f ive subcomplexes, including subcomplex M (NdhANdhG), subcomplex A (NdhH-NdhO), subcomplex B(PnsB1-PnsB5), subcomplex L (PnsL1-PnsL5), and subcomplex EDB (Electron donor binding) (NdhSNdhV) (Ifuku et al., 2011).

Since the defectivendhBwas verif ied to inf luence photosystem I cyclic electron transport (PSI-CEF) inSynechocystissp. PCC 6803 (Ogawa, 1991),chloroplast NDH was proven to be involved in PSICEF (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 PSICEF 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 f lux via Fd and eff ectively 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 inEscherichiacoli, cyanobacteria, and land plants (Efremov et al., 2010). Although plastidencodedndh(ndhA-ndhK) are found in some ancestral algae, e.g.,NephroselmisolivaceaorMesostigmaviride(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).Marchantiapolymorpha, 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 specif ic 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.

Zosteramarina(Alismatales: Zosteraceae) is one of the most productive types of seagrass and widely distributed in the northern Pacif ic 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 31ndhgenes of the entire NDH complex were identif ied inZ.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 thatZ.marinapossesses a highly effi cient 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 specif ically-evolved and highly effi cientZ.marinaNDH has yet to be investigated to our knowledge. The whole genome sequences that have been reported enabled us to conduct a systematic research ofndhgenes inZ.marina(Olsen et al., 2016;Xing and Guo, 2018). We comprehensively analyzed the phylogeny, gene structures, patterns of expression,and regulation ofndhgenes inZ.marina. This study provided clues to understand how thendhcoordinate with each other to drive the highly effi cient NDHdependent PSI-CEF in response to light exposure and chlororespiration.

2 MATERIAL AND METHOD

2.1 Plant material and treatment

Zosteramarinawith 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 identif ied by experienced taxonomists on the basis of their morphology. The identif ications were also conf irmed by sequence similarity with the whole genome sequencing ofZ.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 °C with a photoperiod of 10 h:14 h (light:dark) in minimum saturation light intensity(100 μmol photons/(m2·s)).

The leaves obtained from dark-adapted overnight plants were exposed to 300 μmol photons/(m2·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, f lowers, stems, and rhizomes were collected at the f lowering stage for tissue-specif ic 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 theZ.marinadatabase from the NCBI (https://www.ncbi.nlm.nih.gov/) and Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) databases using theArabidopsisthalianaNdh 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 eachndhgene were conf irmed using Pfam and SMART (http://smart.embl-heidelberg.de/). The exon-intron compositions of thendhgenes were constructed using DNAMAN software by comparing the cDNA with genomic sequences. The duplication modes of the gene pairs were identif ied 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, usingE.colias 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 identif ied 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, includingThermosynechococcuselongatusBP-1,Zeamays,Sorghumbicolor,A.thaliana,Populustrichocarpa, andZ.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 quantif ied using NanoQuant (TECAN Group Ltd., Männedorf, Switzerland). After wiping off the residual DNA, the cDNA was synthesized with a HiScript®II 1st Strand cDNA Synthesis Kit (Vazyme)using 1 μg of total RNA.

2.6 PCR amplif ication and sequencing

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

2.7 Analysis of splicing patterns

The RT-PCR products ofndhgenes from f ive tissues and leaves exposed to diff erent light periods were visualized using a Gel Doc XR+ system (Bio-Rad, Hercules, CA, USA). The specif ic electrophoretic bands corresponded to diff erently spliced products.The possible protein isoforms were predicted by DNAMAN.

2.8 RNA editing analysis of chloroplast-encoded ndh genes

To conf irm the RNA editing sites in theZ.marinandhgenes, the sequenced cDNA and DNA from diff erent light periods were aligned using Vector NTI.The editing sites forSpirodelapolyrhiza(Wang et al.,2015),A.thaliana(Tillich et al., 2005),Z.mays(Maier et al., 1995), andAmborellatrichopoda(Hein et al., 2016) were obtained from respective publications.

2.9 Analysis of gene expression

Quantitative real time-PCR assays of thendhgenes obtained from f ive tissues and diff erent light periods were conducted on a Bio-Rad CFX96 Real Time PCR System with AceQ Universal SYBR qPCR Master Mix (Vazyme). The housekeeping genegapdhfromZ.marinawas used as an internal control. The qRTPCR was programmed as follows: 95 °C for 10 s,followed by 40 cycles of 56 °C for 10 s and 72 °C for 30 s. The qRT-PCR data were calculated using the 2-ΔΔCtmethod. Sequences of the primers used in qRTPCR are described in Supplementary Table S1. The tissue-specif ic expression ofA.thalianawas downloaded as the publicly available RNA-Seq data in the CoNekT-Plants database (www.conekt.plant.tools). The RNA Seq ofZ.marinaunder light stress was analyzed in our previous study (Tan et al., 2020b).The diff erentially expressed genes were selected according toP＜0.05 and fold change ＞1.3 (for upregulation) or fold change ＜0.75 (for downregulation). 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 possiblecis-acting elements in the promoter regions of thendhgenes, 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 withP＜0.05 set to be statistically signif icant.

3 RESULT

3.1 Identif ication of the Ndh subunits

A total of 31 Ndh subunits homologous to theirA.thalianacounterparts were obtained fromZ.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 ofZ.marina. GeneZosma266g00030.1encodes the NdhV subunit andZosma266g00040.1encodes a hypothetical protein that are part of the same transcript (designatedZosma266g0003040.1). Following the conf irmation by RT-PCR and sequencing, in silico translation revealed thatZosma266g0003040.1translated a 191 amino acid protein. Furthermore, the misannotated coding sequence ofndhBwas observed during the analysis of sequenced cDNA fragments in the section of RNA editing. The duplicatedndhBhad a uniform coding sequence with a 509 amino acid protein. To explore the evolutionary clues of duplicatedndhgene pairs, f ive widely accepted duplication modes(whole genome, tandem, proximal, transposed, and dispersed duplication) were systematically scanned.The gene pairsndhB-1/ndhB-2andpnsb3-1/pnsb3-2originated from the expansion of large inverted repeat regions (IRs) and transposed duplication,respectively.

Fig.1 Gene model of ndh genes in Zostera marina

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 diversif ication between monocots and dicots. Unlike other marine species,Z.marinahad 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 ofZ.marinawere at the basal branches of angiosperms (Supplementary Fig.S1b), indicating their earlier evolution status.

Table 1 Features of ndh genes identif ied in Zostera marina

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 nuclearencoded Ndh subunits. Some are clearly unique to theZ.marinaNdh 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.marinalost 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.

Fig.2 Phylogenetic relationship of NADH dehydrogenase-like (NDH) complex among Viridiplantae

Fig.3 Phylogenetic trees and motif compositions of NdhD (a), F (b), G (c), H (d), and I (e) from 19 diff erent species

Fig.3 Continued

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

The alternative splicing (AS) types of thendhgenes that contained introns, including the chloroplastencoded genes (ndhA,ndhB) and most of the nuclearencoded genes (with the exception ofndhSandpnsl2),were analyzed using RT-PCR (Fig.5 and Supplementary Fig.S2). Most types of AS events among these f ive tissues and leaves after exposure to light stress was intron retention. Tenndhgenes (ndhB,ndhL,ndhM,ndhU,ndhV,pnsl1,pnsl4,pnsb2,pnsb3-1, andpnsb3-2) showed tissue-specif ic AS patterns(Fig.5a), indicating that the AS ofndhgenes possibly participated in the regulation of functional specif ication in diff erent tissues. Remarkably, the AS transcript ofndhB(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). ThendhB-AS-2was the primary isoform in stems and rhizomes, whilendhB-1was preferentially expressed in leaves. Similarly,pnsb2produced two transcripts: the functional transcript (pnsb2-1) that was the primary isoform in leaves and f lowers 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 ofndhMwhose AS transcripts (ndhM-2) turned into the dominant isoform after high light treatment (Fig.5b).Subsequently, the gene structures of diff erent transcript isoforms and their inf luence 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, andpnsl1-4), or untranslated mRNAs (ndhL-3andpsnb3-1-2).

Fig.4 ClustalW amino acid alignment of the amino acid sequence of NdhD (a) and NdhF (b) from Z. marina with an additional f ive species

Fig.4 Continued

3.5 RNA editing analysis of the chloroplastencoded ndh genes

A total of 35 editing sites were detected in 11 plastid-encodedndhgenes (Table 2), with partial editing ofndhDandndhFand a lack of editing forndhC. Compared with other species,Z.marinapossessed more editing sites than those inSpirodela polyrhiza(the close relative of the seagrass),A.thaliana(the model dicot) andZ.mays(the modelmonocot) but fewer than those inA.trichopoda(the basal angiosperm).Z.marinashared more editing sites withA.trichopoda(12) than withS.polyrhiza(10),A.thaliana(8) andZ.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 inZ.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 inZ.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 inZ.marinaeither restored evolutionarily conserved amino acid sequences, such as the methionine to isoleucine conversion of thendhB72thcodon, the isoleucine to methionine conversion of thendhD435thcodon and the asparagine to aspartic acid conversion of thendhH50thcodon, or generated lineage-specif ic residues inZ.marinacompared withS.polyrhiza,A.thaliana, andZ.mays, i.e., thendhD220th,ndhI45th, andndhK51thcodons (Fig.7a & b;Supplementary Figs.S3-S4).

Table 2 Comparison of identif ied RNA editing sites in the ndh genes of Z. marina with those of other species

Fig.6 Gene model of the diff erent isoforms of ndh genes occurring AS events

The RNA editing events in diff erent illumination conditions were also detected in thendhgenes ofZ.marina. Interestingly, partial editing of the 361thcodon inndhFthat 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 ref lected the response of RNA editing to the light exposure.

3.6 Expression patterns of ndh in diff erent tissues

Fig.7 Sequencing chromatograms of RNA editing sites that restored evolutionarily conserved amino acid sequences (a),generated lineage-specif ic residues in Z. marina (b), and inf luenced by light stress (c)

The expression of all 31ndhgenes derived from the f ive types of tissues ofZ.marinawere also investigated,with the level of expression of thendhBrepresenting the sum of duplicatedndhBowing to their consensus coding sequences. As revealed in Fig.8, plastidencoded 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 f ive tissues.Five genes (ndhL,ndhU,ndhS,pnsl2, andpnsl4) in the roots, six genes (pnsb1,pnsb2,pnsb3-1,pnsb4,pnsb5, andpnsl1) in the f lowers and seven genes(ndhM,ndhN,ndhO,ndhT,ndhV,pnsb3-2, andpnsl3)in the leaves exhibited relatively high transcript abundances, whilepnsl5was constitutively expressed.It was notable thatpnsb3-1was highly expressed in f lowers. In contrast, its duplicated sister (pnsb3-2) was particularly abundant in leaves, implying a functional divergence in thepnsb3gene pairs. Unlike the otherndhgenes,pnsl4andndhLshowed diff erent patterns of expression with their counterparts inA.thaliana(Supplementary Fig.S5), suggesting their functional divergence in diff erent species.

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

To analyze the dynamic expression ofndhgenes,the 11 undetectable (ndhA-K) and eight diff erentially expressedndhgenes (ndhL,ndhM,ndhS,ndhT,ndhV,pnsb3-1,pnsb3-2, andpnsl5) in the available transcriptome were selected (Supplementary Fig.S6).As shown in Fig.9, fourndhgenes in conserved Subcomplex M (ndhA), Subcomplex A (ndhJandndhK) and Subcomplex EDB (ndhV) were rapidlyinduced and maintained a relatively high level of expression throughout their exposure to light.Moreover, three other genes in Subcomplex M (ndhC,ndhD, andndhE), as well as the two genes in Subcomplex B (pnsb3-2) and Subcomplex L (pnsl5)that exist exclusively in terrestrial plants, were upregulated during the middle and late stages of light exposure. Furthermore, these genes were also upregulated during the subsequent dark recovery period.In contrast,ndhB,ndhH,ndhI,ndhL,ndhM, andpnsb3-1were repressed during light exposure and then signif icantly induced during the recovery period.The dynamic patterns of expression of thendhgenes 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.

Fig.7 Continued

3.8 C is-elements analysis for the ndh promoters

Fig.8 Patterns of expression of ndh genes in f ive tissues analyzed by quantitative real time-PCR (qRT-PCR)

To explore the possible regulatory mechanism of theZ.marinandhgenes, relatedcis-elements were scanned. There were three types ofcis-elements, as shown in Fig.10. Diff erent amounts of lightresponsive elements, including Box4, G-box, the GATA, GT1 and TCT motifs, were present in eachndhgene, with the G-box comprising the most abundant element. Among the 31ndhgenes,pnsl5andndhMhad more light responsive elements. The second type ofcis-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, andpnsl4), indicating that the hormones also played roles in regulatingndhgenes expression. In addition, defense- and stress-responsive and WUN-motif and TC-rich repeats associated with wound responses were identif ied in nine genes (ndhF,ndhH,ndhI,ndhM,ndhU,pnsb3-2,pnsl1,pnsl2, andpnsl5), which contributed to plant tolerance to various stresses. These results indicated that thendhgenes can be regulated by numerous factors related to plant resistance with light serving as the primary factor.

Fig.9 Patterns of expression of ndh genes in light exposure and subsequent recovery analyzed using qRT-PCR

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

4 DISCUSSION

Fig.10 Predicted cis-elements in ndh genes promoters. Promoter sequences (-1 500 bp) of 31 ndh genes are detected by PlantCARE

The loss ofndhgenes 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 thendhmutants 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 overreduction 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 effi cient NDH dependent PSICEF 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 somendhgenes evolved independently. The earlier diversif ication of NdhM, NdhO, and NdhS suggests a possible more primitive evolutionary status of these subunits inZ.marina.

Among the 31ndhgenes, two gene pairs (ndhB-1andndhB-2, chloroplast genome encoded;pnsb3-1andpnsb3-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 onendhBin the outgroup species,while the duplicatedndhBparalogs inCycas(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 ofndhBmight 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 ofndhB(Zhang et al., 2014) is attributed to the loss or contraction of IRs. Nuclear-encodedpnsb3paralogs occurred as “single gene transpositionduplication”. This type of duplication could generate two gene copies that are neither neighboring nor colinear. The non-duplicatedpnsb3in algae,bryophytes, lycophytes, and gymnosperms, together with the independent orthologous groups of the duplicatedpnsb3(pnsb3-1andpnsb3-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. InZ.marina, both paralogs possessed the core angiosperm motifs 1, 2, 3 and 4,while the additional motifs 5, 10, and 15 were specif ic topnsb3-1. The divergent motif composition combined with the diff erent tissue-specif ic pattern of expression suggested that the neo-functionalization and/or sub-functionalization ofpnsb3paralogs could enable plants to adapt to coastal environments (He and Zhang, 2005).

Motifs in the functional domains were conserved in most of theZ.marinaNdh subunits with the exception of NdhD and NdhF. Previous studies(Efremov et al., 2010) showed that there is one nearcontinuous long amphipathic α-helix between the 15thand 16thTM 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 modif ied mechanism for proton transmission inZ.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,ndhare 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 ofndhB,ndhM, andpnsb2, yielded truncated proteins without functional domains and exhibited opposite patterns of expression with their fundamental transcriptions indicating one negative regulation in thesendhgenes.

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 inZ.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 thatZ.marinaretained numerous ancestral editing patterns that had a monophyletic origin but were followed by lineagespecif ic losses and gains. As suggested by Jobson and Qiu (2008), the RNA editing events could change the protein structure or interaction through the nonsynonymous 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 inf luence the functional NDH inZ.marina.

Similar toA.thaliana, most of thendhgenes inZ.marinawere primarily expressed in leaves, which implied their conserved function in photosynthesis of the NDH complex during evolution. However, the prolylcis/transisomerases,pnsl4andpnsl5, 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 f lowering(Romano et al., 2005; Zhang et al., 2013), thereby leading to their unbiased expression in leaves.Additionally, thendhLinZ.marinawas predominantly expressed in roots and f lowers, which diff ered from that inA.thaliana. Enriched bicarbonate in seawater can be converted to CO2in 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 ofndhL, which functioned in the concentration and transport of Ci(Zhang et al., 2005; Shimizu et al., 2008), in roots and f lowers might contribute to the transport of Ci to the leaves inZ.marina.

Most of thendh-encoded components in Subcomplex M, Subcomplex A, and Subcomplex EDB were signif icantly induced or repressed during diff erent 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 ofpnsl5andndhMcould be related to the high and continuous level of expression of these genes that responded to light stress. In addition, the up-regulatedndhgenes diff ered 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 effi cient 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 eff ective photosynthetic performance.

5 CONCLUSION

Thirty-onendhgenes were identif ied, of whichndhBandpnsb3occurred as duplication events during evolution. The long amphipathic helix in NdhF was lost, which ref lects 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 inZ.marinaexhibited the ancestral pattern with many of the primitive editing sites and types. The dynamic prof iles 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.