Developmental Genetic Mechanisms of C4 Syndrome Based on Transcriptome Analysis of C3 Cotyledons and C4 Assimilating Shoots in Haloxylon ammodendron

It is believed that transferring the C4 engine into C3 crops will greatly increase the yields of major C3 crops. Many efforts have been made since the 1960s, but relatively little success has been achieved because C4plant traits, referred to collectively as C4 syndrome, are very complex, and little is known about the genetic mechanisms involved. Unfortunately, there exists no ideal genetic model system to study C4 syndrome. It was previously reported that the Haloxylon species have different photosynthetic pathways in different photosynthetic organs, cotyledons and assimilating shoots. Here, we took advantage of the developmental switch from the C3 to the C4 pathway to study the genetic mechanisms behind this natural transition. We compared the transcriptomes of cotyledons and assimilating shoots using mRNA-Seq to gain insight into the molecular and cellular events associated with C4 syndrome. A total of 2959 differentially expressed genes [FDR≤0.001 and abs (|log2(Fold change)|≥1)] were identified, revealing that the transcriptomes of cotyledons and assimilating shoots are considerably different. We further identified a set of putative regulators of C4 syndrome. This study expands our understanding of the development of C4 syndrome and provides a new model system for future studies on the C3-to- C4 switch mechanism.


Introduction
Photosynthetic CO 2 fixation is a fundamental life process involving the conversion of solar energy into chemical energy that can be later released to fuel the activity of an organism. The ancestral photosynthetic CO 2 -fixation process is C 3 photosynthesis. The first organic product of CO 2 fixation is a three-carbon compound. C 3 photosynthesis and its key enzyme Rubisco (ribulose-1,5-bisphosphate carboxylase/oxygenase) evolved early in the history of life, when there was no oxygen in the atmosphere and atmospheric CO 2 levels were significantly high. As atmospheric CO 2 dropped and O 2 accumulated, Rubisco inefficiency began to limit C 3 photosynthesis because the active site of Rubisco does not completely discriminate between CO 2 and O 2 , leading to the catalysis of two competitive reactions: photosynthetic CO 2 assimilation and photorespiratory CO 2 loss, especially under hot, dry and/or saline conditions that enhance photorespiration. Approximately 30 million years ago, an abrupt drop in atmospheric [CO 2 ] further reduced the efficiency of C 3 photosynthesis and triggered the evolution of C 4 photosynthesis, in which CO 2 is initially fixed into a four-carbon compound and concentrated around Rubisco [1][2][3][4][5][6]. C 4 photosynthesis has evolved >60 times and occurs in approximately 7500 species of flowering plants [7][8][9]. The transition from C 3 to C 4 plants was the green revolution of nature. Although C 4 plants comprise only 3% of land plant species, they account for some 25% of global terrestrial carbon fixation [3,4,7,10,11]. C 4 photosynthesis is a complex trait that combines biochemical, physiological and anatomical characteristics, the so-called C 4 syndrome [12,13]. Other than four single-celled C 4 lineages, the vast majority of C 4 plants possess Kranz anatomy and C 4 syndrome [8,9], indicative of convergent evolution [7]. Traditional biochemical and modern molecular biological studies showed that all of the proteins required for the core C 4 cycle are present in C 3 plants [14]. Based on this evidence, it is hypothesized that no significant genetic changes are required for the transition from C 3 to C 4 photosynthesis [15,16]. The first attempt at introducing C 4 photosynthesis into C 3 plants was made in Atriplex through conventional interspecific hybridization of photosynthetic types to reduce photorespiration and increase photosynthetic capacity [17][18][19][20][21][22]. The results indicated that F1 hybrids of Atriplex rosea (C 4 , NAD-ME type) × Atriplex triangularis (C 3 ) were more similar to their C 3 parents in physiology and failed to form well-developed Kranz anatomy. Similar hybridization studies were conducted in the genera Flaveria, Panicum, Moricandia, and Brassica, but none proved fruitful at converting C 3 species into functional C 4 types [23]. Due to infertility, few hybrids have been developed beyond the F1 generation. In the few advanced generations studied in Atriplex and Flaveria hybrids, correlations among photosynthetic traits were low, indicating that C 4 photosynthesis is a combination of independent biochemical, physiological and anatomical characteristics [23]. Recently, a progress report was released on oat-maize addition lines showing that the addition of individual maize chromosomes to the C 3 species oat caused increases in vein density but did not confer functional C 4 photosynthesis [24]. Therefore, introducing a fully developed C 4 photosynthesis pathway into C 3 plants through interspecific hybridization or even genetic engineering is far more practical than previously thought [25].
In nature, there exist examples of switches from a C 3 pathway to a two-celled C 4 pathway triggered by external or internal signals. The former example includes the freshwater amphibious leafless sedge Eleocharis vivipara, which can switch from a C 3 pathway to a C 4 pathway with Kranz anatomy after induction by environmental changes and exogenous application of abscisic acid (ABA) [26][27][28]. The latter examples were reported in Haloxylon and Salsola species that have different photosynthetic pathways in different photosynthetic organs [29,30]. Haloxylon aphyllum and H. persicum of Chenopodiaceae have C 3 photosynthesis in cotyledons and C 4 photosynthesis in assimilating shoots (the main photosynthetic organs) with a typical Salsoloid-type Kranz anatomy [29]; the same phenomenon has been observed in Salsola gemmascens of the genus Salsola (Chenopodiaceae), in which cotyledons exhibit C 3 -type photosynthesis, while leaves perform NAD-malic enzyme (NAD-ME) C 4 -type photosynthesis with a Salsoloid-type Kranz anatomy [30]. However, this manner of C 3 -to-C 4 switches has not received much attention, and there were no follow-up studies after their discovery. Investigations into the developmental genetic mechanisms controlling the different photosynthetic types in cotyledons and leaves or other photosynthetic organs will illuminate the genetic regulatory network of C 4 syndrome.
With the development of new technologies, more studies have begun to analyze C 4 syndrome at the systems biology level. To understand C 4 formation, mesophyll cells and bundle sheath cells in the leaf blade of maize were used as a model system for C 4 differentiation. The cells were separated and analyzed for transcriptional changes by microarray analysis and nextgeneration sequencing [31][32][33]; 21% of genes were differentially expressed between mesophyll cells and bundle sheath cells [32]. Recently, John et al. sequenced RNA isolated from the mesophyll cells and bundle sheath cells of Setaria viridis and found the significant convergence of cell-specific gene expression in S. viridis and maize [34]. Such studies deepen our understanding of C 4 syndrome. In 2011, two research groups used mRNA-Seq analysis of closely related C 3 and C 4 species for which gene expression is altered, and these groups identified genes associated with the C 4 pathway [35,36]. Up to 603 and 3582 transcripts differed in abundance between C 3 and C 4 leaves in the genera Cleome and Flaveria, respectively. While these two experiments were designed to identify C 4 -related transcriptomic gene expression changes, it is difficult to tell if the observed variation in transcript abundance was associated with differences between the species or C 4 photosynthesis. However, the comparative transcriptomics of C 3 cotyledons and C 4 assimilating shoots in this study will identify more C 4 -specific genes than before because there is no genetic variation between these different species; more importantly, the natural developmental C 3 -to-C 4 transition may give clues concerning its master switch.

Plant Growth and Harvesting
Seeds of Haloxylon ammodendron were provided by the Turpan Eremophyte Botanic Garden, Chinese Academy of Sciences in Turpan, Xinjiang, China (http://english.egi.cas.cn/rs/sr/tdbg/). The seeds were incubated and germinated on moist filter paper in Petri dishes. After germination, seedlings were planted in sand in a greenhouse maintained at 30/20°C day/night, 70% relative humidity and a photoperiod of 12 h light/12 h dark under a light intensity of 1000 µE m -2 s -1 . Cotyledons were fully expanded after 2 days growing in sand and sampled for all analyses. Assimilating shoots were collected from plants at 10 days of age. For mRNA-Seq, samples were taken from 10-15 individual plants during the middle of the light period, immediately frozen in liquid nitrogen, and stored at-80°C until use.

Light Microscopy
Samples of fully expanded cotyledons and assimilating shoots were fixed in Formalin-acetic acid-alcohol (FAA) for 24 hours, dehydrated through an alcohol series, cleared with xylene, and embedded in paraffin wax. Cross-sections were obtained using a microtome. For light microscopy, semi-thin sections were stained with safranin O solution and studied under DIC microscope (BX51, Olympus, Japan) equipped with an LM Digital Camera (DP70, Olympus).

Stable Carbon Isotope Analysis
Stable carbon isotope ratios ( 13 C/ 12 C) were quantified in cotyledons and dried assimilating shoots from plants grown in the greenhouse. Then, 1-2 cm segments of the middle regions of fully expanded cotyledons or assimilating shoots were collected. All samples were oven-dried at 65°C for 48 h to a constant weight.
The measurements of stable carbon isotope ratios were carried out at the Chinese Academy of Forestry's Stable Isotope Laboratory (Beijing, China) using a Flash EA1112 HT elemental analyzer (Thermo Scientific) coupled with a Delta V advantage isotope ratio mass spectrometer (Thermo Scientific). Stable carbon isotope ratios were expressed as δ 13 C (‰), calculated as follows: where R sample and R standard are the 13 C/ 12 C ratios for an individual sample and the reference standard (Pee Dee Belemnite), respectively.

RNA Preparation and Sequencing
Total RNA was prepared with TRIzol according to the manufacturer's instructions (Invitrogen Life Technologies, Shanghai, China). Following extraction, total RNA was purified using the RNeasy Mini Kit from Qiagen (Shanghai, China), including on-column DNase digestion (Qiagen, Shanghai, China). Purified RNA was checked for integrity and quality using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The cDNA library was constructed for sequencing as described in the Illumina TruSeq RNA sample preparation v2 guide (Catalog # RS-930-1021). Sequencing was performed on an Illumina HiSeq 2000.

Mapping and Quantification of the Sequence Reads
We filtered and examined the quality of the raw sequence reads as described by Xu et al. [37]. The first 10 bases in each read were trimmed off before the mapping process.
Mapping and Quantification of the Sequence Reads from Cotyledons and Leaves of Arabidopsis. Clean reads were mapped onto the latest Arabidopsis thaliana genome assembly (http://www.phytozome.net/arabidopsis.php) using the Bowtie2 (http://bowtie-bio.sourceforge. net/bowtie2/index.shtml) [38]. The best hit of each read with a maximum of three nucleotide mismatches was used. For each Arabidopsis Genome Initiative (AGI) code, the number of matching reads was counted, and the raw digital gene expression counts were normalized using the RPKM (Reads Per Kilobase per Million mapped reads) method [39,40].
Mapping and Quantification of the Sequence Reads from the Cotyledons and Assimilating Shoots of H. ammodendron. Clean reads were mapped onto the coding sequences of the latest A. thaliana genome assembly (http://www.phytozome.net/arabidopsis.php) using BLAT [41]. Alignments were performed in the protein space, and the best hit for each read was retained. For each Arabidopsis Genome Initiative (AGI) code, the number of matching reads was counted, and the raw digital gene expression counts were normalized using the RPKM method [39,40].
The differential expression between samples was statistically accessed by the R/Bioconductor package edgeR [42]. Genes with FDR 0.001 and |log 2 (Fold change)|!1 were considered significant.

Overrepresentation Analysis
To identify functional categories with significant differences between cotyledons and assimilating shoots, we performed an overrepresentation analysis using the GO Term Enrichment of AmiGO (http://amigo.geneontology.org/amigo) [43]. We used all detected transcripts in cotyledons and assimilating shoots as the background set and TAIR as the filter.

Photosynthetic Features of Assimilating Shoots and Cotyledons
Haloxylon species have an unusual photosynthetic apparatus. The true leaves are reduced, and the young annual cylindrical shoots (assimilating shoots) are the main photosynthetic tissues. Fifteen years ago, Pyankov et al. [29] discovered that two Haloxylon species, H. aphyllum and H. persicum, have a C 4 type of photosynthesis in assimilating shoots with Kranz anatomy, whereas leaf-like cotyledons lack Kranz-anatomy and incorporate CO 2 via C 3 photosynthesis, as observed through analyses of stable carbon isotope ratios, anatomy, primary photosynthetic products, and activities of carbon metabolism enzymes. We examined the anatomy assimilating shoots and cotyledons in H. ammodendron and confirmed that H. ammodendron, similarly to the other two Haloxylon species, uses different types of photosynthesis in assimilating shoots and cotyledons, as shown in Fig. 1. The cotyledons have no Kranz-type anatomy and several layers of mesophyll cells around only a few vascular bundles (Fig. 1A). The assimilating shoots have Salsoloid-type Kranz anatomy with two continuous layers of chlorenchyma (a layer of palisade mesophyll cells and an inner layer of bundle sheath cells) on the periphery and waterstorage parenchyma in the center (Fig. 1B). The main vascular bundle occupies the central position, and only the small, peripheral vascular bundles are in contact with bundle sheath cells (Fig. 1B).

Major Transcriptional Changes
To identify differences in transcript abundance related to C 4 syndrome, the transcriptomes of H. ammodendron assimilating shoots and cotyledons were compared. cDNA libraries of assimilating shoots and cotyledons were constructed and sequenced using the Illumina HiSeq 2000 platform, resulting in 30,287,044 and 36,971,687 reads, respectively, with a mean read length of 101 nucleotides. After trimming adapters and filtering out low-quality reads, 29,558,368 reads from assimilating shoots and 36,070,605 reads from cotyledons were retained for further analysis. Clean reads were mapped onto the coding sequences of the latest A. thaliana genome assembly (http://www.phytozome.net/arabidopsis.php) using BLAT [41]. 27037741 reads (~75.0%) from cotyledons and 23188916 reads (~78.5%) from assimilating shoots could be mapped onto the Arabidopsis transcriptome.
mRNA-Seq analysis comparing the transcriptomes of H. ammodendron assimilating shoots and cotyledons yielded 2959 differentially expressed genes [FDR 0.001 and abs (|log 2 (Fold change)|!1)], with 1852 and 1107 more abundant transcripts in assimilating shoots and cotyledons, respectively (see Table A in S1 File). To test whether these differentially expressed transcripts are enriched in functional categories, we performed overrepresentation analysis using the GO Term Enrichment of AmiGO. The significantly up-regulated transcripts in assimilating shoots were enriched in several fundamental biological process categories, including methylation, cytokinesis, DNA replication, cell wall organization or biogenesis, biosynthetic process, anatomical structure morphogenesis, signaling pathway and developmental process (see Table B in S1 File). Down-regulated GO categories are less abundant than up-regulated ones, mainly in response to endogenous stimulus (see Table C in S1 File). Because we used TAIR as a filter and Arabidopsis is a typical C 3 plant, no functional C 4 class was detected.

C 4 Cycle Genes Were Up-regulated in Assimilating Shoots
Transcript analysis of known C 4 genes showed that all genes necessary for the core C 4 cycle of NADP-ME type plants were significantly up-regulated in assimilating shoots compared with cotyledons (Table 1). Among these genes, NADP-ME was most significantly different, with a 12.9-(At5g11670), 12.7-(At2g19900) or 9.3-(At1g79750) fold higher transcript abundance in assimilating shoots. The second biggest difference came from PEPC, with a 12.1-(At2g42600) or 10.5-(At1g53310) fold up-regulation. Additionally, transcripts encoding AspAT (At4g31990) and the chloroplastidic MDH (At5g58330) were up-regulated 8.5-and 2.1-fold, respectively. Transcripts encoding PPDK (At4g15530) were increased 2.7-fold.
All genes required for the NAD-ME type of C 4 photosynthesis were also up-regulated in assimilating shoots compared with cotyledons ( Table 1). The transcripts encoding NAD-ME (At4g00570) were up-regulated 1.5-fold, and those encoding AlaAT (At1g17290) and mtNAD-MDH (At1g53240) were up-regulated 3.8-and 1.5-fold, respectively.
Transcripts encoding the regulatory factors PEPC kinase (PEPC-K) increased 11-fold, which was significant (Table 1). Additionally, the transport proteins required for C 4 photosynthesis were significantly up-regulated (see Table 1), such as pyruvate sodium symporter (BASS2), phosphoenolpyruvate/phosphate translocators (PPT1 and PPT2) and chloroplast dicarboxylate transporters (DiT1 and DiT2). Photorespiratory Genes Were Down-regulated in Assimilating Shoots A major advantage of the C 4 pathway is the reduction in photorespiration because high CO 2 concentration around Rubisco in bundle sheath cells effectively suppresses photorespiration. Detailed analysis of gene expression of all photorespiration genes showed that the transcripts of nearly all genes related to photorespiration were lower in abundance in assimilating shoots than in cotyledons ( Table 2). The genes AtAGT1 (At2g13360), AtGGT1 (At1g23310), AtGLDT1 (At1g11860), AtSHM1 (At4g37930) and AtGLYK (At1g80380), all of which play major roles in photorespiration, were significantly down-regulated (log 2 (Fold change) -1).

Genes Controlling Vein Density Were Up-regulated in Assimilating Shoots
Kranz anatomy is accompanied by high vein density [16], so we assessed the transcript abundance of known genes controlling vein density. All of these genes were up-regulated in assimilating shoots (Table 3). Among them, the transcription factors MONOPTEROS/AUXIN RESPONSE FACTOR5 (MP/ARF5) MP and ATHB8, which are involved in the auxin signal transduction pathway controlling leaf vascular development, were up-regulated 3.9-and 1.9fold, respectively.

Differentially Expressed Transcription Factors Encoding Genes
A total of 248 transcription factor-encoding genes were identified, showing differential expression [FDR 0.001 and abs (|log 2 (Fold change)|!1)] in assimilating shoots compared with cotyledons (Table 4). In total, 107 genes were up-regulated, and 141 genes were repressed in assimilating shoots. To avoid capturing variation in transcript abundance associated with differences between the different developmental stages that do not relate to C 4 photosynthesis, we also sequenced the transcriptomes of Arabidopsis cotyledons and leaves in parallel to minimize the influence of developmental stage-specific effects. Excluding the genes that had the same developmental expression pattern in H. ammodendron and Arabidopsis, there were 96 putative positive regulators and 130 putative negative regulators. We tested the likelihood that SCARECROW/ SHORTROOT components [47] of the C 4 regulatory network were among these identified genes (see Table 5). Scarecrow plays a role in establishing Kranz anatomy in maize leaves, and its mutation results in abnormal Kranz anatomy [48]. SHR was included among the identified genes, and SCR was up-regulated by 1.3-fold in assimilating shoots, although it was not identified.

Discussion
The Haloxylon genus comprises three closely related species, including H. ammodendron (saxaul), H. aphyllum (black saxaul) and H. persicum (white saxaul). H. aphyllum and H. persicum were reported to have different types of photosynthesis in assimilating shoots and cotyledons. Similarly, H. ammodendron assimilating shoots had Salsoloid-type Kranz anatomy (Fig. 1B), indicative of C 4 photosynthesis, while the structure of cotyledons was of the non-Kranz type (Fig. 1A). The developmental transition from a C 3 pathway to a two-celled C 4 pathway in Haloxylon species in nature is a good system with which to study the genetic regulatory network of C 4 syndrome. Stable carbon isotope analysis is used as a screening method to determine the photosynthetic pathway when it is unknown in a species [46]. We measured the δ 13 C values of cotyledons and assimilating shoots. The H. ammodendron cotyledons had a δ 13 C value of-15.58±0.72 ‰ (mean ± SE, n = 3), falling into the range of typical values for C 4 plants of-6 to-19 ‰ [49]. Although this finding is contrary to the anatomical results, it is consistent with a previous report  more negative and closer to that of C 3 due to a mixture of assimilates from C 3 (cotyledons) and C 4 photosynthesis (assimilating shoots).
In this study, we performed deep mRNA-Seq using the Illumina HiSeq 2000 platform to analyze the transcriptomes of H. ammodendron assimilating shoots and cotyledons, which exhibit different types of photosynthesis. Using Arabidopsis as the reference genome, 2959 differentially expressed genes [FDR 0.001 and abs (|log 2 (Fold change)|!1)] were identified, with 1852 and 1107 transcripts being more abundant in assimilating shoots and cotyledons, respectively (see Table A in S1 File).
It was reported that the assimilating shoots of H. aphyllum and H. persicum mainly perform NADP-ME-type C 4 shuttling and a small fraction of NAD-ME-type C 4 shuttling, as shown by photosynthetic enzyme activity analysis and immunoblot analysis [29]. The mRNA-Seq analysis presented here confirmed this discovery and further showed that up-regulation occurs at the level of transcript abundance. All genes necessary for the core C 4 cycle of NADP-ME type plants were significantly up-regulated, and all genes required for the NAD-ME type of C 4 photosynthesis were also up-regulated, but to a lesser extent, in assimilating shoots compared with cotyledons (Table 1). This suggests that NADP-ME-type C 4 photosynthesis is predominant, and NAD-ME-type C 4 photosynthesis makes a small contribution to photosynthetic CO 2 fixation. Likewise, nearly all genes encoding photorespiratory proteins had lower steady-state transcriptional levels in assimilating shoots (Table 2).
We identified a list of positive and negative regulators of C 4 syndrome (Table 4). Although not all of these transcription factors are known to be components of the C 4 regulatory network, our observations suggest that at least a subset of these factors are very likely involved in vein density and the development of Kranz anatomy. Within the list, the transcription factor ATHB8 (Table 4) and the upstream gene MP (Table 3) in the auxin signal transduction pathway controlling leaf vascular development were both up-regulated. This evidence is supportive of a role for at least a subset cohort in vein density. Very recently, Wang et al. (2013) performed a genome-wide comparative analysis of developmental trajectories in Kranz (foliar leaf blade) and non-Kranz (husk leaf sheath) leaves of the C 4 plant maize to look for regulators of Kranz anatomy and identified 48 putative positive regulators, of which, 40 genes were assigned to 38 Arabidopsis orthologous genes [47]. Although maize has a Panicoid-type (classical NADP-ME type) Kranz anatomy and Haloxylon species have a Salsoloid-type Kranz anatomy, there remains high overlap between the list of positive regulators in each study. Within our list of putative positive regulators of C 4 syndrome, 11 of 38 Arabidopsis transcription factor-encoding genes were also significantly up-regulated in H. ammodendron assimilating shoots. Another two genes, AT3G54220 and AT3G13960, were increased by 1.3-and 1.9-fold (Table 6). This high overlap ratio confirms that the methods we used identified Kranz anatomy regulators.
The differences in transcript abundance between the different photosynthetic organs of H. ammodendron, cotyledons and assimilating shoots, reflect the development of C 4 syndrome. Previous studies suggested C 4 cycle genes, trans-factors and even cis-elements were recruited from ancestral C 3 plants [50][51][52][53], but the mechanism behind is not clear. With the same genomic context of the cotyledons and assimilating shoots, this natural and developmental transition from C 3 to C 4 would be an ideal model system to study the molecular mechanism of recruitment, if whole genome sequence and genetic transformation are available.