Data from: Convergence, stability, and thermal adaptation in the rubisco large subunit in plants
Data files
Jun 04, 2026 version files 42.11 MB
-
leung_et_al_rubisco_evolution.zip
42.05 MB
-
README.md
65.87 KB
Abstract
Enzymes are adapted to perform optimally in different thermal regimes that would otherwise alter kinetics and stability. Whether adaptive evolution in the photosynthetic enzyme ribulose-1,5-bisphosphate carboxylase/oxygenase (rubisco) also compensates for thermal variation remains uncertain. We examined molecular evolution and modelled the change in folding free energy (ΔΔG, where negative values indicate stabilization) of the rubisco large subunit (RbcL) in four phylogenetically distant plant genera: wood ferns (Dryopteris), sea lavenders (Limonium), pines (Pinus), and viburnums (Viburnum). Using codon evolutionary models in each genus, we observed widespread positive selection and parallel substitution in the catalytic α/β barrel domain. Species with warmer growing seasons had derived amino acids with stronger hydrogen bond contributions to ΔΔG. Protein structure-based modelling showed that the hydrogen bond contribution to stability tracked the growing season temperature of species carrying the derived amino acid. Stronger hydrogen bonds were offset by weaker contributions from hydrophobic solvation interactions, such that total ΔΔG showed no relationship with growing season temperature. In Viburnum, the strength of positive selection differed among biomes, with cold temperate and cloud forest clades showing stronger positive selection. These patterns are consistent with environmental tuning of non-covalent interactions within the enzyme. However, modest effect sizes indicate that other components of the rubisco holoenzyme likely also contribute to its thermal evolution.
DOI:10.5061/dryad.h44j0zq1z
This leung_et_al_rubisco_evolution.zip archive contains the scripts, configuration files, and key analysis outputs for the analyses described in the manuscript. To keep the archive small, several large intermediate files (per-occurrence GBIF and bioclim records, raw NCBI GenBank data, IQ-TREE/CSUBST state files, rendered figure tiffs, and software executables) are not included; all can be regenerated by re-running the relevant scripts on the inputs that are included. The protein structure PDB 8RUC is available from the RCSB Protein Data Bank.
Overview
We investigate positive selection and convergent amino acid evolution in the large subunit of Rubisco (rbcL) across Viburnum (Adoxaceae), with comparisons to Pinus, Limonium, and Dryopteris. The analysis combines phylogenetic selection tests (PAML), convergent substitution detection (csubst), protein stability modelling (FoldX), and species distribution/climate data (GBIF, WorldClim).
Software Requirements
The scripts may run on other versions of the software listed, but below are the versions we used.
| Software | Version | Citation |
|---|---|---|
| R | ≥ 4.2 | R Core Team |
| Python | ≥ 3.9 | |
| PAML (codeml) | 4.10 | Yang (2007) Mol Biol Evol 24: 1586–1591 |
| MUSCLE | 5 | Edgar (2022) Nat Commun 13: 6968 |
| MEGA | 11 | Tamura, Stecher & Kumar (2021) Mol Biol Evol 38: 3022–3027 |
| IQ-TREE | 2 | Minh et al. (2020) Mol Biol Evol 37: 1530–1534 |
| CSUBST | 1.4.0 | Fukushima & Pollock (2023) Nat Ecol Evol 7: 155–170 |
| FoldX | 5.1 | Delgado et al. (2019) Bioinformatics 35: 4168–4169 |
| PyMOL | 3.0.4 | Schrödinger LLC |
| GNU Parallel | 20250322 | Tange (2025) |
R packages
R packages used:
- tidyverse
- here
- readr
- dplyr
- tidyr
- stringr
- purrr
- ape
- phytools
- phangorn
- rgbif
- CoordinateCleaner
- plyr
Key R package citations: ape — Paradis & Schliep (2019) Bioinformatics 35: 526–528; phytools used for Pagel (1994) tests of correlated evolution; rgbif and CoordinateCleaner for GBIF occurrence retrieval and cleaning.
Python packages
Python packages used:
- biopython
- pandas
- numpy
Repository Structure
.
├── download/ # Step 1 — Download rbcL sequences from NCBI
├── align/ # Step 2 — Multiple sequence alignment (MUSCLE)
├── dambe/ # Step 3 — Substitution-saturation analysis (IQ-TREE)
├── prune_tree/ # Step 4 — Prune Landis et al. (2021) species tree
├── paml/ # Step 5 — PAML site models (M0–M8)
├── paml_site_reconstruction/ # Step 6 — PAML ancestral sequence reconstruction
├── paml_cmc/ # Step 7a — Branch-site Clade Model C (CmC)
├── paml_cmd/ # Step 7b — Clade Model D (CmD)
├── paml_cmd_high_K/ # Step 7c — CmD with high kappa constraint
├── paml_cmd_high_w/ # Step 7d — CmD with high omega constraint
├── paml_cmd_high_K_w/ # Step 7e — CmD with both constraints
├── paml_cmc_high_K/ # Step 7f — CmC with high kappa constraint
├── paml_cmc_high_w/ # Step 7g — CmC with high omega constraint
├── paml_cmc_high_K_w/ # Step 7h — CmC with both constraints
├── paml_m2a_rel/ # Step 7i — M2a_rel model
├── csubst/ # Step 8 — Convergent substitution analysis (Viburnum)
├── download_gbif/ # Step 9 — GBIF occurrences + WorldClim climate data
├── landis_biome_range/ # Step 10 — Biome and range data (Landis et al. 2021)
├── analyze/ # Step 11 — Parse outputs, Pagel's tests, FoldX results
├── structure/ # Step 12 — Structural analysis (PDB: 8RUC)
├── foldx_subs/ # Step 13 — FoldX ΔΔG for positively selected sites
├── foldx_species/ # Step 14 — FoldX ΔΔG per substitution on ancestral backgrounds
├── combined_stability/ # Step 15 — Combined stability figures across the 4 genera
│
├── pine_rubisco/ # Pinus rbcL analysis (mirrors Viburnum workflow)
├── limonium_rubisco/ # Limonium rbcL analysis
├── dryopteris_rubisco/ # Dryopteris rbcL analysis
│
└── README.md
Missing and Unavailable Data
The following files are not included in this archive due to file size, licensing, or regenerability:
- WorldClim rasters — Monthly temperature (
wc2.1_30s_tmin,wc2.1_30s_tavg) and CGIAR Global Aridity Index v3 (ai_v3) tiles must be downloaded separately from worldclim.org and the CGIAR-CSI database, and placed indownload_gbif/worldclim/. The same applies to the equivalentworldclim/directories withinpine_rubisco/,limonium_rubisco/, anddryopteris_rubisco/. - Per-occurrence GBIF records — Raw occurrence data downloaded by
get_gbif_data.Rare not included. Only per-species summary tables (download_gbif/supplemental/species_sample_counts.csv,*_summary.csvfiles in each genus'sgbif_data/folder, andanalyze/bioclim_summary.csv) are retained. Re-runget_gbif_data.R(requires a GBIF account andrgbifcredentials) to regenerate. - Raw NCBI GenBank files — The raw file
sequence.gbin each genus directory contains the GenBank-formatted download. These are large; re-runget_sequence.pywith the relevantaccessions.csvto regenerate. - IQ-TREE and CSUBST state files — Large binary checkpoint and state files (e.g.,
.ckp.gz,.model.gz) are included where present but may be absent for some runs; re-run IQ-TREE or csubst to regenerate. - FoldX executables and rotabase — The FoldX binary and
rotabase.txtmust be obtained from the FoldX website (https://foldxsuite.crg.eu/) and placed in the working directory before running any FoldX steps. A copy ofrotabase.txtis present infoldx_species/for reference. - Rendered figure TIFFs — High-resolution rendered figures are not included. SVG and JPG versions of figures are provided in
analyze/plots/andcombined_stability/. - PDB 8RUC — The spinach Rubisco holoenzyme CIF file (
structure/8ruc-assembly1.cif) is included. The original is publicly available from the RCSB Protein Data Bank (https://www.rcsb.org/structure/8RUC).
Analysis Workflow
Part 1 — Viburnum rbcL Analysis
Step 1. Download sequences — download/
python download/get_sequence.py # Download rbcL sequences from NCBI GenBank
python download/clean_sequence.py # Trim to coding sequence, remove stop codons
Input: download/accessions.csv (GenBank accession numbers from Landis et al. 2021)
Output: download/clean.fasta
File descriptions:
| File | Contents |
|---|---|
accessions.csv |
GenBank accession numbers and associated species names used to download rbcL sequences, taken from Landis et al. (2021) Systematic Biology |
get_sequence.py |
Downloads rbcL sequences from NCBI GenBank using accessions in accessions.csv |
clean_sequence.py |
Trims each sequence to the coding region and removes stop codons |
clean.fasta |
Cleaned, trimmed rbcL sequences in FASTA format, ready for alignment |
compare_supp_tree.R |
Compares the downloaded species list against the Landis et al. supplemental species tree to identify mismatches |
landis_species.csv |
Full species list from Landis et al. (2021) used for comparison |
missing_from_download.csv |
Species present in one dataset but absent from the other |
missing_from_download_check.csv |
Manual follow-up verification of identified mismatches |
metadata.csv |
Additional metadata (e.g., collection locality, sequence length) for downloaded sequences |
sequence.gb |
Raw GenBank-format file from NCBI; large, regenerable via get_sequence.py |
Step 2. Align — align/
run_muscle.bat # MUSCLE v5 alignment
run_disperse.bat # Manual-check helper
Manual inspection of the alignment was performed in MEGA v11 (Tamura et al. 2021); in viburnums the first 31 codons were trimmed (because of partial CDS) and the stop codon removed.
Then convert to PHYLIP:
python align/fasta_to_phylip.py
Output: align/clean_aligned.fasta, align/clean_aligned_checked.fasta
File descriptions:
| File | Contents |
|---|---|
run_muscle.bat |
Windows batch file to run MUSCLE v5 alignment on download/clean.fasta |
run_disperse.bat |
Windows batch file helper for visual inspection during alignment checking |
fasta_to_phylip.py |
Converts the checked FASTA alignment to PHYLIP format for PAML |
clean_aligned.fasta |
Raw MUSCLE v5 alignment output |
clean_aligned_checked.fasta |
Alignment after manual inspection in MEGA (first 31 codons trimmed, stop codon removed) |
clean_aligned_checked.fas |
MEGA-format copy of the checked alignment |
Step 3. Substitution-saturation analysis — dambe/
dambe/run_iqtree.bat # IQ-TREE under a GTR+G model on the rbcL alignment
Used to assess substitution saturation in the rbcL alignment before downstream codon-model analyses. The same procedure is run for each comparative genus in <genus>/dambe/.
File descriptions:
| File | Contents |
|---|---|
run_iqtree.bat |
Windows batch file to run IQ-TREE on the PHYLIP alignment under GTR+G |
clean_aligned.phylip |
rbcL alignment in PHYLIP format (IQ-TREE input) |
clean_aligned.phylip.iqtree |
IQ-TREE run log (model parameters, substitution rates) |
clean_aligned.phylip.treefile |
Maximum-likelihood tree from IQ-TREE |
clean_aligned.phylip.mldist |
ML pairwise distance matrix |
clean_aligned.phylip.bionj, .contree, .uniqueseq.phy |
Auxiliary IQ-TREE output files |
clean_aligned.phylip.ckp.gz, .model.gz |
IQ-TREE checkpoint files (may be absent; regenerate by re-running IQ-TREE) |
Step 4. Prune species tree — prune_tree/
python prune_tree/prune_tree.py
Prunes the Landis et al. (2021) RadSeq + cpDNA phylogeny (radseq_cpdna_stage2.tre) to the taxa with rbcL sequences.
Output: prune_tree/radseq_pruned.tre
File descriptions:
| File | Contents |
|---|---|
prune_tree.py |
Prunes the full Landis et al. (2021) species tree to the subset of taxa with rbcL sequences |
radseq_cpdna_stage2.tre |
Full Landis et al. (2021) RadSeq + cpDNA phylogeny (input) |
radseq_pruned.tre |
Pruned species tree containing only taxa with rbcL sequences (output; used by PAML and csubst) |
clean_aligned.phylip, clean_aligned_paml.phylip |
Alignment copies in PHYLIP format used in downstream PAML runs |
Step 5. PAML site models — paml/
Rscript paml/unroot_tree.R # Unroot tree for PAML
paml/run_paml.bat # Run codeml with site models M0, M1a, M2a, M3, M7, M8
Config: paml/codeml.ctl
Output: paml/mlc (log-likelihoods, parameter estimates), paml/rst (ancestral reconstructions)
Parse results:
Rscript analyze/paml_parse.R # → analyze/paml_site_mlc.csv, paml_site_omega.csv, paml_site_arc.csv
File descriptions:
| File | Contents |
|---|---|
codeml.ctl |
PAML codeml configuration file specifying site models (M0, M1a, M2a, M3, M7, M8), input alignment, and tree |
run_paml.bat |
Windows batch file; unroots tree with unroot_tree.R then calls codeml |
unroot_tree.R |
Unroots the pruned species tree for PAML (PAML requires an unrooted tree for site models) |
mlc |
PAML log-likelihood output: model parameters, κ, ω, BEB posterior probabilities for each site under each model |
rst |
PAML ancestral reconstruction output: marginal probability distributions of amino acids at each internal node |
radseq_pruned.tre, radseq_pruned_unrooted.tre |
Input trees (rooted and unrooted versions) |
clean_aligned.phylip, clean_aligned_paml.phylip |
rbcL alignment in PHYLIP format |
Note on Windows batch files:
run_paml.bat(and the analogous batch files in all otherpaml_*/directories) referenceC:\Program Files\R\R-4.2.2\bin\Rscript.exedirectly. Users running a different R version or installation location will need to update this path, or replace it with the bare commandRscriptif R is on the system PATH.
Step 6. Ancestral sequence reconstruction — paml_site_reconstruction/
paml_site_reconstruction/run_paml_reconstruction.bat # codeml model M8 reconstruction
Output: paml_site_reconstruction/rst — ancestral amino acid probabilities at each node
File descriptions:
| File | Contents |
|---|---|
run_paml_reconstruction.bat |
Windows batch file to run codeml ancestral reconstruction under model M8 |
codeml.ctl |
Configuration file for M8 reconstruction run (sets RateAncestor = 1) |
mlc |
PAML output log for reconstruction run |
rst |
Ancestral amino acid probability distributions at each internal node under M8 (used by FoldX ancestral-background pipeline) |
clean_aligned_paml.phylip |
Input alignment |
radseq_pruned_unrooted.tre |
Input tree |
Step 7. Branch-site clade models — paml_cmc/, paml_cmd/, variants
Runs PAML CmC and CmD models to test for episodic positive selection on convergent lineages. The foreground branches are defined by csubst/branch/branches_of_interest.txt.
paml_cmc/run_paml.bat # Clade Model C
paml_cmd/run_paml.bat # Clade Model D
Parse results:
Rscript paml_cmc/paml_CmC_parse.R
Rscript paml_cmd/paml_CmD_parse.R
Sensitivity analyses with fixed κ and ω starting values are in paml_cmd_high_K/, paml_cmd_high_w/, paml_cmd_high_K_w/, paml_cmc_high_K/, paml_cmc_high_w/, and paml_cmc_high_K_w/. The null model is in paml_m2a_rel/.
Each clade-model directory contains the same set of files:
| File | Contents |
|---|---|
codeml.ctl |
Configuration for the specific clade model; foreground branches marked with #1 |
run_paml.bat |
Windows batch file to run codeml; references the R executable (see Step 5 note) |
mlc |
PAML log-likelihood output for this model and branch partition |
rst |
Ancestral reconstructions under this model |
radseq_pruned_ur_clades.tre |
Unrooted tree with foreground branches labelled (#1) |
clean_aligned.phylip, clean_aligned_paml.phylip |
Alignment copies |
paml_CmC_parse.R / paml_CmD_parse.R |
Parses codeml output and extracts log-likelihoods and parameter estimates |
run_paml.command |
macOS/Linux equivalent of run_paml.bat (present in some directories) |
The high-K/high-w/high-K-w variants set initial values of κ (kappa) and/or ω (omega) to higher starting points to check that the primary results are not sensitive to starting parameter values. paml_m2a_rel/ runs the M2a_rel neutral model as a null for the CmC likelihood ratio test.
Step 8. Convergent substitution analysis — csubst/
Prepare inputs:
python csubst/clean_fasta_csubst.py # Format FASTA for csubst
Rscript csubst/clean_tree.R # Format tree for csubst
Rscript csubst/branch/branches_of_interest.R # Define foreground branches
Run csubst (external; see csubst documentation):
# csubst cb --alignment csubst/clean_aligned_csubst.fasta \
# --rooted_tree csubst/radseq_pruned_csubst.tre \
# --foreground csubst/branch/branches_of_interest.txt \
# --aamodel empirical
Analyse results:
Rscript csubst/analysis/tree.R
Rscript csubst/analysis/get_clades.R
Rscript csubst/analysis/get_convergent_clades.R
Rscript csubst/analysis/merge_site_data.R
File descriptions — csubst/ (root):
| File | Contents |
|---|---|
clean_fasta_csubst.py |
Reformats the checked FASTA alignment to the header format required by csubst |
clean_tree.R |
Reformats the pruned species tree to the Newick format required by csubst |
clean_aligned_csubst.fasta |
Reformatted alignment for csubst input |
clean_aligned_csubst.fasta.iqtree |
IQ-TREE log from the tree csubst uses internally |
clean_aligned_csubst.fasta.treefile |
IQ-TREE tree used by csubst for branch-length estimation |
clean_aligned_csubst.fasta.rate |
Site-specific rate file from IQ-TREE |
radseq_pruned_csubst.tre |
Species tree reformatted for csubst |
radseq_pruned.tre |
Pruned species tree (copy from prune_tree/) |
clean_aligned.fasta, clean_aligned_checked.fasta |
Alignment copies (FASTA format) |
clean_aligned_paml.phylip |
Alignment in PHYLIP format |
csubst_branch.command, csubst_site.command |
Example shell commands for running csubst branch-pair and site analyses |
File descriptions — csubst/branch/:
| File | Contents |
|---|---|
branches_of_interest.R |
Defines the foreground branch groups for csubst (convergent lineages) |
branches_of_interest.txt |
Foreground branch specification file passed to csubst --foreground; also used to label foreground branches in the PAML clade-model runs |
branches_of_interest_.txt |
Draft/working version of the foreground file (retained for reference) |
csubst_b.tsv |
csubst branch-level output: substitution counts and convergence statistics per branch |
csubst_cb_stats.tsv |
csubst convergent-branch-pair statistics: pairwise convergence metrics for each foreground branch pair |
csubst_alignment_aa.fa |
Amino-acid alignment generated internally by csubst |
csubst_alignment_codon.fa |
Codon alignment generated internally by csubst |
csubst_tree.nwk |
Tree with csubst internal node labels in Newick format |
csubst_marginal_branch.txt |
Marginal ancestral state probabilities from csubst |
csubst_instantaneous_rate_matrix.tsv |
Amino-acid substitution rate matrix used by csubst |
branch_id_nodes.csv |
Mapping of csubst branch IDs to tree node identifiers |
File descriptions — csubst/analysis/:
| File | Contents |
|---|---|
tree.R |
Annotates and visualises the csubst tree with convergence information |
get_clades.R |
Extracts clade membership for each branch from the csubst tree |
get_convergent_clades.R |
Filters csubst output to identify branch pairs with significant convergence |
merge_site_data.R |
Merges csubst site-level convergence data with PAML BEB probabilities |
convergent_clades.csv |
Convergent lineage pairs identified by csubst (used in downstream analyses) |
convergent_clades_subset.csv |
Subset of convergent clades meeting additional filtering criteria |
landis_clades.csv |
Clade assignments from Landis et al. (2021) mapped to tree nodes |
branch_id_nodes.csv |
Branch ID to node mapping (copy of csubst/branch/branch_id_nodes.csv) |
Step 9. GBIF occurrences and climate data — download_gbif/
Rscript download_gbif/get_gbif_data.R
Downloads GBIF occurrence records for all Viburnum species, runs coordinate cleaning, and extracts WorldClim monthly temperature data and CGIAR Global Aridity Index values at occurrence points. Requires a GBIF account and rgbif credentials. WorldClim rasters (Fick & Hijmans 2017) must be downloaded separately from worldclim.org and placed in download_gbif/worldclim/ (not deposited here due to size). Aridity index rasters (Zomer et al. 2022, Global Aridity Index v3) must be downloaded separately from the CGIAR-CSI database.
Growing season is defined as months with monthly mean temperature ≥ 5 °C and aridity index > 0.05.
Per-occurrence GBIF and bioclim records are not included in this archive due to size; only the per-species summary tables (download_gbif/gbif_data/*_summary.csv, equivalent files in each comparative-genus gbif_data/ folder, and analyze/bioclim_summary.csv) are retained. Re-run get_gbif_data.R to regenerate the per-occurrence files.
File descriptions:
| File | Contents |
|---|---|
get_gbif_data.R |
Downloads GBIF occurrences, cleans coordinates, extracts WorldClim and aridity index values, and produces per-species climate summaries |
supplemental/GBIF_citation.txt |
Full data citation for the GBIF occurrence download, as required by GBIF's terms of use |
supplemental/species_keys.csv |
GBIF species keys for all Viburnum species queried |
supplemental/taxon_keys.csv |
GBIF taxon keys used to resolve accepted names |
supplemental/species_sample_counts.csv |
Number of GBIF occurrence records retrieved per species after coordinate cleaning |
Step 10. Biome and range data — landis_biome_range/
Rscript landis_biome_range/parse_biome_range.R
Parses biome and geographic range reconstructions from Landis et al. (2021).
File descriptions:
| File | Contents |
|---|---|
parse_biome_range.R |
Parses biome and range NEXUS files from Landis et al. (2021) and outputs tidy CSV tables |
viburnum.biome.n4.nex |
NEXUS file of biome character states for Viburnum species (from Landis et al. 2021) |
viburnum.range.n6.nex |
NEXUS file of geographic range character states (from Landis et al. 2021) |
biomes.csv |
Parsed biome assignments per species |
landis_biomes.csv, landis_ranges.csv |
Biome and range data tables ready for joining with other analyses |
Step 11. Main analysis — analyze/
python analyze/aa_at_sites.py # Extract amino acid states at positively selected sites
Pagel's tests for correlated evolution between amino acid states and ecological variables are recorded in analyze/pagels_test_*.xlsx.
File descriptions:
| File | Contents |
|---|---|
aa_at_sites.py |
Extracts the amino acid at each positively selected site for each taxon from the alignment |
paml_parse.R |
Parses PAML mlc and rst files from paml/; outputs paml_site_mlc.csv, paml_site_omega.csv, and paml_site_arc.csv |
analyze.Rmd |
Main analysis R Markdown: Pagel's tests for correlated evolution, bioclim correlations, figure generation |
paml_site_mlc.csv |
PAML site model log-likelihoods and likelihood ratio test statistics for each model comparison |
paml_site_omega.csv |
BEB posterior probabilities of positive selection (ω > 1) for each site |
paml_site_arc.csv |
Most likely ancestral amino acid at each positively selected site for each node |
paml_site_mlc_unfill.xlsx |
Unformatted version of the PAML LRT table (working file) |
pos_selected_conv_sites.csv |
Sites inferred to be both positively selected (BEB > 0.95) and convergently substituted by csubst |
aa_at_sites.csv |
Amino acid states at positively selected sites across all extant Viburnum taxa |
aa_at_sites_ancestral.csv |
Amino acid states at positively selected sites at all ancestral nodes |
bioclim_summary.csv |
Per-species summary of WorldClim bioclim variables (means and ranges across GBIF occurrences) |
bioclim_aa_summary.csv |
Summary of bioclim values grouped by amino acid state at each selected site |
climond_aa_summary.csv |
Same as bioclim_aa_summary.csv using CliMond climate variables |
lma_aa.csv |
Leaf mass per area (LMA) data joined with amino acid states |
phy_anova.csv, phy_anova_edwards2014.csv, phy_anova_edwards2017.csv |
Phylogenetic ANOVA results for amino acid state vs. ecological variables under different comparative methods |
convergent_subs.xlsx |
Summary table of convergent substitutions identified at positively selected sites |
h_bonds.xlsx |
Hydrogen bond counts at convergent sites from structural analysis |
pglmm_mesophyll.xlsx |
Phylogenetic generalised linear mixed model results for mesophyll-related variables |
pagels_test_cold.xlsx |
Pagel's (1994) test for correlated evolution between amino acid state and cold-temperature adaptation |
pagels_test_cold_reverse.xlsx |
Same test with the direction of correlation reversed |
pagels_test_habit.xlsx |
Pagel's test for correlated evolution with growth habit (not in the manuscript but retained for reference) |
pagels_test_margin.xlsx |
Pagel's test for correlated evolution with leaf margin type (not in the manuscript but retained for reference) |
pagels_test_tropical.xlsx |
Pagel's test for correlated evolution with tropical biome occupancy (not in the manuscript but retained for reference) |
pagels_test_tropical_reverse.xlsx |
Same test with correlation direction reversed (not in the manuscript but retained for reference) |
Rplot001.svg |
Exploratory plot (intermediate output) |
plots/ |
Directory of SVG figures output by analyze.Rmd, including dN/dS plots (dnds_viburnum.svg) and stochastic character mapping plots for each convergent site (simmaps/) |
Step 12. Structural analysis — structure/
Rscript structure/get_polar_residues.R # Identify polar residues at convergent sites
Input: structure/8ruc-assembly1.cif (spinach Rubisco holoenzyme; PDB: 8RUC, Andersson 1996)
Visualisation in PyMOL v3.0.4. PyMOL commands: structure/pymol.txt
File descriptions:
| File | Contents |
|---|---|
get_polar_residues.R |
Identifies polar residues at convergent sites and calculates polarity statistics |
8ruc-assembly1.cif |
Spinach Rubisco holoenzyme structure in mmCIF format (PDB: 8RUC; Andersson 1996); also available from the RCSB PDB |
pymol.txt |
PyMOL command script for colouring and visualising convergent and positively selected sites |
h_bonds_all.xlsx |
Complete hydrogen bond count table for all residues at convergent sites, from structural analysis |
macrostructure.png, macrostructure_unlabelled.png |
Rendered structural images of the Rubisco holoenzyme |
pymol pics/ |
Directory of PyMOL-rendered structural images |
Step 13. FoldX ΔΔG — positively selected sites — foldx_subs/
Computes ΔΔG (change in folding free energy) for each amino acid substitution observed at positively selected sites, using the spinach Rubisco structure as background.
python foldx_subs/mutations.py # Generate mutation list
python foldx_subs/split_mutations.py # Split into 8 parallel batches
python foldx_subs/constraints.py # Generate constraints file
Repair the structure, then run FoldX BuildModel in 8 batches (3 replicates each):
# foldx --command=RepairPDB --pdb=8ruc-assembly-clean.pdb
# foldx --command=BuildModel --pdb=8ruc-assembly-clean_Repair.pdb \
# --mutant-file=individual_list_mutations_1.txt --numberOfRuns=3
Parse output:
Rscript foldx_subs/read_foldx.R # → foldx_subs/ddG_values.csv
File descriptions:
| File | Contents |
|---|---|
mutations.py |
Generates the complete FoldX mutation list for all substitutions at positively selected sites |
split_mutations.py |
Splits the mutation list into 8 batches (individual_list_mutations_1.txt through _8.txt) for parallel FoldX runs |
constraints.py |
Generates the FoldX constraints file (individual_list_constraints.txt) |
aa_at_WT_site.py |
Extracts the wild-type (spinach) amino acid at each site of interest |
read_foldx.R |
Parses FoldX BuildModel output files and computes mean ΔΔG across replicates |
__buildmodel_SUBS.sh |
Shell script to run FoldX BuildModel for all 8 batches |
__echo_progress.sh |
Helper script to report FoldX progress |
individual_list_mutations_1.txt – individual_list_mutations_8.txt |
FoldX mutation input files (one per batch) |
individual_list_mutations.txt |
Combined (unsplit) mutation list |
individual_list_constraints.txt |
FoldX constraints file |
individual_listdraft.txt |
Draft mutation list (working file) |
individual_list.txt |
Alternative mutation list format |
ddG_values.csv |
Final ΔΔG values for each substitution at positively selected sites (spinach Rubisco background); key output for stability analyses |
repaired_pdb/ |
Directory containing the FoldX-repaired spinach Rubisco PDB |
molecules/ |
Directory of FoldX BuildModel intermediate structural output files |
Step 14. FoldX ΔΔG — ancestral backgrounds — foldx_species/
Computes per-substitution ΔΔG on the reconstructed ancestral Rubisco structures of each group (Viburnum, Pinus, Limonium, Dryopteris), allowing assessment of background-dependent stability effects.
python foldx_species/00_extract_ancestor_from_rst.py # Extract ancestral sequences from PAML rst
python foldx_species/01_generate_ancestor_mutlist.py # Generate ancestor→species mutation lists
python foldx_species/03_generate_mutations.py # Generate FoldX input files
Build ancestral structures and run ΔΔG calculations:
bash foldx_species/01_build_ancestor.sh # Build ancestor PDB from spinach structure
bash foldx_species/04_run_ddg.sh # Run FoldX BuildModel (all ancestors, ~species)
bash foldx_species/run_4genera_full.sh # Master wrapper for all 4 genera
Run single-mutation ΔΔG for each site substitution:
python foldx_species/06_generate_single_mutations.py
bash foldx_species/07_run_single_mutations.sh
Combine results:
python foldx_species/05_combine_results.py
python foldx_species/08_combine_single_mutations.py
File descriptions — foldx_species/ (root):
| File | Contents |
|---|---|
00_extract_ancestor_from_rst.py |
Extracts most-likely ancestral sequences from PAML rst files for each genus |
01_generate_ancestor_mutlist.py |
Generates the list of mutations from the inferred root ancestor to each extant species |
03_generate_mutations.py |
Generates per-genus FoldX BuildModel input mutation files |
01_build_ancestor.sh |
Applies the ancestor-to-spinach mutations to the repaired spinach PDB to build each ancestral PDB |
04_run_ddg.sh |
Runs FoldX BuildModel for each species on the appropriate ancestral background |
05_combine_results.py |
Combines FoldX BuildModel output into per-genus ΔΔG results tables |
06_generate_single_mutations.py |
Generates single-mutation FoldX input files for each substitution at each site |
07_run_single_mutations.sh |
Runs single-mutation FoldX BuildModel for each site substitution |
08_combine_single_mutations.py |
Combines single-mutation ΔΔG results across genera |
__repair.sh |
Runs FoldX RepairPDB on each ancestral structure |
run_4genera_full.sh |
Master wrapper that coordinates all 4-genera FoldX runs |
run_4genera_ddg.sh, run_4genera_phase2.sh |
Sub-steps of the 4-genera wrapper |
rotabase.txt |
FoldX rotamer library (required by FoldX; also available from the FoldX website) |
repaired_pdb/ |
Repaired ancestral PDB files |
Per-genus ancestor subdirectories (viburnum_ancestor/, dryopteris_ancestor/, limonium_ancestor/, pine_ancestor/):
| File | Contents |
|---|---|
*_ancestor_wt.pdb |
Reconstructed ancestral Rubisco large subunit PDB built from the spinach structure |
*_ancestor_wt_Repair.pdb |
FoldX-repaired version of the ancestral PDB (input for BuildModel) |
*_root_ancestral_species.fasta |
Ancestral sequence at the root node of this genus (from PAML rst) |
individual_list_mutations.txt |
FoldX mutation list (ancestor → each extant species) |
*_ddg_results.tsv |
Per-species ΔΔG values computed on the ancestral background |
*_single_mutation_ddg.tsv |
Single-mutation ΔΔG values for each substitution at each site |
species_to_sequence.tsv |
Mapping of species names to their reconstructed amino acid sequences (Viburnum only) |
species_mutations/ |
Per-species mutation list files |
logs/, runs/, single_mutations/, single_runs/, molecules/ |
FoldX intermediate output directories |
Part 2 — Comparative rbcL Analysis (Pinus, Limonium, Dryopteris)
Each comparative group has its own directory (pine_rubisco/, limonium_rubisco/, dryopteris_rubisco/) following the same numbered-script workflow as above:
| Script | Step |
|---|---|
0. get_taxa.R |
Retrieve taxon list |
1. get_sequences.py |
Download rbcL sequences from NCBI |
1a. clean_sequences.py |
Clean sequences; was not necessary for Pinus, because GenBank data were complete coding sequences |
4. fasta_to_phylip.py |
Convert alignment to PHYLIP |
5. unroot_tree.R |
Unroot tree for PAML |
6. run_paml.bat |
Run PAML |
7. paml_parse.R |
Parse PAML output |
7a. aa_at_sites.py |
Extract amino acids at sites |
8. run_paml_reconstruction.bat |
Ancestral reconstruction |
10. branches_of_interest.R |
Define csubst foreground |
12. merge_site_data.R |
Merge csubst site data |
13. get_gbif_data.R |
GBIF occurrences |
14. analyze.Rmd |
Genus-specific analysis |
All comparative analyses (PAML site/clade models, csubst convergence, FoldX ddG) are run on the published species trees of each genus, not on rbcL gene trees, for better resolution of clades. Each genus directory contains the same output files as described for the Viburnum workflow above, including per-species PAML outputs (mlc, rst), per-species GBIF and bioclim summaries (gbif_data/), csubst outputs (branch/), and IQ-TREE saturation results (dambe/). Key genus-specific files:
| File | Contents |
|---|---|
accessions.csv |
GenBank rbcL accession numbers for this genus; except for Dryopteris, for which accessions taken from the original paper for multiple genes are in sessa_tree_accession.csv |
clean.fasta |
Cleaned rbcL sequences |
clean_aligned.fasta, clean_aligned_checked.fas |
Alignment files |
clean_aligned_paml.phylip, clean_aligned_csubst.fasta |
PAML- and csubst-formatted alignments |
*_paml.tre, *_csubst.tre |
Published species trees formatted for PAML and csubst |
mlc, rst |
PAML site-model output |
paml_site_mlc.csv, paml_site_omega.csv, paml_site_arc.csv |
Parsed PAML output (same format as Viburnum) |
aa_at_sites.csv, aa_at_sites_ancestral.csv |
Amino acid states at selected sites |
pos_selected_conv_sites.csv |
Sites identified as positively selected and convergent |
convergent_subs.xlsx |
Convergent substitution summary for this genus |
codeml.ctl |
PAML configuration file |
paml_reconstruction/ |
Subdirectory with PAML M8 reconstruction run |
gbif_data/bioclim_summary.csv, gbif_data/bioclim_aa_summary.csv |
Per-species and per-amino-acid climate summaries |
gbif_data/gbif_citation.txt |
GBIF data citation for this genus's download |
gbif_data/species_keys.csv, gbif_data/taxon_keys.csv, gbif_data/species_sample_counts.csv |
GBIF metadata |
pagels_test_*.xlsx |
Pagel's tests for this genus |
bisse.xlsx |
BiSSE model results (where applicable) |
plots/ |
Figures for this genus |
branch/ |
csubst branch outputs (same structure as csubst/branch/) |
dambe/ |
IQ-TREE saturation analysis outputs (same structure as dambe/) |
h_bonds.xlsx |
Hydrogen bond analysis at convergent sites |
Species-tree files specific to each genus:
- Dryopteris:
sessa_tree_accessions.csv,sessa_tree_labels.csv— Taxon mapping for the Sessa et al. species tree;dryopteris_paml.tre,dryopteris_csubst.tre - Pinus:
Pinus_time.tre,Pinus_time_paml.tre,Pinus_time_csubst.tre - Limonium: species tree files as provided by the source publication
Combined stability — combined_stability/
Combines per-substitution ddG results across the four genera and produces the cross-genus stability figures.
File descriptions:
| File | Contents |
|---|---|
cross_genus_analysis.Rmd |
R Markdown that combines FoldX ΔΔG results across the four genera and generates cross-genus stability figures |
cross_genus_analysis.html |
Rendered HTML output of cross_genus_analysis.Rmd |
dated_species_tree.tre |
Time-calibrated species tree grafting all four genera, used for PGLS analyses |
single_mutation_wls_combined.csv |
Weighted least squares (WLS) results for single-mutation ΔΔG across all genera combined; columns: response (energy term), r, r², p_wls, n, label |
single_mutation_wls_ps_only.csv |
Same WLS analysis restricted to positively selected sites |
species_pgls_grafted_tree.csv |
PGLS results for per-species cumulative ΔΔG on the grafted species tree |
species_wls_per_genus.csv |
Per-genus WLS results for per-species ΔΔG |
figure_*.jpg |
Cross-genus stability figures output by cross_genus_analysis.Rmd |
Root-level files
Note, the *.xlsx files are included for convenience but can also be generated as output files from R scripts (see below).
| File | Contents |
|---|---|
rubisco_tables.Rmd |
R Markdown that generates the supplemental data tables (PAML LRT tables, FoldX summary, Pagel's test tables) for the manuscript |
table_PAML_LRTs.xlsx |
Summary PAML likelihood ratio tests (site models) |
table_PAML_CmC_LRTs.xlsx |
CmC branch-site LRTs across all foreground partitions |
table_PAML_CmD_LRTs.xlsx |
CmD branch-site LRTs |
table_CSUBST.xlsx |
csubst convergent substitution summary table |
table_Pagel_tests.xlsx |
Pagel's test results (all traits) |
supp_table_FoldX.xlsx |
FoldX ΔΔG supplemental table |
Key Output Files
| File | Contents |
|---|---|
analyze/paml_site_mlc.csv |
PAML site model log-likelihoods and LRTs |
analyze/paml_site_omega.csv |
BEB posterior probabilities per site |
analyze/paml_site_arc.csv |
Ancestral reconstruction at positively selected sites |
analyze/pos_selected_conv_sites.csv |
Convergent sites under positive selection |
csubst/analysis/convergent_clades.csv |
Convergent lineages from csubst |
foldx_subs/ddG_values.csv |
ΔΔG for substitutions at selected sites (spinach background) |
combined_stability/single_mutation_wls_combined.csv |
WLS results for single-mutation ΔΔG (all genera combined) |
combined_stability/species_pgls_grafted_tree.csv |
PGLS results for per-species ΔΔG |
table_PAML_LRTs.xlsx |
Summary PAML likelihood ratio tests |
table_PAML_CmC_LRTs.xlsx |
CmC branch-site LRTs |
table_PAML_CmD_LRTs.xlsx |
CmD branch-site LRTs |
table_CSUBST.xlsx |
csubst convergent substitution summary |
table_Pagel_tests.xlsx |
Pagel's test results |
supp_table_FoldX.xlsx |
FoldX ΔΔG supplemental table |
Notes
- All R scripts use the
herepackage; set the project root withhere::i_am()as shown in each script. Run from the repository root. - PAML must be run before
analyze/paml_parse.Rand before the FoldX ancestral structure pipeline. - The
csubstforeground branch files (branches_of_interest.txt) are generated bybranches_of_interest.Rafter PAML and csubst have been run. - WorldClim rasters are not deposited. Download 30-arcsec monthly temperature (
wc2.1_30s_tmin,wc2.1_30s_tavg) and CGIAR aridity index (ai_v3) tiles and place them indownload_gbif/worldclim/. - PAML site-model runs use
cleandata = 0because the rbcL alignment is a partial CDS; usingcleandata = 1would discard substantial data. - Windows batch files (
.bat):run_paml.batandrun_paml_reconstruction.batin the PAML directories, and6. run_paml.batin the comparative genus directories, callC:\Program Files\R\R-4.2.2\bin\Rscript.exedirectly. Users with a different R version or installation path will need to update this line, or replace it with the bare commandRscriptif R is on the system PATH. All other file paths within the batch files are relative to the script location (%~dp0). - The CSV output files in
combined_stability/(single_mutation_wls_combined.csv,single_mutation_wls_ps_only.csv,species_pgls_grafted_tree.csv,species_wls_per_genus.csv) are generated directly bycross_genus_analysis.Rmd.
