Improving the efficiency of Rubisco by resurrecting its ancestors in the family Solanaceae.

Plants and photosynthetic organisms have a remarkably inefficient enzyme named Rubisco that fixes atmospheric CO2 into organic compounds. Understanding how Rubisco has evolved in response to past climate change is important for attempts to adjust plants to future conditions. In this study, we developed a computational workflow to assemble de novo both large and small subunits of Rubisco enzymes from transcriptomics data. Next, we predicted sequences for ancestral Rubiscos of the (nightshade) family Solanaceae and characterized their kinetics after coexpressing them in Escherichia coli. Predicted ancestors of C3 Rubiscos were identified that have superior kinetics and excellent potential to help plants adapt to anthropogenic climate change. Our findings also advance understanding of the evolution of Rubisco's catalytic traits.

INTRODUCTION

Rubisco (ribulose-1,5-bisphosphate carboxylase/oxygenase; EC 4.1.1.39) catalyzes the first step of the reductive pentose phosphate cycle by fixing CO2 into ribulose-1,5-bisphosphate (RuBP) (). The catalytic mechanism of Rubisco first arose more than 2.5 billion years ago before the Great Oxidation Event at a time when there was no need to distinguish CO2 from oxygen (O2) (, ). As the O2 level rose, evolution resulted in an increase in Rubisco’s specificity for CO2, but the enzyme could no longer eliminate its oxygenase activity, which leads to a counterproductive process called photorespiration and lowers the photosynthetic efficiency (). In addition, Rubisco is a slow enzyme with a typical turnover number (kcat) of about 2 to 5 s−1 in terrestrial plants, necessitating investment of immense plant resources to produce Rubisco in abundance (). Since Rubisco is a major bottleneck in photosynthesis, understanding how its kinetics evolved in response to changing CO2 and O2 levels is crucial to improving its catalysis in crops (–).

Form I Rubiscos, found in most oxygenic photosynthetic organisms such as cyanobacteria, algae, and plants, are most adapted to aerobic environments and use eight small (S) subunits to stabilize four homodimers of large (L) subunits as hexadecameric L8S8 complexes (, ). In plants and most algae, the L8S8 Rubisco is assembled with the L subunit encoded from a single rbcL gene located in the chloroplast genome and the S subunits produced from the RBCS multigene family in the nucleus and imported into the chloroplast. Considerable progress has been made to engineer Rubisco with superior kinetics into plants by modifying the L subunit (, , ), the S subunit (–), or both subunits simultaneously (–). However, biogenesis of L8S8 complexes in the chloroplast stroma of algae and plants is an elaborate process and involves the chaperonins and multiple chaperones (–). Consequently, evolutionarily distinct foreign Rubisco subunits are poorly compatible with the host chaperones, leading to either no or insufficient production of functional enzymes (, ).

Identifying closely related Rubisco enzymes with superior kinetics is therefore a priority to improve photosynthesis in plants (–). Biochemical analyses of Rubisco from a wide variety of species indicate that Rubisco enzymes with greatly varying kinetic traits exist in nature (–). Periodic reductions in atmospheric CO2 concentrations starting at ~30 million years (Ma) ago have triggered convergent evolution of a CO2-concentrating mechanism (CCM) called C4 photosynthesis in multiple plant families (). A typical Rubisco in a C4 plant has a lower affinity for CO2 and a higher kcat compared to that found in a C3 plant, which has no CCM (, , ). Because of the rapidly increasing atmospheric CO2 levels in the past 200 years, the Rubisco enzymes in C3 plants are likely no longer optimized to the current and future CO2 levels. Although carbon fixation in C3 plants would increase at higher CO2 levels, the increase would be limited by the relatively low kcat of their Rubiscos. Biochemical models predicted that installing selected C4 Rubiscos in C3 plants could improve photosynthesis by more than 25% (, ). Previous attempts to capture kinetic signatures of C4 Rubiscos were mostly performed through evolutionary analyses of the L subunits, with limited success (–, , , ). Despite multiple lines of evidence showing the influence of both subunits on catalysis (, , –), it is still challenging to carry out large-scale phylogenetic analyses of the S subunits in plants due to the lack of available sequences except in a relatively small number of model species.

In this study, we focus on deep phylogenetic analyses of both Rubisco subunits to understand the evolution of C3 Rubiscos in the family Solanaceae. We used the family Solanaceae because any Rubisco modified from a Solanaceous enzyme can be readily expressed in Escherichia coli for characterization of its kinetic properties (, ) and then introduced into a model Solanaceous plant, Nicotiana tabacum (tobacco), for subsequent investigation of its performance in plants (). We developed a computationally efficient workflow to assemble Rubisco sequences de novo from transcriptomics data generated with next-generation sequencing technologies. Data from the workflow markedly expanded the known sequences of both subunits and allowed us to predict sequences at multiple ancestral nodes within the Solanaceae from phylogenetic analyses. We resurrected these predicted ancestral Rubisco enzymes using a recently developed E. coli expression system (, ). Many of these enzymes have kcat values similar to those from C4 Rubiscos and exhibit significantly higher catalytic efficiency than both C3 and C4 Rubiscos. We hypothesize that some of these ancestors could predate the emergence of C4 photosynthesis in several other families and illustrate the evolutionary mechanism of C3 Rubisco through past climate changes. These ancestral Rubisco enzymes appear to be particularly promising candidates to improve photosynthesis in C3 plants.

RESULTS

De novo assembly of Rubisco sequences

We began with Sequence Read Archives (SRAs) containing raw sequences from Solanaceous species at the National Center for Biotechnology Information (NCBI) public repository (www.ncbi.nlm.nih.gov/sra), which were previously generated with next-generation sequencing technologies. Trinity is one of most frequently used bioinformatic programs for de novo assembly of transcript sequences from SRA files (, ). A typical SRA file’s size is several gigabytes with millions of reads derived from thousands of transcripts. As a result, using entire SRA files for de novo assembly is computationally intensive. Since our targets include sequences only from the two Rubisco subunits, we first extracted relevant reads using the BBMap program (https://sourceforge.net/projects/bbmap/) (Fig. 1A). Next, for each set of reads extracted from an SRA, we performed de novo assembly of the Rubisco transcripts with Trinity () under different configurations. Generation of chimeric sequences with de novo assembly is inevitable especially when multiple paralogs with high-sequence homology are present. Thus, we used the known transcript sequences of the S subunits from several model Solanaceae species, such as tobacco, tomato (Solanum lycopersicum), and pepper (Capsicum annuum), as benchmarks to evaluate the accuracy of the assemblies. We found that most of the assemblies were chimeras due to pervasive overlaps among the rbcS paralogs. Thus, we implemented two sequential clean-up steps to identify and remove potential chimeras: (i) Chimeras with overlaps shorter than the read length can be readily recognized from the gaps in their read coverages of starting bases, and (ii) chimeras having long overlaps were found to be assembled much less frequently than the authentic transcripts over multiple Trinity runs and were excluded from the final assemblies (Fig. 1A). We also removed assemblies with extremely low read coverages since they are unlikely to be physiologically important. We tested the workflow with multiple SRAs from each model species and found that we were able to remove all chimeric sequences reliably, although some authentic rbcS transcripts from tobacco were not assembled even after multiple Trinity runs.

Fig. 1.

De novo assembly of Rubisco transcripts from RNA-seq data.

(A) The workflow processes one SRA at a time by downloading it with the SRA Toolkit, extracting the reads aligned to Rubisco subunit sequences with the BBMap program, assembling them de novo with the Trinity program (), and removing potential chimeric assemblies in two clean-up steps. Chimeras with gaps in read coverages of starting bases were identified and removed in the first clean-up step. More potential chimeras with long overlaps with other assemblies were removed in the second clean-up step. The steps automated with Python scripts (available at https://github.com/myattlin/de-novo-assembly) are indicated with green arrows. (B) The numbers of unique L and S subunit protein sequences in Solanaceae assembled in this study and previously available. The numbers of SRAs and assemblies before and after the two clean-up steps are summarized in table S1.

We automated most of the de novo assembly workflow starting from fetching each SRA file from the online repository up to generating images of read coverages used in the first clean-up step with Python scripts that can be executed in Windows Subsystem for Linux (Fig. 1A). Our approach is computationally efficient and can assemble Rubisco sequences from dozens of SRAs a day simply with a modern personal computer equipped with the Windows 10 operating system and high-speed internet. We assembled sequences from 119 publicly available SRA files to obtain 44 unique L subunits and 134 unique S subunits from 15 Solanaceae genera (Fig. 1B and table S1). Trinity was able to assemble complete L subunit sequences from most SRAs, although the data were typically generated from samples enriched with nuclear transcripts. Few chimeras were generated in the assemblies of L subunit sequences, requiring only minimal post-assembly quality control.

Since species belonging to the Solanum and Nicotiana genera were overrepresented in the publicly available sequences, we aimed to expand the number of sequences from a more diverse range of genera from the Solanaceae, with a particular focus on those genera that diverged early in the family’s evolution such as Fabiana, Browallia, Schizanthus, and Vestia, as well as those that emerged from the common ancestor of Solanum and Nicotiana such as Anthocercis, Nicandra, and Jaborosa. We performed additional RNA sequencing (RNA-seq) experiments on complementary DNAs (cDNAs) enriched with S subunit sequences using leaf samples from those seven additional genera and added the sequences for 14 S subunits (table S1).

Predicting ancestral Rubisco sequences

Next, we applied two widely used methods for phylogenetic inference, namely, Bayesian inference and maximum likelihood, with the newly expanded protein sequences of L and S subunits from Solanaceae generated both from mining existing sequences and from the additional RNA-seq experiments (Fig. 2, A and B, and figs. S1 to S4). Since the Rubisco subunits are extremely conserved and not suitable for deriving phylogenetic history, we placed constraints at all major nodes that are consistent with the consensus Solanaceae phylogeny (). Figure 2A also displays three nodes within the family with fossil-calibrated divergence time points and the historical CO2 levels estimated for a similar time frame showing periodic reductions in the CO2 levels that presumably gave rise to C4 photosynthesis in many other families (, ). We named eight ancestral nodes (for example, CaWi for the clade including Capsicum, Lycianthes, Physalis, and Withania genera; SoJa for the clade including Solanum and genera) and separated them into four colored groups based on the similarity among the predicted residue substitutions (Fig. 2, A and B, and table S2). Both Bayesian inference and maximum likelihood generally produced similar predictions, from which we derived 20 and 23 highly probable L and S subunit sequences, respectively, at these nodes, giving rise to 98 predicted ancestral Rubiscos for further characterization (Fig. 2C and Table 1). Compared to the tobacco subunits, the ancestral L and S subunits have up to 12 and 11 mutations, respectively. Notably, the L subunits contain fewer changes than the S subunits except for the Sofa and SoCe ancestors. All three Nico L subunits and four of six Sola and SoDa L subunits are identical to extant Solanaceae L subunits, while only 1 of 23 ancestral S subunits, SoNi2, is found in the extant sequences (Table 1). These findings suggest that the evolution of C3 Rubiscos in response to the climate change in the past 30 Ma has been driven more by changes in the S subunits than in the L subunits.

Fig. 2.

Prediction of ancestral Rubisco sequences in Solanaceae.

(A) A simplified phylogenetic tree for Solanaceous L subunits obtained from Bayesian inference. The fossil-calibrated divergent times for three ancestral nodes () and the names of eight ancestral nodes in four color groups selected in this study are indicated. The inset displays the history of atmospheric CO2 levels estimated from sea surface pH (), with the arrows indicating periodic CO2 reductions that likely resulted in evolution of C4 photosynthesis in several other families. ppm, parts per million. (B) A simplified phylogenetic tree for Solanaceous S subunits obtained from Bayesian inference. The names of eight ancestral nodes in four color groups are indicated. (C) Summary of L and S subunits and Rubiscos predicted for different ancestral nodes of Solanaceae. The residue substitutions in each predicted ancestral subunit are listed in Table 1.

Table 1.

Summary of residue substitutions in the L and S subunits of 98 predicted ancestral Rubisco enzymes.

The identities of extant subunits with the same sequences are also listed.

Predicted ancestral L subunits		Predicted ancestral S subunits		Number of ancestral Rubiscos
Name	Residue substitutions compared to the L subunit of tobacco	Name	Residue substitutions compared to the S-T2 subunit of tobacco
Nico1	L225I K429Q (Sola2 L, Nicotiana acuminata L)	Nico1	N8G V30I E88Q	36
		Nico2	I7Y N8G V30I E88Q
		Nico3	I7Y N8G V30I E88G
		Nico4	I7Y N8G V30I N55H E88G
Nico2	V145I L225I K429Q (Nicotiana undulata L)	SoNi1	K9M V30I E88G
		SoNi2	K9M E23D R28K V30I E88G (Lyciumbarbarum RBCS1)
		SoNi3	K9M V30I E88Q
		SoNi4	N8G K9M V30I E88Q
Nico3	K429Q (Nicotiana tomentosiformis L)	SoNi5	V30I E88Q
		SoNi6	N8G K9M E23D R28K V30I E88Q (Sola2 S)
		SoNi7	N8G E23D R28K V30I E88Q
		SoNi8	N8G E23D R28K V30I K57R E88Q
SoCe1	V145I L225I K429Q C449S V466R A470E V472MV474T	SoCe1	N8G K9M S22T E23D R28K V30IN36K N56H E88Q Q96N	20
SoCe2	V145I L225I K429Q E443D C449S V466R A470EV472M V474T
		SoCe2	N8G S22T E23D R28K V30I N36KN56H E88Q Q96N
Sofa1	V91I V145I L225I K429Q E443D C449S V466RA470E V472M V474T
		SoCe3	N8G S22T E23D R28K V30I K35NN36K N56H E88Q Q96N
Sofa2	V91I V145I L225I V354I K429Q E443D C449SV466R A470E V472M V474T
		SoCe4	N8G S22T E23D R28K V30I K35NN36K N56H K57R E88Q Q96N
Sofa3	V91I V145I L225I V354I K429Q E443D C449SV466R A470E V472M V474T K477GEKK
Sola1	L225I (Przewalskia tangutica L)	Sola1	K9M E23D R28K V30I E88Q	24
Sola2	L225I K429Q (Nico1 L, N. acuminata L)
		Sola2	N8G K9M E23D R28K V30I E88Q (SoNi6 S)
SoDa1	Y226F
SoDa2	Y226F S279T Q439R (Solanum pennellii L)	Sola3	N8G K9M E23D R28K V30I K57RE88Q
SoDa3	None (Atropa belladonna L, Nicotiana sylvestris L)
		SoJa1	K9M E23D R28K V30I K57R E88Q
SoDa4	Y226F S279T
CaWi1	V145I (Salpichroa origanifolia L)	CaWi1	K9M E23D R28K V30I K35R A85NE88Q	18
CaWi2	V145I S279T
CaWi3	V145I L219C	CaWi2	K9M E23D R28K V30I K35R K57RA85N E88Q
CaWi4	V145I L219C E443Q
CaWi5	V145I S279T Q439R C449S	CaWi3	K9M E23D R28K V30I K35R N36SK57R A85N E88Q
CaWi6	V145I L219C E443Q C449S

Ancestral Rubiscos are more efficient

We produced the 98 predicted ancestral Rubisco enzymes of Solanaceae using two expression plasmids that had been previously adapted to produce tobacco Rubisco in E. coli by coexpressing essential chaperonins and chaperones (, ). We screened the RuBP carboxylation activities of these enzymes at a saturating [CO2] using their soluble E. coli extracts. None of the residue substitutions led to a total loss of activity as all samples displayed robust carboxylation activities. Their activities, when normalized with the Rubisco active sites, ranged from about 65 to 128% of the control sample expressing tobacco wild-type (WT) L and S-T2 subunits, with more than half of the predicted ancestors having similar or higher carboxylation rates (Fig. 3). We tested multiple sequences for both L and S subunits at each node because of the nature of ambiguity associated with predicting ancestral sequences and biases arising from incomplete data, which represented only 36 and 22 genera of 92 known genera in Solanaceae for L and S subunits, respectively. As a result, we observed different catalytic rates for the predicted ancestral Rubiscos at each node, likely due to differences in either the S subunits or the L subunits. For example, the Nico ancestors with Nico2, Nico3, and Nico4 S subunits displayed markedly lower carboxylation rates than those with Nico1 S subunits regardless of the L subunits. Among the Sola ancestors, those with Sola1 and Sola2 L subunits have consistently higher carboxylation rates than those with SoDa1 to SoDa4 L subunits (Fig. 3).

Fig. 3.

Initial screening of RuBP carboxylation rates from the predicted ancestral Rubiscos.

The RuBP carboxylation rates were measured at a saturating [CO2] of 108 μM at 25°C under N2 and normalized to the numbers of active sites. Each bar in the chart shows the ratio of the mean of two technical replicates from each sample to that from the tobacco Rubisco with S-T2 subunit expressed in E. coli. Carboxylation kinetics at 25°C were measured for samples marked with * or ** (Figs. 4 and 5). Native PAGE analysis was carried out for samples marked with † (Fig. 6). Carboxylation kinetics at 30°C and SC/O at 25°C were measured for samples marked with ** (Table 2 and Fig. 7). WT, wild type.

Fig. 4.

Carboxylation kinetics of the predicted ancestral Rubiscos of Solanaceae.

(A) A scatterplot for KM,air versus kcat at 25°C. (B) A scatterplot for kcat/KM,air versus kcat at 25°C. RuBP carboxylation rates were measured for 38 predicted ancestors, three tobacco Rubiscos with different S subunits, and seven native Rubiscos from leaf tissues at six [CO2]s, and KM,air and kcat were obtained from nonlinear least square fitting to the classical Michaelis-Menton equation. The means from (n = 3) three E. coli soluble extracts or three leaf soluble extracts from each sample were plotted. The identities of native Rubiscos are as follows: Nb = Nicotiana benthamiana, Np = Nicandra physalodes, Nt = Nicotiana tabacum (Petit Havana), Ph = Petunia hybrida, Sl = Solanum lycopersicum (M28), Ss = Solanum sarrachoides, and St = Solanum tuberosum (Russett Burbank). The SDs and P values compared to the tobacco enzyme are summarized in table S3.

Fig. 5.

Comparison of Rubisco’s RuBP carboxylation kinetics.

Catalytic turnover numbers (kcat), Michaelis-Menten constants for CO2 in air (KM,air), and catalytic efficiencies (kcat/KM,air) at 25°C reported in the literature for Rubiscos from C3 plants, C4 plants (), and those measured in this study from Solanaceous plants and predicted ancestral Rubiscos expressed from E. coli are plotted.

Fig. 6.

Native PAGE analysis of Rubisco complexes in the soluble extracts of tobacco leaf tissue and E. coli cultures.

The immunoblot was performed with an antibody that recognizes form IB Rubisco.

Table 2.

Summary of RuBP carboxylation kinetics at 25° and 30°C for six representative ancestral Rubiscos predicted for different Solanaceae nodes and WT tobacco enzymes with different S subunits.

Means ± SD (P values) of kcat, KM,air, and kcat/KM,air obtained from three E. coli or leaf soluble extracts (n = 3) for each sample are shown. The P values compared to the measurements from the tobacco enzyme with L and S-S1 subunits were determined with two-tailed heteroscedastic t tests. ΔHa values are in kJ−1 mol−1.

Rubisco sample	kcat (s−1)			KM,air (μM)			kcat/KM,air (μM−1 s−1)
	25°C	30°C	ΔHa	25°C	30°C	ΔHa	25 °C	30 °C	ΔHa
Native (Nicotiana tabacum)	3.4 ± 0.2(0.760)	5.1 ± 0.2(0.214)	60.0	18.8 ± 0.7(0.150)	24.0 ± 1.0(0.102)	36.6	0.183 ± 0.004(0.006)	0.214 ± 0.017(0.785)	23.8
Nt-L + Nt-S-S1(reference)	3.5 ± 0.3	4.9 ± 0.2	50.2	17.0 ± 1.4	22.6 ± 1.5	42.6	0.206 ± 0.006	0.219 ± 0.012	9.3
Nt-L + Nt-S-T1	2.4 ± 0.2(0.008)	3.8 ± 0.2(0.001)	67.2	16.0 ± 1.6(0.440)	22.8 ± 1.7(0.902)	53.1	0.151 ± 0.012(0.007)	0.167 ± 0.019(0.035)	14.5
Nt-L + Nt-S-T2	3.4 ± 0.2(0.563)	4.9 ± 0.2(0.889)	54.9	17.7 ± 2.1(0.761)	22.5 ± 0.6(0.772)	35.6	0.193 ± 0.020(0.563)	0.217 ± 0.003(0.984)	18.1
#1 Nico1 L + Nico1 S	4.1 ± 0.2(0.052)	5.6 ± 0.6(0.162)	47.2	18.3 ± 1.3(0.333)	23.5 ± 3.2(0.688)	37.6	0.225 ± 0.004(0.013)	0.241 ± 0.021(0.176)	10.3
#5 Nico2 L + Nico1 S	4.4 ± 0.1(0.030)	5.9 ± 0.3(0.013)	46.2	18.9 ± 0.5(0.141)	22.9 ± 1.6(0.835)	28.4	0.231 ± 0.002(0.012)	0.261 ± 0.007(<0.001)	18.0
#18 Nico1 L + SoNi6 S	4.2 ± 0.0(0.049)	6.0 ± 0.2(0.005)	53.5	18.4 ± 0.6(0.243)	24.3 ± 1.6(0.198)	42.2	0.230 ± 0.005(0.007)	0.248 ± 0.013(0.035)	11.5
#37 Sofa1 L + SoCe1 S	3.7 ± 0.1(0.299)	5.6 ± 0.3(0.034)	61.4	17.7 ± 2.5(0.711)	24.2 ± 1.0(0.101)	46.5	0.214 ± 0.022(0.632)	0.233 ± 0.004(0.002)	13.3
#49 Sola1 L + Sola1 S	4.1 ± 0.1(0.040)	5.8 ± 0.5(0.082)	50.3	19.2 ± 1.4(0.141)	22.4 ± 2.7(0.885)	23.3	0.217 ± 0.011(0.218)	0.260 ± 0.009(0.003)	27.1
#80 CaWi2 L + CaWi2 S	3.7 ± 0.2(0.591)	5.6 ± 0.3(0.025)	61.0	17.2 ± 0.6(0.878)	22.6 ± 0.5(0.997)	41.3	0.218 ± 0.008(0.120)	0.248 ± 0.006(0.001)	19.7

Fig. 7.

CO2/O2 specificity factors (SC/O) of the predicted ancestral Rubiscos of Solanaceae.

(A) The specificity factors were measured at three [CO2]/[O2] ratios at 25°C, and the means and SDs of five or six (n) technical replicates are plotted. The P values compared to the measurements from the tobacco enzyme with L and S-S1 subunits were determined with two-tailed heteroscedastic t tests. (B) Comparison of SC/O at 25°C reported in the literature for Rubiscos from C3 plants, C4 plants (), and those measured in this study for predicted ancestral Rubiscos expressed from E. coli.

Since one of our main interests was to identify Rubisco enzymes with improved catalysis, we selected 38 predicted ancestors, 34 of which displayed higher RuBP carboxylation activities in the initial screening, for measurement of their RuBP carboxylation rates at six different [CO2] levels under air at 25°C along with native Rubisco extracted from leaf tissues of seven Solanaceae species and three E. coli control samples expressing tobacco WT L and either S-S1, S-T1, or S-T2 subunits. The kcat values obtained from these measurements are consistent with their carboxylation activities at the saturating [CO2] (Figs. 3 and 4 and table S3). Several enzymes assembled with Nico L + Nico/SoNi S or Sola L + Sola S subunits have kcat values that are substantially higher than those from the controls and similar to the reported kcat values of typical C4 Rubiscos (Fig. 5) (, ). All the ancestors displayed a similar range of Michaelis constants for CO2 (KM,air) as the control samples. Notably, there appears to be positive correlation between the catalytic efficiencies (kcat/KM,air) and kcat with many of the ancestors with high kcat also having elevated catalytic efficiencies (Figs. 4B and 5 and table S3). Several of those predicted ancestors with extant L subunits such as Nico1, Nico2, Sola1, and Sola2 L subunits (#1, #5, #18, #19, #23, #61, #62, and #67 in table S3) display higher kcat and carboxylation efficiency than the extant Solanaceae and C3 Rubiscos in general (Table 1 and Fig. 5). Hence, S subunits likely play crucial roles in improving the kinetics of these ancestral enzymes.

Just as in our previous study (), the tobacco L + S-T1 Rubisco produced from E. coli displayed a markedly lower kcat likely due to the non-optimal E. coli environment for its assembly. Native polyacrylamide gel electrophoresis (PAGE) analysis of 11 predicted ancestors with both high and low catalytic rates from each of the four ancestral nodes shows that most had similar migration as the tobacco leaf control and L + S-S1 or L + S-T2 enzyme produced in E. coli (Fig. 6). Only ancestor #2 with Nico1 L and Nico2 S subunits and the tobacco L + S-T1 Rubisco migrated at a slightly slower rate. Both the Nico2 S subunit and the tobacco S-T1 subunit share I7Y mutation, which could explain the reduced mobility of the ancestor #2. This did not lead to poor carboxylation catalysis for the ancestor #2 as in the tobacco L + S-T1 Rubisco (table S3). A recent study on Arabidopsis Rubisco expressed in E. coli found that incomplete N-terminal processing of its L subunit led to about 20% lower kcat (). The status of the N-terminal processing of the L subunits in the enzymes expressed from E. coli in our study is not known, but we did not observe any negative impact on the kcat for the Rubiscos expressed from E. coli except for the enzyme with the tobacco S-T1 subunit.

Next, we measured the RuBP carboxylation rates at 30°C for six representative ancestors and the same control samples. Both kcat and KM,air values of all samples were higher at 30° than at 25°C as expected (Table 2). All six ancestors displayed similar or higher activation energies (ΔHa) for kcat/KM,air than the reference WT L + S-S1 control, indicating that their catalysis potentially has a higher optimal temperature. This is not unexpected since these enzymes should be adapted to a hotter climate associated with elevated CO2 more than 20 Ma.

C4 Rubiscos typically have lower CO2/O2 specificity factors (SC/O) compared to C3 versions (, , ). Since many ancestors predicted here have similar kcat as C4 Rubiscos, we tested whether they are also associated with similar SC/O as C4 enzymes. We partially purified six representative ancestral enzymes and measured their SC/O at 25°C. Unexpectedly, the SC/O values of five ancestors are statistically similar to that of the tobacco WT L + S-S1 control. Only one predicted ancestor (#80 CaWi2 L + CaWi2 S) and the tobacco WT L + S-T2 sample had somewhat lower SC/O (Fig. 7A). Comparison to the previously reported SC/O values of C3 and C4 enzymes also indicates that these six ancestors were able to distinguish CO2 from O2 as efficiently as the C3 enzymes (Fig. 7B).

DISCUSSION

In this study, we overcome the lack of available Rubisco sequences, especially for the S subunits, with de novo assembly from transcriptomics data. Our workflow is computationally efficient and capable of removing most, if not all, chimeric assemblies and can generally be applied to any gene of interest. We have identified errors in several NCBI records mostly generated from early periods when DNA sequencing was tedious and had low accuracy (data S1).

The ancestral Rubiscos of Solanaceae predicted in this study appear to be robust, thermally stable, and represent excellent candidates for evolutionary studies. We found several enzymes with higher kcat and efficiency in each of the four ancestral groups, indicating that all these enzymes probably evolved at higher CO2 levels. The best enzymes were identified among Nico and Sola ancestral groups, potentially due to higher accuracy in their predicted sequences enabled by the overrepresentation of extant Solanum and Nicotiana sequences used in our phylogenetic analyses. Despite the relatively small numbers of residue substitutions with no apparent alteration in their overall polarity or electrostatic properties, the subtle mutations in many of these predicted ancestors were able to capture important kinetic traits likely had by the actual ancestors. Notably, most of the predicted ancestors have more mutations in the S subunits than in the L subunits, although the S subunits are only one-fourth the size of the L subunits and are not directly involved in catalysis. A recent study found that the kinetics of potato Rubisco expressed in tobacco were significantly affected by the identity of the S subunit (). This is consistent with our findings that show that many of the predicted ancestors have extant L subunits and yet are able to perform the catalysis more efficiently than the extant enzymes, indicating that the ancestral S subunits in them likely influence the kinetics positively. However, none of the predicted ancestors with enhanced carboxylation abilities contains either of the two unique amino acid residues identified in the S subunit of the potato Rubisco with higher kcat and efficiency (). This highlights the difficulty of predicting the key residues that might control the kinetic properties and the importance of considering both subunits simultaneously to optimize the assembly and overall rigor of the enzyme.

Residue substitutions at 145, 219, 225, 279, 439, and 449 in the L subunits of our predicted ancestors were previously identified to be positively selected during the evolution of Rubiscos in plants (), and the L225I substitution in most of our predicted ancestral L subunits is consistent with the I225L substitution previously found to be associated with the evolution of C3 Rubiscos (). It is not unexpected that none of the substitutions in our predicted ancestors was found to be involved in the transition from C3 to C4 photosynthesis () since C4 photosynthesis is not present in Solanaceae. Because the residues altered in both subunits of the ancestors are not directly associated with those at the active site, it is challenging to decipher how these residue substitutions were able to influence the kinetic properties without further structural studies.

In some families with both C3 and C4 photosynthesis, the C3 Rubiscos have lower SC/O than the average SC/O of typical C3 Rubiscos, which likely facilitated the evolution of C4 photosynthesis in those families (). In contrast, the ancestral C3 Rubiscos of Solanaceae predicted here have similar SC/O as typical C3 Rubiscos. Recent structural analyses indicated a correlation between SC/O and positively charged cavities close to the active site (). On the basis of the residue substitutions, most of the predicted Solanaceae ancestors are expected to have similar electrostatic profiles as typical C3 Rubiscos. Nevertheless, our findings support the hypothesis that the catalytic behavior of C3 Rubiscos in ancient plants before the emergence of C4 photosynthesis may be more similar to the present day C4 Rubiscos in having higher kcat. The evolution of C4 photosynthesis likely shifted their Rubiscos’ SC/O and affinity for CO2 lower, while the enzymes remaining in C3 plants shifted their kcat lower during their adaptation to decreasing CO2 levels. A previous study on the C3 and C4 L subunits in Flaveria species identified residue 309 as the catalytic switch, which is specific to the Flaveria species and incompatible with the tobacco L subunit background (). Multiple ancestral L and S subunits of Solanaceae characterized in this study were able to achieve the high catalytic rates of C4 enzymes without sacrificing affinity for CO2. It is also noteworthy that these ancestral subunits are highly similar to the tobacco sequences and are expected to be compatible with the Rubisco assembly system of tobacco chloroplasts. Our approach can be applied to study Rubiscos in other families of higher plants, especially the ones that include C4 members, to investigate whether their ancestral Rubiscos display comparable features.

Higher catalytic efficiency of Rubisco is beneficial not only for growth but also for water and nitrogen use efficiency in plants. The ancestral Rubiscos predicted in this study also appear adapted to hotter and drier environments based on their catalysis at a higher temperature and SC/O values that are similar to the current C3 Rubiscos. The next step will be to introduce these ancestral Rubiscos into plants and assess their performance. Although the technology to replace both Rubisco subunits was recently reported for tobacco (), transformed plants must be able to produce sufficient amount of Rubisco to take advantage of improved kinetics. Emerging technologies such as targeted base editing of chloroplast genes () should expand the engineering of Rubisco to other plants where generation of stable chloroplast transformation is not available. The procedure in this study can be a blueprint to identify superior Rubiscos in other families to eventually enhance carbon fixation in agricultural crops such as rice and wheat.

MATERIALS AND METHODS

De novo assembly of sequences encoding Rubisco subunits

Each SRA file was downloaded with fastq-dump 2.8.0 program available from SRA Toolkit (https://hpc.nih.gov/apps/sratoolkit.html). The SRA file’s reads aligned to sequences encoding Rubisco L or S subunits were selected with BBMap 38.22-1 program (by Bushnell B; https://sourceforge.net/projects/bbmap/) using the DNA sequences encoding tobacco L subunit or the mature S subunit S1 as references in “vslow” and “local” modes and “maxindel” set to 100. Next, the paired reads in the fastq file exported by BBMap were separated into two fastq files with BBMap’s bbsplitpairs scripts. Reads in the two fastq files were then assembled de novo by Trinity 2.8.5 three separate times as follows: (i) --KMER_SIZE 32; (ii) stringent setting, which includes “--min_kmer_cov 4 –min_glue 4 –min_iso_ratio 0.2 –glue_factor 0.2 –jaccard_clip”; and (iii) both –KMER_SIZE 32 and stringent setting. If there were more than 10,000 reads in each fastq file, the first 5000 reads extracted by seqtk 1.3-r106 program (https://github.com/lh3/seqtk) were assembled in two more Trinity runs with –KMER_SIZE 32 with or without the stringent setting. The read coverages of starting bases for coding sequences were then obtained for assemblies that covered at least 90% of the reference sequences with alignment scores greater than 350 using BBMap scripts with “perfectmode” and “startcov = t” settings. We automated the above process with Python scripts (Fig. 1A) and executed them in Windows Subsystem for Linux from a shell script file, which can be supplied with multiple SRA IDs for high-throughput assembly. The Python scripts and example shell scripts are available at https://github.com/myattlin/de-novo-assembly/. The scripts were written for the paired-end format of SRA files, although they can be adapted for single-end format with slight modifications. The automated process wrote SRA IDs, reference files used in BBMap, assembled sequences, sequences encoding the L and S subunits of Rubisco, and locations for the read coverage files of all assemblies to a csv file. In addition, it also saved read coverage files and PNG format images of read coverage profiles for the assemblies. In the first clean-up step, the read coverage images were visually inspected for gaps to remove chimeric assemblies. In the second clean-up step, assemblies generated for each species were compared against one another for the presence of long overlaps, and those that have long overlaps and were assembled at lower frequencies were removed.

RNA-seq of partial rbcS transcripts

The seeds for Browallia viscosa (Bv), Nicandra physalodes (Np), Schizanthus coccineus (Sc), Schizanthus grahamii (Sg), and Vestia lyciodes (Vl) were obtained from www.plant-world-seeds.com, and Anthocercis littorea (Al), Fabiana imbricata (Fi), and Jaborosa sativa (Js) were obtained from http://b-and-t-world-seeds.com. DNA oligonucleotides were synthesized by Integrated DNA Technologies Inc. (Coralville, IA, USA) and are included in data S4. An Invitrogen PureLink RNA mini kit (Thermo Fisher Scientific Inc.) was used to prepare RNA samples from leaf tissues of plants grown under 100 photosynthetically active radiation (μmol/m2 per second) with a 16-hour photoperiod in Lambert LM-111 all-purpose mix. Invitrogen SuperScript III First-Strand Synthesis Supermix (Thermo Fisher Scientific Inc.) was used to synthesize cDNA with the Not I–dT-R oligonucleotide according to the manufacturer’s instructions. Partial rbcS transcripts were amplified from each cDNA sample by Phusion high-fidelity DNA polymerase with Not I–Adpr-R and Mau BI–SSU-D-F oligonucleotides, and ~650–base pair (bp) amplicons were extracted from agarose gels with an EZ-10 spin-column polymerase chain reaction (PCR) product purification kit (Thermo Fisher Scientific Inc.). Bv, Np, Sc, Sg, and Vl samples were fragmented with Covaris E220 followed by reparation and adenylation of ends and adapter ligation with a TruSeq DNA PCR-Free kit (Illumina Inc.) before they were pooled and sequenced with NextSeq 550 (Illumina Inc.) in 2 × 150–bp runs. Np, Al, Fi, and Js samples were fragmented and indexed with a Nextera DNA library prep kit (Illumina Inc.) and sequenced with MiSeq nano (Illumina Inc.) in 2 × 250–bp runs. All new RNA-seq data are available at www.ncbi.nlm.nih.gov/sra/PRJNA759940.

Predicting ancestral rubisco sequences

Multiple sequence alignments of the Rubisco L and S subunits were performed with Clustal Omega 1.2.4 (). Bayesian inference was performed separately with MrBayes 3.2.7a () using the amino acid sequences of the L and S subunits with the following parameters: lset nst = mixed rates = invgamma, prset aamodelpr = mixed, mcmc ngen = 600,000 for L subunits or 800,000 for S subunits, temp = 0.06 for L subunits or 0.04 for S subunits, startparams = reset, and starttree = random. The topology was fixed at multiple nodes based on the reported consensus tree (), and the probabilities of the ancestral states at those nodes were generated with the setting “report applyto= (1) ancstates = yes.” The average SDs of split frequencies from Metropolis-coupled Markov chain Monte Carlo sampling bottomed at about 0.02. The nexus files executed in MrBayes are available as data S2. The ancestral states were also estimated with RAxML 8.2 () with PROTGAMMAAUTO for model configuration, autoMRE for rapid bootstrapping with automatic criteria, “-g” option with a constraint tree file to ensure the topology remained consistent with the established tree (), and “-f A” setting with the resulting best tree rooted with FigTree program v1.4.3 (https://github.com/rambaut/figtree/releases). The phylogenies of L and S subunits reached convergence after 650 and 750 bootstrap replicates, respectively. The input sequence files, constraint files, and the command lines for the RAxML runs are available as data S3. From the predicted probabilities at each residue position of eight selected nodes (table S2), we selected 98 combinations of ancestral L and S subunits (Table 1).

Expressing the predicted ancestral rubiscos in E. coli

DNA oligonucleotides were purchased from Integrated DNA Technologies Inc. (Coralville, IA, USA), and their sequences are available in the Supplementary Materials. Phusion high-fidelity DNA polymerase, FastDigest restriction enzymes, and T4 DNA ligase were purchased from Thermo Fisher Scientific Inc. and used to amplify, digest, and ligate DNA fragments. Mau BI site was inserted before T7P-lacO-RBS-Nt-rbcL operon by amplifying the operon with Mlu I–Age I–Mau BI–for and BJFEseqR oligonucleotides from BJFE-T7P-lacO-RBC-Nt-rbcL plasmid (), which was then digested with Mlu I and Not I and ligated into the Mlu I and Not I sites of a holding vector to obtain pHD-T7P-NtL vector. Next, T7P-lacO-RBC-NtrbcL operon digested from pHD-T7P-NtL with Age I was ligated into the Age I site of pAtC60αβ/C20 () vector to obtain pET-AtC60AB20-T7P-NtL-v2 vector. The L subunit gene was separated into three fragments based on the two internal restriction sites: Bam HI at residue 155 and Nde I at residue 387. The mutations in the predicted ancestral L subunits (Table 1) were introduced with overlapping PCRs by corresponding oligonucleotides and accumulated in each of the three fragments, which were then simultaneously ligated into Mau BI and Not I sites of pET-AtC60AB20-T7P-NtL-v2 vector to generate the final expression vectors. The tobacco S subunit T2 gene was separated into two fragments at Eco RI restriction site located at residues 43 to 44 and used as the template to generate the predicted ancestral S subunits (Table 1). Substitutions at residues 23, 28, 30, 85, 88, and 96 were achieved by overlapping PCRs, while the remaining substitutions were generated with a Q5 site-directed mutagenesis kit (New England Biolabs) with the corresponding oligonucleotides. The mutations accumulated in each of the two fragments were combined by ligation into Nco I and Not I sites of pCDF-NtXT2R1AtR2NtB2 vector () to obtain the final expression vectors. The sequence of each ligated DNA in the expression vectors was confirmed by Sanger sequencing. The pET-AtC60AB20-T7P-NtL-v2 and pCDF-NtXT2R1AtR2NtB2 vectors were cotransformed into BL21*(DE3) E. coli, and each Rubisco sample was expressed from the E. coli culture grown in ZYP-5052 autoinduction medium as described previously ().

Enzyme kinetics of the predicted ancestral Rubiscos

Soluble extracts from 6-ml E. coli cultures lysed in 400 μl of 50 mM tris-HCl (pH 8), 10 mM MgCl2, 1 mM EDTA, 20 mM NaHCO3, 2 mM dithiothreitol (DTT), and Pierce protease inhibitor minitablet (Thermo Fisher Scientific Inc.) were used to measure RuBP carboxylation activities of the Rubisco samples. For leaf extracts, about 5 cm2 of leaf tissue each suspended in 500 μl of 100 mM Bicine-NaOH (pH 7.9), 5 mM MgCl2, 1 mM EDTA, 5 mM ε-aminocaproic acid, 2 mM benzamidine, 50 mM 2-mercaptoethanol, protease inhibitor cocktail, 1 mM phenylmethanesulfonyl fluoride, 5% (w/v) poly(ethylene glycol) 4000, 10 mM NaHCO3, and 10 mM DTT was crushed in a 2-ml Wheaton homogenizer for about 1 min on ice, and insoluble materials were removed by centrifugation at 16,000 rcf at 4°C for 5 min. Each supernatant of leaf extracts was then applied to a 2-ml Zeba spin desalting column with 40,000 molecular weight cutoff preequilibrated with 100 mM Bicine-NaOH (pH 8), 20 mM MgCl2, 1 mM EDTA, 1 mM benzamidine, 1 mM ε-aminocaproic acid, 1 mM KH2PO4, 2% (w/v) poly(ethylene glycol) 4000, 20 mM NaHCO3, 10 mM DTT, and each eluate following centrifugation at 1000 rcf at 4°C for 2 min was incubated at 23°C for 30 min for full activation of Rubisco active sites. RuBP carboxylation experiments were performed as described previously with NaH14CO3 solutions with different concentrations and specific activities, such that 14C activities of acid-stable compounds in the vials following the termination of the reactions gave a similar range of values (). For initial screening of the 98 predicted ancestral enzymes, RuBP carboxylation activities were measured in vials equilibrated with N2 gas at 25°C and 108 μM [CO2], and 14C fixed to stable organic compounds was counted with Tri-Carb 2810TR Scintillation counter (PerkinElmer). The same Rubisco samples were used for quantification of Rubisco active sites on the same day with 14C-carboxyarabinitol bisphosphate (CABP) bound to each sample as described previously (). The specific activity of 14C CABP was precalibrated with a soluble extract from spinach leaf tissue, where the Rubisco concentration was determined from an immunoblot along with a commercial spinach RbcL standard (Agrisera, part no. AS01 017S) using a polyclonal antibody against wheat Rubisco (). To measure kcat and KM,air, the RuBP carboxylation activities of E. coli soluble extracts with 38 predicted ancestral Rubiscos and three tobacco Rubiscos and soluble extracts from tobacco leaf tissue were measured at six different [CO2] concentrations ranging from 5.5 to 90 μM at pH 8 in vials equilibrated with CO2-free air at 25°C, and the Rubisco active sites were subsequently quantified with 14C CABP. kcat and KM,air were obtained from nonlinear least square fitting to the classical Michaelis-Menton equation as described previously (). Three biological replicates were performed for each sample from three separate E. coli cultures or leaf extracts. The same measurements were repeated at 30°C for six predicted ancestral Rubisco samples and the same control samples of tobacco Rubiscos.

Specificity factors of the predicted ancestral rubiscos

CO2/O2 specificity factors (SC/O) of six predicted ancestral Rubiscos and tobacco Rubiscos were measured with partially purified Rubisco samples. First, E. coli pellets from 1.5- to 2-liter cultures were each resuspended in ~20 ml of extraction buffer [25 mM triethanolamine (pH 8), 5 mM MgCl2, 0.5 mM EDTA, 1 mM KH2PO4, 1 mM benzamidine, 5 mM ε-aminocaproic acid, 10 mM 2-mercaptoethanol, 5 mM NaHCO3, 2 mM DTT, and 1 mM phenylmethylsulfonyl fluoride] and sonicated with eight 10-s pulses over 5 min at 4°C. Insoluble materials were separated with centrifugation at 35,000g at 4°C for 30 min. The supernatant was applied to a 5-ml HiTrap Q HP anion exchange column (GE Healthcare) connected to the ÄKTA P-900 Fast Protein Liquid Chromatography System equipped with an Inv-907 valve and a Frac-950 fraction collector and equilibrated with Q buffer [25 mM triethanolamine (pH 8), 5 mM MgCl2, 0.5 mM EDTA, 1 mM benzamidine, 1 mM ε-aminocaproic acid, 5 mM NaHCO3, 2 mM DTT, and 12.5% (v/v) glycerol]. NaCl in the buffer applied to the column was then increased from 0 to 0.5 M over 75 ml of volume at 2 ml min−1, and the eluents were collected in 2-ml fractions. The Rubisco-containing fractions were identified by bound 14C CABP, concentrated to ~500 to 700 μl with Amicon Ultra-15 centrifugal filter units, and stored at −80°C before use. Rubisco was also purified with the 5-ml HiTrap Q HP column from ~500 cm2 of tobacco leaf tissue broken in ~200 ml of extraction buffer in a blender, precipitated with PEG at a final concentration of ~20% (w/v), and resuspended in ~10 ml of Q buffer. Total protein concentration in the samples was estimated with Bradford assays. The Rubisco purified from tobacco leaf tissue represented about 90% of the total soluble protein, while the Rubisco samples from E. coli represented about 25 to 30% of the total soluble protein.

The SC/O values were calculated with the formula (RuBP carboxylated / RuBP oxygenated) / ([CO2] / [O2]) after measuring RuBP carboxylated at three different ratios of [CO2] / [O2] (). The amount of RuBP oxygenated was derived from the total RuBP consumed in each experiment. After ~25 nmol of RuBP was entirely catalyzed by ~140 pmol of Rubisco active sites at three [CO2] concentrations in each reaction vial equilibrated with CO2-free air at 25°C, the 14C fixed to stable organic compounds was counted. Each reaction was also repeated in a second vial with 2 min of additional incubation period to ensure that all RuBP was consumed in both measurements. In addition, each reaction was repeated in a vial equilibrated with N2 gas, from which the total amount to RuBP consumed in each vial was obtained, since all RuBP was carboxylated in these vials.

Native PAGE and immunoblot

Soluble extracts were prepared from either E. coli cultures or tobacco leaf tissue in the same procedure as in the determination of Rubisco kinetics as described above. The total soluble protein concentrations were determined with Bradford assays, and 4 μg of total soluble proteins from each E. coli extract or 0.1 μg from tobacco leaf extract was mixed with the loading buffer made up of 50 mM bis-tris (pH 7.2), 50 mM NaCl, 0.001% Ponceau S, and 10% glycerol. The electrophoresis was carried out in an Invitrogen 3 to 15% bis-tris protein gel from Thermo Fisher Scientific with 50 mM bis-tris and 50 mM tricine (pH 6.8) anode buffer and 0.002% Coomassie Brilliant Blue G250, 50 mM bis-tris, and 50 mM tricine (pH 6.8) cathode buffer at 150 V and 4°C for 30 min followed by 250 V for 60 min. The samples were then transferred to a polyvinylidene difluoride membrane with 0.45-μm pore size in 25 mM tris, 192 mM glycine, and 20% methanol at 100 V and 4°C for 1 hour. The membrane was blocked with 5% milk in TBST (tris-buffered saline with Tween 20) buffer [20 mM tris (pH 7.5), 150 mM NaCl, and 0.1% Tween 20] at 23°C for 1 hour, incubated with an antibody against Rubisco (from P.J. Andralojc from Rothamsted Research, raised in a rabbit) in 5% milk in TBST buffer at 4°C overnight, and detected with horseradish peroxidase–conjugated secondary antibody in 2.5% milk in TBST buffer at 23°C for 1 hour. The chemiluminescent signals from enhanced chemiluminesence substrate were captured with a ChemiDoc MP imaging system from Bio-Rad.