Cross-species gene redesign leveraging ortholog information and generative modeling

Domain: chemistry
Published: 2026-03-03
Editorial outcome: accept

Abstract
Conventional approaches to heterologous gene expression rely on codon optimization, which is limited to swapping synonymous codons and fails to capture deeper adaptive changes. In contrast, naturally evolved orthologous genes include non-synonymous mutations, insertions, and deletions that confer functional adaptation to different host contexts. Here we present OrthologTransformer, a Transformer-based deep learning model that converts orthologous genes between species by learning from large-scale orthologous gene datasets. The model recapitulates evolutionary differences—from synonymous codon swaps to amino acid-changing mutations and indels—to predict coding sequences optimized for target species while preserving protein function. In extensive tests across diverse bacterial species pairs, the model’s context-aware gene designs more closely resembled native host orthologs, preserved protein functionality, and achieved superior expression yields compared to codon-optimized sequences. As proof of concept, an OrthologTransformer-redesigned PETase expressed in Bacillus subtilis showed robust activity, producing approximately 10-fold more reaction product than the codon-optimized enzyme, and achieving higher expression levels, thereby demonstrating improved enzyme performance via AI-guided gene design.

Introduction
Heterologous gene expression—expressing a gene in a non-native host organism—has become a cornerstone of modern biotechnology, enabling the production of recombinant proteins and enzymes in convenient expression systems. A fundamental challenge in this process is that different organisms exhibit distinct codon usage biases; the frequency with which synonymous codons are used can vary greatly between a gene’s native species and the host species1,2. To overcome this, codon optimization is routinely employed to modify a gene’s coding sequence without altering the amino acid sequence, replacing rare codons with those preferred by the host’s translational machinery. Traditional codon optimization methods—such as algorithms maximizing the codon adaptation index (CAI) or simply swapping out low-frequency codons—have been successful in many cases3–5. However, these conventional approaches have notable limitations: they inherently freeze the protein sequence1,4 and often neglect other critical sequence features that influence protein expression, such as mRNA secondary structure, codon context, regulatory motifs, and translation kinetics6–9. As a result, a sequence optimized purely on codon usage frequencies may still perform suboptimally in the host, and in some instances over-optimization can even reduce protein yield7,10,11.

Orthologous genes—genes in different species that encode the same or similar proteins—provide valuable insights into how nature resolves the challenge of adapting a coding sequence to different organismal contexts12. Even between moderately distant species, orthologs typically accumulate numerous non-synonymous substitutions (amino acid-changing mutations) and sometimes small insertions or deletions (indels) over the course of evolution13. For instance, when comparing the coding sequences of Bacillus thuringiensis genes to their orthologs in Bacillus subtilis, only about 57% of codons encode the same amino acid (25% of codons are identical and about 32% are synonymous substitutions where the amino acid is conserved). The remaining 43% of codons differ in amino acid, and in addition, insertions and deletions account for roughly 4% of positions. Similarly high levels of amino acid divergence are observed for more distantly related species (e.g., Escherichia coli vs B. subtilis orthologs have over 50% non-synonymous codon differences). These observations underscore that nature’s solution to cross-species gene transfer is not limited to swapping codons for synonymous alternatives—it also involves adjusting the protein sequence itself when necessary.

In recent years, data-driven approaches have begun to address some shortcomings of purely frequency-based approaches14,15. For example, deep learning models trained on large collections of synonymous gene variants and their expression levels have captured subtle sequence features (e.g., codon pair context, GC content distribution, mRNA structure propensity) associated with high protein yields. Deep generative models can design synthetic coding sequences that not only match a host’s codon usage bias but also minimize problematic structures or motifs, effectively learning the language of codons for a given organism. Indeed, several studies have demonstrated improved protein expression using neural network–designed gene sequences compared to those optimized by conventional algorithms16. Nevertheless, existing machine-learning-based gene design tools still typically restrict changes to synonymous substitutions, leaving the protein sequence itself unaltered.

In this study, we introduce a problem formulation at the DNA level: orthologous gene conversion framed as a sequence-to-sequence (seq2seq) translation task between species. We develop OrthologTransformer, a general deep learning framework based on a seq2seq architecture17 that expands codon optimization into cross-species gene adaptation for prokaryote. By training on orthologous gene pairs from diverse bacteria, OrthologTransformer learns to introduce nucleotide changes (including amino acid substitutions and indels) in a biologically informed manner, aiming to produce a gene sequence that looks native to the target host while preserving the protein’s function. Thus, OrthologTransformer is trained in a supervised seq2seq manner on naturally occurring ortholog pairs that already embody the balance between host adaptation and function preservation. Consequently, the model does not decide this balance via an explicit objective or hand tuned weights; rather, it reproduces the patterns present in orthologs (synonymous changes, conservative amino acid substitutions, occasional indels). We show that OrthologTransformer-generated sequences more closely resemble true orthologs of the target species and improve heterologous expression outcomes. We then showcase a practical case study: the PETase enzyme, which enables bacteria to degrade PET plastic18. We use OrthologTransformer to convert the PETase gene from its native bacterium (Ideonella sakaiensis) to a target host species (B. subtilis), producing a gene sequence encoding a putative orthologous enzyme. We demonstrate that an OrthologTransformer-designed PETase gene enables B. subtilis to produce and secrete an active PET-degrading enzyme at levels surpassing those achieved by standard codon optimization. Our study focuses on DNA sequence generation and design; consequently, direct comparisons with protein-level generative models such as ProGen19 or ESM20 (which generate amino-acid sequences, often de novo or under high-level conditioning) are not appropriate. OrthologTransformer conditions on an input gene and performs cross-species rewriting at DNA resolution on prokaryotic genes, a fundamentally different objective from unconditional or tag-conditioned protein generation.

Results
A deep learning model for orthologous gene conversion
To enable gene sequence adaptation beyond synonymous codon changes, we developed OrthologTransformer, a sequence-to-sequence deep learning model that converts a coding DNA sequence from a source species into a predicted orthologous sequence in a target species. The model was trained on a large collection of known orthologous gene pairs from diverse bacteria. We implemented OrthologTransformer using the Transformer architecture, as illustrated in Fig. 1a, which is well-suited for modeling long sequences with complex dependencies. In essence, OrthologTransformer learns to edit the input gene, inserting or removing codon tokens and changing codons as needed to produce an orthologous sequence for the target species. Because the training ortholog pairs were naturally aligned by protein function, the model learned to make changes that preserve function (as true orthologs do) rather than introducing random, disruptive mutations.

During training, we additionally conditioned the model on the identity of the target species. In practice, a special token indicating the target species was prepended to each input sequence. This enabled a single model to handle translations into many possible target species (Fig. 1b). Through this mechanism, OrthologTransformer learns each species’ particular codon usage biases and typical amino acid adaptations. The final trained model takes an input gene (from a specified source species) and a desired target species, and it generates a new coding sequence that could function as that gene’s ortholog in the target species.

Benchmarking OrthologTransformer’s performance
We first evaluated how well the trained OrthologTransformer could generate known orthologous sequences across diverse species. In these tests, the model was given a gene from one species and tasked with predicting its ortholog in another species, and we measured the similarity between the AI-generated sequence and the actual native ortholog. Our benchmark datasets ranged from one-to-one ortholog pairs between 2 species (658 sequence pairs) to many-to-many ortholog relationships across 2138 species (over 4.9 million sequence pairs from the OMA database21), as shown in Supplementary Table 1. We found that the model’s performance improved dramatically with the breadth of training data available. This progressive expansion of training data yielded remarkable improvements in ortholog conversion accuracy, as demonstrated in Fig. 2, where the codon sequence identity between generated sequence and target sequence dramatically increased from 0.15 when using a very limited dataset (ortholog pairs from 2 species) to 0.40 when leveraging the full 2138-species dataset.

As part of our computational validation, we conducted a large-scale benchmark designed to match the practical upper limit of existing baseline methods. Specifically, we extended the performance comparison to the maximum possible scale, covering all 45 bacterial species included in both the OMA database (2138 species) and the set of species supported by the existing deep-learning method CodonTransformer22 (164 species), and evaluating 450 source-to-target species combinations (10 source species per target species). CodonTransformer, a recent deep learning model, uses a Transformer architecture to optimize synonymous codon choices in a context-aware manner, and is trained on large multispecies datasets22. The extensive comparison result is summarized in Fig. 3 using codon-level sequence identity to the target ortholog as the evaluation metric. In this large-scale benchmark, OrthologTransformer outperformed conventional frequency-based codon optimization and CodonTransformer across all evaluated source-to-target species combinations. OrthologTransformer consistently achieved significant improvements (typically with p values less than 1e-5) over synonym-focused approaches, indicating that capturing ortholog-like sequence patterns goes beyond what can be achieved by synonymous substitutions alone.

To provide concrete and biologically interpretable examples, we highlight nine representative species pairs spanning diverse genomic and physiological characteristics, where each species pair contains several thousand orthologous genes (Table 1). As shown in Table 1, the codon-level identity between the generated sequences and the true target sequences increases substantially: especially for conversions between species with markedly different genomic contexts. For example, when converting between B. subtilis (43.5% GC content) and I. sakaiensis (66.7% GC content), the identity to the target sequence doubled from 0.221 (original source sequence) to 0.424 (generated sequence). Likewise, between L. lactis and T. thermophilus, which differ substantially in optimal growth temperatures (30°C versus 65°C), the sequence identity increased more than three-fold—from 0.157 (original source sequence) to 0.467 (generated sequence).

Importantly, improvements are not limited to synonymous patterns. When examining sequence similarity at the amino acid level, OrthologTransformer preserves functional integrity while introducing biologically plausible amino-acid substitutions (that is, non-synonymous substitutions): Table 1 also reports BLOSUM-based amino-acid similarity, showing consistent improvements across multiple pairs. Many changes correspond to conservative substitutions (chemically similar amino acids), suggesting that the model preferentially introduces subtle, function-preserving protein-level edits when supported by ortholog supervision. In addition, OrthologTransformer captures species-specific sequence characteristics beyond raw identity: during conversion, GC content shifts toward the target genome, and CAI approaches the distribution of native target genes, as illustrated for representative pairs in Fig. 4: the distribution of CAI and GC content for OrthologTransformer-designed sequences (green) aligned much more closely with native target genes than did the original source sequences (blue).

Finally, these results translate into clear advantages over existing methods. Across the representative pairs, OrthologTransformer achieves approximately a 1.7-fold improvement in codon sequence identity compared to conventional codon optimization (Table 1). This suggests that synonymous-only optimization is insufficient to meet the contextual and evolutionary demands of host adaptation. In the subset of five species pairs where CodonTransformer was evaluated, OrthologTransformer still maintains a clear advantage, delivering an average improvement of ~1.8-fold. For example, when considering sequence conversion to E. coli, CodonTransformer achieves only slight improvements, whereas OrthologTransformer demonstrates more than a two-fold enhancement. These findings highlight the limitations of relying solely on synonymous substitutions and demonstrate that OrthologTransformer, by incorporating non-synonymous substitutions and indels, provides a more advanced and effective approach to gene adaptation than existing synonym-focused methods. Furthermore, we conducted statistically rigorous test (two-sided paired t-test) against conventional methods with three clearly defined evaluation metrics: codon-level sequence identity, CAI (codon adaptation index) proximity and GC alignment between source, target, and model generated sequences, and report statistical significance of improvements (Supplementary Table 2). Across those three metrics, OrthologTransformer consistently outperformed baseline codon optimization and, where comparable, CodonTransformer, with significant gains across pairs.

Designing a PETase enzyme for B. subtilis
Having validated the model’s general performance, we next applied OrthologTransformer to a specific biotechnologically relevant challenge: adapting a plastic-degrading enzyme (PETase) from Ideonella sakaiensis (a Gram-negative bacterium) to function in Bacillus subtilis (a Gram-positive host). I. sakaiensis PETase naturally breaks down PET plastic18,23,24, but I. sakaiensis grows slowly and is not an ideal organism for industrial use. In contrast, B. subtilis is a fast-growing, spore-forming bacterium amenable to large-scale fermentation25, making it an attractive host for PETase26–28—if the enzyme’s gene can be successfully expressed.

Using the model, we generated a set of candidate B. subtilis-adapted PETase sequences. To ensure we explored a broad design space, we employed additional computational optimization steps in conjunction with the model. In particular, we performed a multi-objective search using Monte Carlo Tree Search (MCTS) to refine the model’s outputs (Fig. 1c). This search aimed to optimize two key properties of the PETase coding sequence: (1) GC content around 36–37%, which is closer to B. subtilis genome composition, and (2) mRNA secondary structure stability, to maintain sufficient RNA structure for mRNA longevity29. We also evaluated a variant strategy where the OrthologTransformer model was fine-tuned specifically on orthologous gene data from I. sakaiensis and B. subtilis (to further specialize it for this species pair) before generating sequences. In total, we designed 12 distinct PETase gene variants by systematically varying (i) the breadth of training species (23, 54, or 2138), (ii) alignment processing (±), (iii) pair-specific fine-tuning (±), and (iv) MCTS multi-objective refinement (±). For interpretability, we group the variants by training breadth and data source: OrthoFinder-based (23/54 species; AI-S1, AI-M1–AI-M6) and OMA-based (2138 species; AI-L1–AI-L5). These twelve variants are summarized in Table 2, and the corresponding DNA sequences are provided in Supplementary Data 1. As illustrated in Fig. 5a, these modifications involved varying degrees of insertions, deletions, synonymous substitutions, and non-synonymous substitutions across the 12 AI-designed sequences. The extent of these changes varied dramatically among the different versions: from minimal modifications in AI-L1 (no changes) to extensive remodeling in AI-M3 (160 insertions, 139 deletions, 72 synonymous substitutions, and 30 non-synonymous substitutions). Despite such extensive sequence modifications, structural predictions showed that key functional domains remained intact.

Several AI-designed sequences were predicted to exhibit markedly improved properties in B. subtilis. As shown in Fig. 5b, c, one design in particular, AI-L2, achieved an optimal balance of characteristics: the highest predicted structural stability, a TM-score of 0.98 for the modeled 3D structure, a GC content of 37.0% (the wild-type I. sakaiensis PETase gene has substantially higher GC), and a favorable mRNA secondary structure (predicted folding free energy ΔG ≈ −281 kcal/mol for the full-length mRNA). This high TM-score supports the conclusion that the amino-acid differences in AI-L2 relative to wild-type PETase do not disrupt the overall fold. To provide complementary local assessments, we also examined backbone RMSD and per-residue pLDDT. As summarized in Fig. 5b, c, the redesigned sequences preserve the PETase fold (consistently high TM-scores), show small RMSD deviations, and display high pLDDT confidence across the fold. Notably, AI-L2 combines a TM-score ≈ 0.98 and the highest predicted structural stability with RMSD/pLDDT patterns consistent with near-identity to the wild-type structure, in line with its top functional performance.

Experimental validation of AI-designed PETase in B. subtilis
We synthesized a selection of the AI-designed PETase genes and tested their expression and activity in B. subtilis. Twelve constructs (AI-S1, AI-M1–AI-M6, AI-L1–AI-L5) were assembled, each encoding the designed PETase variant (or controls) under an inducible promoter on a B. subtilis shuttle plasmid. To facilitate secretion of the enzyme (since PET is an extracellular substrate), all constructs included an N-terminal signal peptide for Sec pathway secretion and a C-terminal 6×His-tag for detection (Supplementary Fig. 1). We included two control genes: WT, the wild-type I. sakaiensis PETase coding sequence, and CO, a purely codon-optimized sequence (identical amino acids, every codon replaced by the most preferred B. subtilis synonymous codon). All PETase constructs were transformed into B. subtilis, and expression was induced in shaking flask cultures.

mRNA transcription. All engineered B. subtilis strains successfully transcribed the full-length PETase mRNA. PCR amplification yielded the expected around 0.9 kb PETase transcript (including signal peptide and His-tag regions) in every AI-designed strain as well as in the WT and CO controls (Supplementary Fig. 2). Quantitative real-time PCR (qPCR) confirmed that all designs produced substantial PETase transcripts, although transcript levels varied by construct (Supplementary Fig. 3). Several variants (e.g., AI-S1, AI-M1, AI-M6, AI-L1, AI-L4) showed PETase mRNA levels comparable to or higher than the WT and CO controls, while a few (AI-M3, AI-M5, AI-L2) had somewhat lower transcript levels. Nonetheless, all AI-designed genes were robustly transcribed in B. subtilis.

Protein expression and secretion. Western blot analysis of culture supernatants (anti-His detection) showed that the PETase protein (~30 kDa) was present in many of the AI-designed strains (Supplementary Fig. 4), indicating successful expression and secretion. No PETase band was detected in the empty-vector negative control. Notably, variants AI-L1, AI-L2, and AI-L5 had especially strong PETase bands in the supernatant, comparable to or exceeding the CO control, indicating that those designs achieved particularly high secreted enzyme levels. The presence of PETase in the media confirms that the signal peptide functioned and the enzyme was exported out of the cells (a crucial feature for PET degradation, since PET is extracellularly located as a solid substrate).

PET degradation activity of AI-designed PETase. Finally, we tested the functional activity of each enzyme using a PET degradation assay. In this assay, B. subtilis cells expressing PETase were incubated with a film of additive-free PET plastic, and the breakdown products (terephthalic acid (TPA), mono (2-hydroxyethyl) terephthalate (MHET) and Bis(2-hydroxyethyl) Terephthalate (BHET)) were measured over time by HPLC method (Supplementary Fig. 5). The results confirmed that the AI-designed PETase is functional. B. subtilis cells harboring the AI-S1—AI-L5 genes exhibited clear PET degradation activity, comparable to cells with the codon-optimized PETase CO. PET hydrolysis was monitored by measuring its soluble breakdown products (TPA, MHET, and BHET) in the culture supernatant on days 1, 2, 3, and 7. Due to PETase’s known endo-activity, which predominantly generates MHET, the accumulation of MHET serves as the expected indicator of PETase activity in this assay system. MHET was detected as early as the 2nd day in some engineered strains (AI-M2, AI-M3, AI-L2, AI-L3, AI-L5, and WT), and by the 3rd day MHET had accumulated in all PETase-expressing cultures (Fig. 6a). Notably, the AI-L2 variant stood out by producing roughly three-fold more MHET than any other strain by day 3, reflecting significantly higher PET-degrading activity (p < 0.05; Fig. 6b). On the other hand, TPA and BHET remained below detectable (Supplementary Table 3), consistent with accumulation of MHET in the absence of MHETase, as expected for a PETase-only system.

Consistent with the MHET measurements, scanning electron microscopy (SEM) of the PET films after 7 days revealed the most extensive surface erosion in the AI-L2 treatment. The PET film incubated with the negative control (PET incubated with the empty-vector strain, denoted EV) remained mostly intact, showing a smooth surface with no obvious damage (Fig. 6d). In contrast, the film recovered from the AI-L2 culture displayed large cavities and holes approximately the size of bacterial cells (1–2 μm) penetrating into the film (large, cell-sized pits perforating the surface) (Fig. 6d), suggesting that colonies of B. subtilis on the film had extensively degraded the plastic directly beneath them. The other PETase strains also caused some surface roughening and etched patterns (Supplementary Fig. 6), whereas the negative control (EV) remained smooth. These qualitative observations align with the quantitative data, reinforcing that the AI-L2 strain achieves the most pronounced PET degradation among the strains examined.

Furthermore, we measured in vitro enzyme kinetics for WT, CO, and a redesigned PETase variant AI-L2. The enzyme kinetics (kcat and Km) were measured using BHET as substrate 30°C by quantifying the degradation of BHET (Bis (2-hydroxyethyl) Terephthalate) to the breakdown product MHET (mono (2-hydroxyethyl) terephthalate) over short time intervals (initial-rate assay). As shown in Supplementary Table 4, AI-L2 exhibited a pronounced kinetic advantage, driven by a substantial increase in kcat while maintaining a favorable Km. Consequently, its catalytic efficiency (kcat/Km ≈ 8.7 × 10³ M⁻¹s⁻¹) clearly exceeded that of WT and CO tested. Taken together, these observations indicate that PETase performance in vivo is governed by multiple interacting, non-linear factors—including mRNA abundance/stability/structure, translation efficiency/codon adaptation, and beneficial amino acid substitutions/codon usage adaptation, which in turn can affect mRNA expression, protein folding, and enzyme kinetics—rather than by any single determinant. The OrthologTransformer-designed sequence likely optimizes several of these axes simultaneously.

Another gene-host experiment in a Gram-negative host (Escherichia coli)
To assess portability across bacterial clades, we redesigned a gene for E. coli and compared it to a frequency-based codon-optimized control. In this experiment, OrthologTransformer was used to convert a gene from I. sakaiensis for expression in E. coli, using the same model configuration as the AI-L2 variant, which was pretrained on the OMA ortholog database with alignment preprocessing and MCTS-based optimization for GC content and mRNA stability. OrthologTransformer designs yielded higher expression and PET degradation activity in our initial test, supporting generality across distinct GC/regulatory backgrounds (see Supplementary Figs. 7–9 for expression data and Supplementary Data 2 for constructs data).

mRNA transcription: Full-length transcription of PETase mRNA was confirmed in all engineered E. coli strains (Supplementary Fig. 7a). qPCR analysis detected PETase transcripts in all designed strains (Supplementary Fig. 7b). Among them, the AI-E2 strain exhibited the highest transcription level, ~3.5-fold higher than the WT and about 10-fold higher than the CO control. Western blot analysis of culture supernatants (anti-His detection) showed that the PETase protein (~30 kDa) was present in WT, CO, and all the AI-designed strains (Supplementary Fig. 7c), indicating successful expression and secretion. Notably, AI-E2 showed strong band compared with CO.

PET hydrolysis was monitored by measuring MHET in the culture supernatant on days 1, 2, 3, and 7 (Supplementary Fig. 8a). MHET was detected as early as the 2nd day in AI-E2, and by the 3rd day the AI-E2 strain exhibited the strongest PET degradation (Supplementary Fig. 8b), while MHET had accumulated across all PETase-expressing cultures by the 7th day. After 7 days of incubation, the PET films treated with the AI-E2 strain exhibited distinct surface degradation traces, confirming that the cells actively degraded PET in E. coli (Supplementary Fig. 8c).

Notably, AI-E2 contains no amino-acid substitutions (Supplementary Fig. 9a); its divergence from WT and CO arises solely at the codon level (i.e., synonymous recoding). Consistent with this, qPCR and Western blot measurements showed that AI-E2 achieved the highest mRNA and protein levels, respectively, among the constructs tested. Thus, even without any protein-level changes, the OrthologTransformer-designed codon scheme delivers a measurable advantage over frequency-based codon optimization, suggesting that the model’s choices capture subtle host-specific constraints that are not fully addressed by conventional synonym-only approaches.

Finally, overall activity in B. subtilis exceeded that in E. coli; a plausible contributing factor is that B. subtilis forms abundant biofilm that promotes adhesion to PET film, whereas E. coli exhibits weaker biofilm formation under our conditions, reducing effective enzyme-substrate contact on the polymer surface (see SEM images in Supplementary Fig. 10).

dN/dS-informed analysis
To evaluate the evolutionary context of OrthologTransformer’s proposed amino acid substitutions, we performed a site-wise dN/dS analysis on alignments of homologous sequences using the FEL method30. For each site, the nonsynonymous (dN) and synonymous (dS) substitution rates were estimated; sites showing significant evidence of purifying selection (p < 0.01 and dN/dS < 1) were classified as “purifying”, whereas sites with dN/dS ≥ 1 or without significant purifying signal were classified as “non-purifying”.

Using this classification scheme, we compared each OrthologTransformer-generated sequence with its source ortholog sequence to determine the fraction of substitutions occurring at purifying versus non-purifying sites. We applied this analysis to nine representative species pairs. For each species pair, a Fisher’s exact test was used to assess whether substitutions were enriched at non-purifying sites compared to purifying sites (Table 3). In Table 3, the column “Significant Groups (p < 0.05)” indicates the number of orthologous sequence pairs that showed a statistically significant difference in substitution rates between site categories, out of the total pairs analyzed for that species comparison. Across all species pairs, non-purifying sites exhibited significantly higher substitution rates than purifying sites (Fisher’s exact test, p < 10⁻¹⁴). These findings indicate that OrthologTransformer preferentially introduces substitutions at evolutionarily variable sites, while changes at highly conserved positions are comparatively suppressed.

Methods
Datasets construction
Development of the OrthologTransformer model required large-scale, high-quality ortholog datasets to accurately capture diverse genetic variations, including synonymous and non-synonymous substitutions, as well as insertions and deletions (indels). A strict separation between training and test datasets was also essential to rigorously evaluate model generalization. Initially, we processed sequence data from 2 to 132 bacterial species retrieved from RefSeq using OrthoFinder42. To ensure high data quality, sequences were restricted to lengths under 2100 bp, and only ortholog pairs with length ratios between 0.97 and 1.03 were selected. For constructing a larger and more comprehensive dataset, we leveraged the OMA database, which provides highly reliable ortholog data identified through rigorous algorithms. Distinct orthologous groups were allocated separately to training and testing sets, enabling accurate evaluation of the model’s predictive versatility. We partitioned ortholog groups such that no ortholog pair in the test set shares a protein with any pair in the training set. This means if a particular gene (or any genes from its species pair) was used for training, none of its orthologs or that gene itself appear in testing. The finalized training dataset comprised 1,071,976 orthologous groups and 4,971,060 ortholog pairs derived from 2138 bacterial species.

To further enhance dataset quality, we applied gap-processing by removing alignment gaps following amino acid alignments, thus preserving conservation in functionally important regions. Both gap-processed and unprocessed datasets were prepared to maintain a balance between evolutionary information retention and alignment quality.

We selected the test set, we deliberately selected species pairs with diverse GC contents and optimal growth temperatures to examine how biological characteristics impact conversion accuracy. This multifaceted approach enables OrthologTransformer to achieve robust and highly accurate predictions even for genes it has not encountered previously.

OrthologTransformer architecture
We present OrthologTransformer, a sequence-to-sequence (seq2seq) model specifically designed for cross-species DNA sequence conversion (Fig. 1a). The model employs a Transformer architecture43 with a critical innovation: species-conditioning tokens that encode the direction of conversion.

Species-conditioning tokens: We prepend species identity tokens to the inputs to make the conversion direction explicit and to enable a single model to handle many source-target pairs with distinct sequence preferences. Concretely, the source token (s_src) is prepended to the encoder input, and the target token (s_tgt) is prepended to the decoder input, thus conditioning the model to learn directionality (e.g., converting a gene from species X to species Y) and is a deliberate design choice that allows a single model to handle multiple source-target species pairs.

Tokenization and embeddings: Input/output sequences are tokenized at codon resolution (three nucleotides per token). The vocabulary comprises 64 sense codons + start/stop + N species tokens (N is the number of distinct species used in training and test), and each token is mapped to a 512-dimensional learnable embedding. Embeddings are randomly initialized and optimized end-to-end during training. The context length is 700 tokens, which corresponds to ~2.1 kbp when using codon tokens, since each token is 3 bases.

Model configuration and optimization. The architecture comprises 20-layer encoder and 20-layer decoder stacks (total 114,095,940 parameters) with 16-head multi-head attention and position-wise feed-forward networks in each layer (Supplementary Method 1). We train with cross-entropy loss and the Adam optimizer using the following hyperparameters: embedding dimension 512, feed-forward dimension 1024, dropout 0.0, batch size 64, and initial learning rate 0.0001.

Allowance of non-synonymous substitutions and indels: Enabling amino-acid-changing edits and insertions/deletions (indels) is central to the biological task: many adaptive changes observed among natural orthologs require more than synonymous codon swaps, which synonym-only codon optimization cannot capture. Accordingly, OrthologTransformer’s output space and training objective explicitly accommodate length and amino-acid changes. During inference, the decoder generates codon tokens autoregressively and terminates upon an end-of-sequence condition, so insertions/deletions relative to the input emerge naturally.

Two-stage learning strategy
To obtain a single model that is both broadly transferable across bacteria and highly accurate for a specific engineering direction, we train OrthologTransformer in two stages (Fig. 1b).

Pretraining (many-to-many): We treat each ortholog pair as a directed sequence-conversion task from a source species (s_src) to a target species (s_tgt). Coding sequences are tokenized at codon resolution, and the conversion direction is explicitly encoded by prepending a source-species token to the encoder input and a target-species token to the decoder input; the vocabulary, therefore, consists of the standard codon tokens (plus start/stop) and one token per species. We optimize the model parameters by minimizing the token-level negative log-likelihood (cross-entropy) of the target coding sequence given the source sequence and the source/target species tokens across many-to-many source-target species combinations. Pretraining is performed on large multispecies ortholog corpora (ranging from small many-to-many sets up to the full OMA-based bacterial corpus). The largest pretraining set comprises 2138 bacterial species and 4,971,060 ortholog pairs (Supplementary Table 1).

Fine-tuning (one-to-one): We specialize the pretrained model for a particular source-to-target direction (e.g., I. sakaiensis to B. subtilis) by continuing training on ortholog pairs restricted to that specific species pair only. Fine-tuning uses the same codon tokenization, species-conditioning scheme, and cross-entropy objective, but concentrates learning on the pair-specific genomic and translational context of the target host.

This two-stage regime addresses the practical limitation that any single species pair typically provides far fewer ortholog examples than a multispecies corpus, while still enabling strong performance for the specific conversion direction used in downstream design.

PETase sequence design
To generate B. subtilis-compatible PETase sequences from I. sakaiensis, we strategically initialized the decoder with the target species token s_tgt concatenated with a codon-optimized signal sequence derived from the original PETase. This initialization guides the model toward generating sequences suitable for expression in the target organism. We implemented a comprehensive design strategy that explored multiple dimensions of sequence generation to create a diverse set of PETase variants (Table 2). Our approach systematically varied several key parameters: the source database (OrthoFinder or OMA) with corresponding variation in species diversity (ranging from 23 to 2138 species), whether training data underwent alignment preprocessing, whether the model underwent fine-tuning for the specific species pair, the application of MCTS for sequence optimization, and the choice of signal sequence (either codon-optimized or AI-designed). Regarding alignment preprocessing, we prepared both gap-retained and gap-removed datasets derived from amino-acid multiple-sequence alignments to balance conservation with edit visibility. This setup explicitly exposes the positions where natural orthologs introduce amino-acid substitutions and indels (gaps), enabling the model to learn those patterns. This multifaceted approach allowed us to generate twelve distinct sequences that collectively represent different combinations of biological knowledge integration, computational optimization, and AI-driven design, enabling us to evaluate which factors most significantly influence successful heterologous expression.

Monte Carlo tree search optimization
We developed a MCTS algorithm to optimize generated sequences for both GC content and RNA secondary structure stability. This addresses the competing requirements of synthetic accessibility (favoring low GC content around 0.36) and functional mRNA stability (requiring stable secondary structures). Specifically, while lowering GC content facilitates DNA synthesis and cloning, maximizing translational output requires mRNA secondary structures stabilized to an optimal degree, extending transcript half-life without hindering ribosomal access44,45. Our MCTS implementation builds a search tree where each node represents a partial DNA sequence, with child nodes adding additional codons based on Transformer-generated probability distributions. We select the top three most probable tokens for expansion at each step. Node evaluation employs a reward function that balances GC content optimization with secondary structure stability:

Reward = (10 / |100 * GC_generated - GC_target|) + |MFE| (1)

where GC_generated represents the sequence’s GC content, GC_target equals 0.36, and MFE denotes minimum free energy (MFE) of the RNA structure. The algorithm iteratively performs selection using UCB146, expansion based on model probabilities, simulation to terminal states, and backpropagation of rewards, enabling efficient exploration of the sequence space while maintaining biological constraints. The weighting scheme used in the MCTS reward function and the corresponding hyperparameter details are provided in Supplementary Method 2.

Evaluation metrics
To comprehensively evaluate the sequences generated by our approach, we employed multiple metrics spanning sequence identity, structural prediction, and biological compatibility. For sequence evaluation, we calculated codon sequence identity by performing pairwise alignment using the Needleman-Wunsch algorithm with match score 1 and gap penalty 0, followed by normalization against sequence length. Similarly, we computed amino acid sequence similarity by pairwise alignment based on the BLOSUM score, applying gap penalties of 10 for opening and 0.1 for extension, with scores normalized by sequence length. For structural assessment, we utilized AlphaFold3 to compute TM-scores comparing predicted structures of AI-designed PETase variants against the wild-type PETase structure, providing quantitative evaluation of structural conservation. We also analyzed GC content, a critical parameter affecting both synthetic accessibility and thermal stability, ensuring compatibility with the target organism’s genomic characteristics. Translation efficiency was evaluated using the CAI, calculated as:

CAI = exp (1/L * Σ ln w_ij) (2)

ere L represents gene length in codons, and w_ij denotes the relative adaptiveness of each codon, defined as:

w_ij = x_ij / x_imax (3)

Here, x_ij represents the frequency of synonymous codon j for amino acid i, while x_imax indicates the maximum frequency among synonymous codons for the same amino acid i. We computed CAI values using whole-genome codon usage frequencies as reference data for each species, with values ranging from 0 to 1, where values closer to 1 indicate optimal codon usage for the target organism.

Structure prediction tools
Protein structure prediction. We used AlphaFold3 (ver. 3.0.1)47 to predict the 3D structures of the redesigned PETase variants, as well as the wild-type and codon-optimized PETase, for structural comparisons. The default AlphaFold3 pipeline was used for each sequence, and the top-ranked predicted model was selected for analysis. We evaluated the predicted structures by calculating the TM-score and RMSD relative to the wild-type PETase structure. The per-residue pLDDT scores provided by AlphaFold3 were used to assess local model confidence. Additionally, the overall structural stability of each model was estimated using FoldX (ver. 4.0 Beta)48.

mRNA secondary structure prediction. We evaluated mRNA secondary structure for each designed sequence using the ViennaRNA package version 2.6.4 (RNAfold)49. We computed the MFE structure for the full mRNA sequence and recorded the MFE value as a measure of mRNA folding strength. This was done to compare whether our designs had potentially more favorable (less stable, thus more easily unwound) mRNA structures than the originals.

dN/dS analysis procedure
For each orthologous group, we selected close homologs (≥90% identity) using the source and target sequences as anchors and generated an amino-acid multiple sequence alignment, which was mapped to coding sequences to obtain a codon-sequence alignment. A phylogeny was inferred for FEL30, and HyPhy FEL was used to estimate site-wise nonsynonymous (dN) and synonymous (dS) substitution rates; sites showing significant evidence of purifying selection (p < 0.01 and dN/dS < 1) were classified as “purifying”, whereas sites with dN/dS ≥ 1 or without significant purifying signal were classified as “non-purifying”.

Experimental procedures
Strain construction. Synthesized genes and the expression vector (E. coli/B. subtilis shuttle vector pUC980-250, refer to Supplementary Method 3) were amplified by PCR (Oligonucleotides sequences of primers used in this study are provided in Supplementary Data 3). The resulting DNA fragments were mixed at a 1:1 molar ratio (vector:gene) and assembled using NEBuilder HiFi DNA Assembly Master Mix (New England Biolabs, USA). After incubation at 50°C for 30 min, the assembled DNA was transformed into E. coli HST08 (Takara, Japan), and plasmids were subsequently extracted. These plasmids were then transformed into B. subtilis RM12551 using a competence-based transformation protocol52, and transformants were selected on LB agar containing kanamycin (25 μg/mL). The strains and plasmids used or constructed in this study are listed in Supplementary Data 1.

RNA extraction and reverse transcription. Recombinant B. subtilis cells were harvested during the logarithmic growth phase (Supplementary Fig. 11) by centrifugation at 5800 × g for 3 min (refer to Supplementary Method 4 for details). Three biological replicates were prepared for each recombinant strain. The cells were resuspended in 100 μL of a 5 mg/mL lysozyme solution prepared in 1× TE buffer and incubated at 37°C for 10 min. RNA extraction was performed using the NucleoSpin RNA kit (Macherey-Nagel, Germany) according to the manufacturer’s instructions. Total RNA was quantified and assessed for quality using the RNA 6000 Nano Kit on an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). All samples exhibited an RNA integrity number (RIN) greater than 9.8. Reverse transcription of total RNA was performed using the ReverTra Ace qPCR RT Master Mix with gDNA Remover (Toyobo Co., Ltd., Japan) according to the manufacturer’s guidelines.

PCR for full-length transcript. PCR was then performed using a forward primer binding near the 5′ end of the PETase coding sequence (within the signal peptide region) and a reverse primer at the 3′ end (within the His-tag region) (refer to Supplementary Method 5 for primer design and Supplementary Data 3 for oligonucleotide sequences). Amplification was performed using PrimeSTAR Max DNA Polymerase Ver. 2 (Takara, Japan), with the following thermal cycling conditions: 20 cycles of 98°C for 10 s, 55°C for 10 seconds, and 68°C for 5 s. The expected amplicon was visualized on agarose gel. All PETase constructs yielded the correct band from cDNA, confirming presence of full-length mRNA.

Quantitative real-time PCR (qPCR). For quantification of PETase transcripts, we performed qPCR on a StepOne Plus Real-Time PCR System (Thermo Fisher Scientific, USA) with Power SYBR Green PCR Master Mix (Applied Biosystems, USA). Gene-specific primers were designed to amplify an internal region within the signal sequence of PETase gene, which is common to all recombinant strains (refer to Supplementary Method 5 for primer design and Supplementary Data 3 for oligonucleotide sequences). As an internal control, primers for the single-copy gyrA gene of B. subtilis (DNA gyrase subunit A) were used. The ΔCt for PETase was calculated relative to gyrA, and then expression levels were normalized to the WT (wild-type) PETase sample (set as 1.0) using the ΔΔCt method.

Western blotting (immunoblot)
To detect PETase protein, recombinant B. subtilis cells were grown in LB medium supplemented with 1% xylose and 25 μg/ml kanamycin at 30°C for 12 h, and samples of both the culture supernatant and cell pellet were collected. For the supernatant, 1 mL of culture was centrifuged (5800 × g, 5 min, 4°C) to remove cells, and the supernatant was mixed with 6× sample loading buffer. All samples were normalized by culture volume (since cell density was similar across strains at 24 h), and boiled (5 min, 95°C) before SDS-PAGE. Proteins were separated on 12.5% SDS-PAGE gels. For secreted proteins, we loaded 9.0 μL of culture supernatant samples directly. After electrophoresis, proteins were transferred to a PVDF membrane. The membrane was probed with an anti-polyhistidine monoclonal antibody (Proteintech, Cat# 66005-1-Ig, Lot# 1B7G5, 1:5000 dilution) against the 6× His-tag at PETase C-terminus, followed by an Alkaline Phosphate-conjugated secondary antibody (Proteintech, Cat# SA00002-2, 1:2000 dilution). Detection was performed with BCIP/NBT and recorded on an imager. A His-tagged protein ladder was used to verify molecular weight.

Enzyme kinetics (kcat and Km) in vitro. Measurement of the enzymatic activity on the degradation of BHET (Bis (2-hydroxyethyl) Terephthalate) to the breakdown product MHET (mono (2-hydroxyethyl) terephthalate) was performed by HPLC analysis of the decrease of BHET as substrate after 2-h incubation at 30°C. The specific activity was defined as the amount of enzyme needed for hydrolyzing 1 mmol of BHET per hour. Kinetic parameters Km and Vmax were determined using the corresponding substrate at concentrations from 0.1 mM to 3 mM under the same reaction conditions as described for substrate specificity determination. Kinetic parameters were derived from the Michaelis-Menten equation. All enzymes were purified using Capturem TM his-tagged purification maxiprep columns (Takara, Japan) and removing contaminating proteins using Nanosep (PALL, Japan), following the manufacturer’s instructions. The protein concentration was determined using the Bradford kit (ATTO, Japan).

PETase activity assessment
PET film degradation experiment. Additive-free PET film (thickness ~0.25 mm) was cut into ~20 cm² pieces (4 cm × 5 cm rectangles). Before use, the films were thoroughly washed with 1% (w/v) sodium dodecyl sulfate (SDS) for 30 min and rinsed with Milli-Q water, to remove any contaminants and surface treatments. Films were then soaked in 70% ethanol for 1 h for sterilization and air-dried in a laminar flow hood. Each PET film piece was placed in a 50 mL Erlenmeyer flask containing 20 mL of LB medium with kanamycin (25 μg/mL). B. subtilis transformant pre-cultures were grown overnight in LB medium containing kanamycin (25 μg/mL) at 30°C, and 200 μL of the overnight culture was used to inoculate the flask (initial OD ~0.1). 1% xylose was added at inoculation to induce PETase expression from the start. Cultures with PET film were incubated at 30°C, 180 r/min for up to 7 days. Flasks were covered but not sealed (to allow aeration). At predetermined time points (24 h, 48 h, 72 h, and 7 days), 0.5 mL of culture medium was sampled for HPLC analysis of PET degradation products, using the same protocol as for BHET samples. Because PET is solid and only its soluble hydrolysis products are in the medium, we focused on MHET and TPA detection. BHET was not expected unless PET hydrolysis stops at dimer stage (which PETase typically does not). Each assay was done in triplicate flasks for each strain. After 7 days, the PET film was recovered from the culture for imaging.

HPLC quantification of PET hydrolysis products. Culture supernatants from PET film flasks were processed as for BHET (Bis (2-hydroxyethyl) Terephthalate) samples (acetonitrile/formic acid quench and centrifugation). The HPLC analysis was identical. We confirmed linearity of MHET detection in the relevant range (0–3 mM) by running standard solutions. In control flasks (no PETase), no MHET or TPA peaks were observed at any time point. In PETase-expressing samples, a clear MHET peak emerged by day 2–3. MHET concentration was quantified for each sample. TPA was not observed (likely remaining bound in undetectable form or not produced due to lack of MHETase). The results were plotted as MHET (μM) released per flask over time. Statistical comparison of MHET at day 3 among strains was performed by one-way ANOVA with Dunnett’s multiple comparisons. AI-L2 produced significantly more MHET than all other strains (p < 0.0001 for AI-L2 vs each other PETase strain). These quantitative data corroborate the qualitative ranking of enzyme performance.

Scanning electron microscopy (SEM). After 7 days of incubation, PET film pieces were retrieved from the cultures with sterile tweezers. To remove attached cells and biofilm, each film was gently sonicated in 2% SDS solution for 5 min, then rinsed thoroughly with Milli-Q water (several dips). Films were then air-dried completely (Supplementary Method 6). The specimens were then plasma-coated with osmium for 5 s using an osmium coater (HPC-20, Vacuum Device) to prevent charging. Images were acquired using an Inspect S50 scanning electron microscope (FEI) operated with the xT microscope control software. The accelerating voltage was set at 12.5 kV with a working distance of 10 mm. No further image processing or manipulations were performed after acquisition. Images were acquired at 10,000× magnifications. We focused on 10,000× for detailed surface features. Multiple fields were imaged for each film. The control film showed a smooth surface with machining striations but no biological etching. Treated films showed various degrees of erosion: we observed pits, grooves, and in some cases, large holes that appear to correspond to where a bacterial cell was attached (as inferred by size and shape). These SEM results visually confirm PET hydrolysis and were consistent across triplicate samples.

Experiment in a Gram-negative host (Escherichia coli)
For the E. coli experiments, strain construction, RNA extraction and reverse transcription, Western blotting, and PETase activity assessment were performed following the same protocols as those used for B. subtilis. Details are provided in Supplementary Method 7.

Statistics and reproducibility
For the large-scale benchmark shown in Fig. 3, each plotted value represents the mean codon-level sequence identity across 10 source-to-target conversions per target species (n = 10), with error bars indicating the standard deviation. OrthologTransformer was compared against baseline methods using predefined evaluation metrics, including codon-level sequence identity, CAI proximity, and GC alignment. Statistical significance for pairwise method comparisons was assessed using a two-sided paired t-test where applicable, while Fisher’s exact test was used for categorical comparisons, such as enrichment of significant groups. Reported p values are provided in the corresponding Supplementary Tables.

Data were analyzed using ordinary one-way ANOVA to compare across the 15 experimental groups in Fig. 6b. The assumption of equal variance was confirmed by the Brown-Forsythe test (P = 0.3978). Following the ANOVA, post-hoc comparisons were carried out using Dunnett’s multiple comparisons test to determine significant differences between specific groups. Statistical analyses were performed using GraphPad Prism, v10.1.0.

No statistical method was used to predetermine sample size.

All computational analyses were performed using custom Python scripts and are reproducible given the same input data and parameters. Source code and datasets are available as described in the Code Availability and Data Availability statements.

All unique strains and plasmids generated in this study are available from the corresponding author upon reasonable request.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "1.The experimental validation, while elegantly executed, focuses on a single enzyme (PETase) and a single host species (B. subtilis). While this case is illustrative, the generality of the model’s predictions remains underexplored. Either expand experimental validation to additional gene-host pairs (ideally with contrasting GC content or expression challenges) or better contextualize this limitation by discussing in the Discussion section why PETase was chosen and whether this design paradigm can generalize to less robust or more structurally sensitive proteins.\n2.The model introduces non-synonymous mutations, which could potentially affect enzyme kinetics, stability, or interactions. Although the structural modeling (e.g., TM-score) is used as a proxy for function preservation, it is insufficient to rule out functional divergence. Include deeper in silico evaluation of biochemical properties (e.g., stability predictions, hydrophobicity shift, domain conservation) of the introduced mutations. Experimental validation of kinetic properties (e.g., kcat, Km) of PETase variants would be ideal.\n3.While the model leverages ortholog data, the paper lacks discussion of whether the model is indeed learning evolutionary constraints or merely performing large-scale memorization. Are the mutations enriched in evolutionarily variable regions? Do they co-locate with natural ortholog divergence?\n4.The manuscript compares OrthologTransformer to codon optimization and CodonTransformer, but other tools (e.g., ProGen, ESM-based design frameworks) may also support non-synonymous generation or function-guided editing. Clarify how OrthologTransformer differs from such models, particularly in its ability to generate ortholog-like rather than just function-retaining sequences. Are there trade-offs in computational cost (e.g., training time, hardware requirements) between methods?\n5.While the core framework of OrthologTransformer is based on the Transformer architecture, it remains unclear to what extent the model has been specifically designed to address the challenges of cross-species gene redesign beyond simply applying a generic seq2seq framework to biological sequences. For instance, the use of species-conditioning tokens, the two-stage training (pretraining + species-specific finetuning), and the sequence-level editing capabilities (e.g., non-synonymous substitutions, indels) are potentially task-driven innovations. However, the manuscript does not clearly articulate the design motivations behind these choices or how they directly address the biological and functional challenges of adapting genes across evolutionarily distant hosts.\n6.The AI1–AI12 naming is acceptable, but consider grouping variants into meaningful categories based on mutation type or optimization strategy, which would aid interpretation.\n7.While the PET degradation assay is compelling, key functional data are missing. Only MHET (an intermediate) was quantified, not TPA (the final product). Does OrthologTransformer-designed PETase fully degrade PET, or does it stall at MHET? SEM images qualitatively show degradation, but quantitative analysis would strengthen the conclusions. Enzyme kinetics are absent. Are the activity improvements due to higher expression, better folding, or enhanced catalytic efficiency?\n8.The manuscript states that datasets and code are available on GitHub, but critical details are unclear. Is the test set completely independent of training data? Potential data leakage could inflate performance metrics.\n9.Can OrthologTransformer handle eukaryotic genes (e.g., intron-containing sequences)? For longer genes (>2 kb), does performance degrade due to Transformer context limits? Could AlphaFold-predicted structures guide non-synonymous mutation selection?\n10.How dose the model balances sequence adaptation with functional conservation? What features (e.g., conserved domains, structural stability) are prioritized during sequence generation? Are non-synonymous mutations guided by known evolutionary principles (e.g., BLOSUM62 scores)?", "reviewer_2": "In this study, OrthologTransformer, an AI-powered gene design tool, was developed to optimize heterologous gene expression across species. Its superior performance over conventional codon optimization was presented through adaptation of PETase in Bacillus subtilis. While the study addresses an important challenge in synthetic biology, several critical issues should be resolved.\n1.Limited model innovation:\nThe innovative contribution contribution of this algorithmic model appears limited, as the authors have essentially adopted the standard Transformer architecture without substantial modifications tailored to the specific biological problem. The claimed 'critical innovation' of simply adding special tokens to the encoder and decoder inputs does not represent a significant advancement from the original Transformer's application in machine translation.\n2.Insufficient methodological comparisons:\nThe manuscript would benefit from quantitative comparisons between their algorithm and traditional approaches. To properly demonstrate the model's advancement, they should conduct statistically rigorous benchmarking against conventional methods with clearly defined evaluation metrics.\n3.The manuscript should include a detailed description of the tokenization strategy in the main text. In addition, justification for using AlphaFold-predicted TM-scores for structural comparison was recommended rather than establishing structural alignment algorithms (e.g., RMSD calculations from experimentally determined structures)\nSpecific comments on the PETase case study:\n1.In the abstract section, it mentions that the OrthologTransformer-redesigned PETase showed 'robust activity surpassing codon-optimized controls'. However, this claim would benefit from quantitative validation (e.g., fold-improvement in enzymatic activity or expression levels) that would clearly demonstrate the magnitude of enhancement over traditional codon optimization methods.\n2.In this study, the PETase case study only demonstrates adaptation in B. subtilis (a Gram+ bacterium). Have the authors validated the model's performance in Gram-negative hosts or eukaryotic systems (e.g., yeast), given their distinct codon usage and regulatory elements?\n3.In the “PET degradation activity of AI-designed PETase“ section (lines 235-250), the authors described that “the breakdown products (terephthalic acid (TPA), mono(2-hydroxyethyl) terephthalate (MHET) and Bis(2-hydroxyethyl) Terephthalate (BHET)) were measured over time by HPLC method (Supplementary Figure S5)”. However, in Supplementary Figure S5, we only found the PET degrading pathway. Similarly, in Figure 5, the “Quantification of PET degradation products by HPLC” only showed the quantification of MHET hydrolyzed by PETase from PET. Since the HPLC methodology was used to measure TPA/BHET, the corresponding chromatograms and standard curves should be provided to validate the assay sensitivity. The TPA/ BHET data analyzed by HPLC should be provided in a supplementary panel, and discussed why MHET is the most relevant readout, which could further confirm “the AI-designed PETase is functional.”", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "Disclaimer: Please note that this review was done as a co-review, involving an early career researcher from my group in the context of Nature Communications’ commitment to facilitating peer review training and appropriately recognizing contributions.\n---\nThis paper presents a substantial contribution to synthetic gene design. The proposed tool, OrthologTransformer, aims to fill a gap left by codon optimization methods that restrict themselves to synonymous substitutions. The authors define their work as a new category of tool that we might call an ortholog designer. We agree that this is an interesting and worthwhile direction: rather than preserving the amino acid sequence at all costs, the presented model seeks to mimic the kinds of sequence changes observed in orthologs such as non-synonymous mutations and INDELs, while maintaining protein function, and engineering a gene/protein that fits the genetic background of the anticipated target bacteria species.\nIn addition to presenting the method, the paper includes a strong real-world application example, making it well-rounded by combining model development with meaningful experimental validation.\nBelow are some questions and suggestions that may help improve the manuscript further and should be considered before publication.\n# Major\n1) The model uses a large dataset of orthologs from species with varying evolutionary distances, but evolutionary distance is not explicitly included in the modeling. Could this be a limitation? Including information about their distance directly might help the model learn more accurately how orthologs evolve across different species.\n2) Some of the most successful redesigned variants, particularly AI8, AI9, and AI12, are not fine-tuned and contain few or no non-synonymous mutations. In particular, the AI9 variant demonstrated the highest PET-degrading efficiency among tested strains, producing threefold more MHET, and causing the most pronounced PET surface erosion - while only comprising a single non-syn mutation. The section discussing the most successful variants could be expanded to include a discussion of these observations and their implications. Particularly, because you emphasize the ability of OrthologTransformer not only to model synonymous changes. However, the current results give the impression that the most powerful redesigned proteins, at least for your PETase example, are indeed those with few/no non-synonymous mutations. While we agree that it can be beneficial to also model changes beyond synonymous changes, this might question your approach, focusing also on non-syn/INDELs, and can be discussed further.\n3) It is not very clear how the 12 constructs (AI1–AI12) were selected. The manuscript could clarify which parameters and combinations were considered relevant for testing. Was there a larger initial pool of candidate designs? If so, why were certain combinations excluded from experimental validation?\n4) Figure 1A would benefit from a more straightforward graphical layout or explanation, as it is not immediately clear what the input and output of the model are in this figure. Along the line, the figure caption could be more elaborated to help understand the architecture.\n5) In the same context, how the embedding was done is also not entirely clear to us. You mention the level of granularity of the tokens (codons), and the final embedding dimensionality (512) in lines 341-343, but not how the embedding was initialized or derived. Did you do a one-hot intermediate step or direct learnable embedding lookups? How were the embeddings initialized? Maybe you skipped describing that part because the one-hot encoding is implicitly done while embeddings are looked up; however, it's better to describe it explicitly.\n6) How were the protein structures predicted in 4a? This is missing from the methods. Also, you mention mRNA folding kinetics, and according to the code, the ViennaRNA package is used. This should also be described/mentioned in the Methods.\n# Minor\n7) This comment is more out of curiosity and not necessarily something you have to mention in the manuscript: your approach seems to have some methodological similarity with the detection of sites under negative and positive selection (using dN/dS). We wonder if the codon sites changed by the model via non-synonymous mutations would also be detected as under positive selection pressure using approaches such as CODEML or the HyPhy suite (FEL, MEME, …) on the MSAs of orthologous genes. We believe that this could have interesting applications in designing proteins under varying selective pressure constraints.", "reviewer_6": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Dear Reviewers,\nWe thank you for the thorough evaluation of our manuscript “Cross-species gene redesign leveraging\northolog information and generative modeling.” We appreciate all the constructive comments from Reviewer\n#1 to Reviewer #6. We have carefully revised the manuscript to address every point raised. Below we provide\na point-by-point response to each comment. All corresponding changes in the manuscript are highlighted in\nred in the revised Word document as requested.\nIn summary, the main revisions include the following:\n1. Additional host validation (Gram‑negative):\nWe added a new case study in Escherichia coli to complement the B. subtilis experiments, broadening\nthe biological context and demonstrating generality across distinct GC/regulatory backgrounds. The\nResults and Supplementary materials document the constructs, expression data, and PET-degradation\nassays for E. coli.\n2. Rigorous computational benchmarking:\nThe rigorous computational benchmarking was performed across multiple representative species pairs\nspanning diverse genomic or physiological characteristics, using large gene sets per pair (around\nseveral thousands of orthologous genes per pair). We conducted statistically test against conventional\nmethods with three clearly defined evaluation metrics: codon-level sequence identity, CAI proximity\nand GC alignment between source, target, and model generated sequences, and report statistical\nsignificance of improvements.\n3. Clarified novelty and scope:\nWe emphasize that our contribution is a new problem formulation at the DNA level – orthologous\ngene conversion framed as a sequence‑to‑sequence “translation” task (distinct from protein‑level\ngeneration). OrthologTransformer is trained in a supervised seq2seq manner on naturally occurring\northolog pairs, enabling it to learn the evolutionary balance between host adaptation and functional\nconservation. Consequently, the model does not decide this balance via an explicit objective or hand\ntuned weights; rather, it reproduces the patterns present in orthologs.\nReviewer #1\nComment 1 (Single enzyme and host (PETase in B. subtilis) – generality of the model): Reviewer #1 was\nconcerned that our experimental validation focuses on one enzyme-host pair, leaving the broader applicability\nof the model underexplored. They suggested either adding another gene-host example (ideally in a different\norganism such as a eukaryote) or, at minimum, discussing in the manuscript why PETase was chosen and\nwhether the design paradigm can generalize to other proteins.\nResponse: We are actively extending our experiments to additional gene-host pairs. As suggested, we have\ninitiated a new case study in a different organism (a Gram-negative bacterium, E. coli). In this supplementary\nexperiment, OrthologTransformer was used to convert a gene from I. sakaiensis for expression in E. coli, and\nwe observed improved expression and PET degradation activity compared to a codon-optimized control (data\nnow included in new Supplementary Figures S8–S10 and Supplementary Table S7 in the revision). These\nresults provide supporting evidence that our model’s benefits generalize to a distinct host with different GC\ncontent and regulatory context. We have also added a brief discussion on potential application to eukaryotic\nsystems (e.g. yeast): because our model operates on coding DNA sequences, it can in principle handle\neukaryotic genes after removal of introns, but we note that eukaryotes introduce additional complexities (such\nas splicing and chromatin context) that would require further investigation. Moreover, extending activity\nassays beyond PETase would require building enzyme-specific assay systems de novo, which is non-trivial\nand time-intensive; accordingly, we designate systematic wet-lab validation for additional enzymes/genes as\nfuture work.\nIn addition, we have expanded our Discussion to clarify why PETase was selected and to address the potential\ngeneralizability of our approach. PETase (a robust enzyme with a clear activity assay) was chosen as an initial\nproof-of-concept because it allowed us to isolate the effects of our gene redesign in a controlled setting. We\nnow explicitly discuss how OrthologTransformer is a general platform not limited to PETase: because the\nmodel was trained on a broad, multispecies ortholog dataset, it inherently learns patterns of gene adaptation\nthat are applicable across diverse genes and hosts. This is a key novelty of our approach – by leveraging many\northologous gene pairs during training, the model captures general principles of cross-species adaptation,\nrather than being tailored to a single gene. We have added text to emphasize that the concept of using a\nTransformer-based model to perform “orthologous gene conversion” is broadly applicable, in principle\nenabling redesign of any gene for any target species, which has not been explored by prior methods.\nComment 2 (Non-synonymous mutations and functional effects – deeper analysis requested): Reviewer\n#1 pointed out that our model introduces non-synonymous mutations, which could affect enzyme kinetics,\nstability, or interactions. They noted that while we used structural modeling (AlphaFold2 TM-scores) as a\nproxy for function preservation, this alone may be insufficient to rule out functional divergence. The reviewer\nrequested a deeper in silico evaluation of biochemical properties (e.g., stability predictions, hydrophobicity\nshifts, domain conservation) of the introduced mutations, and ideally experimental validation of kinetic\nparameters (kcat, Km) for the PETase variants.\nintroduced non-synonymous changes do not adversely affect protein function. In the revised manuscript, we\npresent a more comprehensive in silico characterization of the mutations introduced by OrthologTransformer\n(in new Figure 4b and 4c). Specifically, we evaluated structure stability for each variant with amino acid\nsubstitution using protein stability prediction tools. As shown in Figure 4b and 4c, variants differ in structural\nconservation (TM-score, backbone RMSD, per-residue pLDDT), predicted structural stability, RNA\nsecondary-structure free energy (MFE), and GC content. The top performer AI9 (renamed to AI-L2 in the\nrevision) achieves an optimal balance of characteristics: the highest predicted structural stability, a TM\nscore of 0.98 for the modelled 3D structure, a GC content of 37.0%, and a favorable mRNA secondary\nstructure (predicted folding free energy ΔG ≈ −281 kcal/mol for the full-length mRNA).\nWith respect to experimental validation of enzyme kinetics (kcat and Km), we agree that this would be ideal.\nWe measured in vitro kinetics (kcat and Km) by quantifying MHET formation over short intervals (initial-rate\nassay) for three representative samples, AI9 (renamed to AI-L2 in the revision), wild-type, and\ncodon-optimized. As summarized in Supplementary Table S6, AI9 (renamed to AI-L2 in the revision)\nexhibited a pronounced kinetic advantage, driven by a substantial increase in kcat while maintaining a\nfavorable Km. Consequently, its catalytic efficiency (kcat/Km) clearly exceeded that of the wild-type enzyme\nand the codon‑optimized control tested.\nWe additionally examined mutation positions with respect to conserved domains and relate these observations\nto the dN/dS analysis presented below, as Reviewer #5 suggested. We have included this analysis in the\nResults (and provided details in new Table 4).\nThese jointly optimized multiple properties contribute in a non-linear manner for enhanced activity.\nComment 3 (Evolutionary constraints vs. memorization – is the model learning biology?): Reviewer #1\nasked whether the model is genuinely learning evolutionary constraints from ortholog data or merely\nperforming large-scale memorization. They inquired if the mutations the model introduces are enriched in\nevolutionarily variable (less conserved) regions, and whether these changes co-locate with divergence\nobserved in natural orthologs.\nResponse: We have conducted additional analysis to confirm that OrthologTransformer is indeed capturing\nevolutionary relevant patterns rather than naively memorizing sequences. In the revision, we report that the\nmodel’s substitution choices strongly correlate with known evolutionary variability. Specifically, we mapped\nthe non-synonymous mutations proposed by OrthologTransformer onto a multiple-sequence alignment of\nhomologous sequences (including the original source and target orthologs and related orthologs from other\nspecies) and did dN/dS-informed analysis (shown in new Table 4). We found that the majority of positions\nwhere our model introduced an amino acid change correspond to sites that are not highly conserved across the\nfamily. In other words, OrthologTransformer tends to alter residues that natural evolution also frequently\nvaries, while leaving highly conserved positions unchanged. This suggests the model has learned which\nregions of the protein can tolerate change (likely reflecting functional flexibility) versus which regions must\nbe preserved for function. We have added these findings to the Results section.\nESM, etc.) – clarify OrthologTransformer’s unique contributions:): Reviewer #1 noted that our manuscript\ncompares OrthologTransformer to traditional codon optimization and to a baseline called\n“CodonTransformer”, but pointed out that other tools (e.g., ProGen, ESM-based design frameworks) can also\ngenerate non-synonymous changes or do function-guided protein design. The reviewer requested clarification\non how OrthologTransformer differs from such models, particularly in its ability to generate ortholog-like\nsequences rather than just function-retaining sequences. They also wondered about potential trade-offs in\ncomputational cost between our approach and others.\nResponse: Our study focuses on the generation and design of DNA sequences, not protein sequences.\nTherefore, it is not appropriate to directly compare our method with protein-level generative models like\nProGen or ESM, which operate on amino acid sequences.\nCompared to such protein-level generative models, our approach tackles a different problem formulation.\nProGen is a large language model trained on millions of protein sequences to generate entirely new protein\nsequences with specified functions or properties. In contrast, OrthologTransformer does not generate proteins\nfrom scratch; it works as a sequence-to-sequence (seq2seq) model that converts an existing gene sequence into\na version optimized for a target species. In other words, given a coding DNA sequence from one organism, it\nleverages the context of that input sequence and knowledge of orthologous variants to “translate” it into a\nDNA sequence better suited for another organism. This seq2seq formulation (analogous to translating a\nsentence from one language to another) is fundamentally different from ProGen’s generation approach. To our\nknowledge, OrthologTransformer is the first work to apply a Transformer-based seq2seq model for direct\northologous gene conversion at the DNA level, whereas prior protein large language models like ProGen\noperate at the amino acid level and do not perform cross-species rewriting of a given input sequence.\nWe have added a brief discussion of computational considerations. Training OrthologTransformer required a\nsizable dataset and GPU resources, comparable to training other Transformer models. However, we note that\nour model’s size is moderate (embedding dimension 512, 20 layers) and much smaller than large protein\ngenerative models such as ProGen and ESM, making training feasible on available hardware. In inference,\nOrthologTransformer can generate a redesigned gene in a matter of seconds on a single GPU, which is a\none-time cost per gene. By contrast, ProGen’s largest models (with billions of parameters) demand more\nmemory and are not tailored to specific input sequences (one would need to filter or fine-tune their outputs for\neach new gene), and traditional codon optimization tools, while extremely fast, lack the sophistication to\nexplore non-synonymous changes. We have included these points to reassure the reviewer that our method’s\ncomputational requirements are reasonable and to further underline the novel capability it delivers relative to\nprior art.\nComment 5 (Model design motivations – tailoring the Transformer to cross-species gene redesign):\nReviewer #1 observed that our framework is based on a standard Transformer, and asked to what extent the\nmodel has been specifically designed for the challenges of cross-species gene adaptation. They listed\n4\nexamples (species-conditioning tokens, two-stage training, allowing non-synonymous substitutions and indels)\nas potentially important innovations, but felt the manuscript did not clearly\narticulate the rationale behind these choices or how they address biological challenges.\nResponse: We appreciate this detailed prompt to clarify our design rationale. In the revised manuscript, we\nhave added subsections in the Methods to explicitly describe our model design choices and their motivations.\n• Species-conditioning tokens: We explain that incorporating species identity tokens was crucial so that the\nmodel could learn the “direction” of conversion (e.g., gene from species A to species B). This addresses\nthe challenge that different target species have different codon biases and sequence preferences – by\ngiving the model an indicator of the source and target context, we enable it to modulate its predictions\naccordingly. Without such tokens, a “single” model would struggle to generalize across diverse\nevolutionary distances. In the revised manuscript, we now emphasize that this choice was made to imbue\nthe model with flexibility and specificity in generating host-tailored genes.\n• Two-stage training (pretraining + fine-tuning): We now clearly state that we adopted a two-stage training\nregime to address data scarcity for any given species pair. During pretraining, the model sees a broad\nrange of orthologs (many-to-many) which teaches it general principles of gene evolution (e.g., typical\nmutations that occur between distant bacteria, basic codon usage patterns, etc.). During fine-tuning, we\nspecialize the model to a particular source-target species pair (one-to-one direction) so it can capture the\nnuances of that pair’s genomic context. Thus, we justify that this design was chosen to ensure both\ngenerality and precision, directly tackling the challenge of varying evolutionary distances. This reasoning\nis now explicitly described in the revised Methods.\n• Allowance of non-synonymous substitutions and indels: We underscore that enabling amino‑acid\nchanging edits and indels is central to the biological problem: many adaptive changes observed among\nnatural orthologs require small non‑synonymous substitutions and occasional insertions/deletions, which\nconventional codon‑optimization (synonym‑only) frameworks cannot capture. Accordingly, we designed\nOrthologTransformer’s output space and training objective to explicitly accommodate length and\namino‑acid changes. Concretely, the model employs an autoregressive decoder that generates codon\ntokens to an end‑of‑sequence condition, so insertions/deletions relative to the input emerge naturally, and\nwe train by maximizing the likelihood of true orthologs. The resulting variants in Figure 4a illustrate that\nthe model indeed proposes edits spanning indels and non‑synonymous substitutions.\nIn addition, we intentionally leveraged alignments during data preparation to expose the model to explicit\npatterns of non‑synonymous changes and indels. The Methods clarify that we prepared both\ngap‑processed and unprocessed datasets after amino‑acid alignment to balance conservation with edit\nvisibility; this was done precisely so the model could learn where natural orthologs introduce amino‑acid\nchanges and gaps. In the conversion from I. sakaiensis to B. subtilis, alignment‑based training sometimes\nimproved conversion accuracy relative to alignment‑free training (see the gray vs. black bars in Figure 2),\nconsistent with the idea that alignments provide a clearer supervisory signal for position‑specific\nsubstitutions/indels. We now make this rationale explicit in the revision.\nReviewer #1 noted that while naming the designed variants as AI1–AI12 is acceptable, it might aid\ninterpretation to group these variants into meaningful categories (e.g., by the type of mutations they contain\nor the optimization strategy used), rather than treating them as an arbitrary list.\nResponse: We agree that clearer categorization of the designed variants would aid interpretation. In the\nrevision, we have adjusted how we present AI1–AI12 (renamed to AI-S1–AI-L5 in the revision) to make their\ndifferences more transparent and aligned this with the underlying design parameters.\n• Moved & expanded table: The former Supplementary Table S2 that concisely summarized each variant’s\ndesign “recipe” (database/source, number of species, alignment processing, fine-tuning, MCTS, signal\nsequence) has been moved into the main text as new Table 3, and expanded with a brief, step-by-step\nnarrative of how each model/design was generated under systematically varied conditions – training\nsource (Orthofinder vs OMA), breadth of training species (23, 54, or 2,138), alignment processing (±),\npair specific fine tuning (±), MCTS (Monte Carlo Tree Search Optimization) to balance GC and RNA\nMFE.\n• Subgrouping by training breadth (major distinction): Following the reviewer’s suggestion to group by\nstrategy and given that the number of species in the training data most strongly differentiates model\nbehavior, we now grouped and renamed variants by training breadth/source in the Results and in a revised\nfigure:\n Orthofinder-based (23/54 species): AI‑S1 (23 sp.), AI‑M1–AI‑M6 (54 sp.)\n OMA-based (2,138 species): AI‑L1–AI‑L5 (2,138 sp.)\nThis grouping is now stated in the text and visually demarcated, so readers can immediately see\nhow training diversity relates to downstream properties and performance.\nThese changes address the request by replacing an essentially alphabetical list with meaningful,\nstrategy-aware subgroups, while keeping per-variant edit/metric summaries prominent for readers who wish to\ncompare designs at the mutation level.\nComment 7 (Missing functional data in PET degradation assay – quantification of final product (TPA)\nand kinetics): Reviewer #1 expressed concern that key functional data were missing from the PETase\nexperiment. They observed that only MHET (an intermediate) was quantified, not TPA (the final monomeric\nproduct of PET hydrolysis), and they requested more quantitative analysis to strengthen conclusions. They\nalso mentioned the absence of enzyme kinetics data, asking whether the observed activity improvements are\ndue to higher expression levels, better folding, or enhanced catalytic efficiency of the enzyme.\nResponse: We have addressed these points thoroughly by augmenting our experimental analysis and\ndiscussion, as Reviewer #2 also pointed out:\n• We have now included data on TPA (and BHET) detection. Although the redesigned PETase did not\nproduce a significant amount of TPA (due to the reasons explained: lacking MHETase to complete PET\ndepolymerization), we wanted to ensure the assay’s completeness. The supplementary table\n(Supplementary Table S5) now contains the chromatographic evidence that TPA was monitored, and we\n6\nexplicitly state that TPA was not detected above background in any samples for either wild-type,\ncodon-optimized, or AI-designed PETase. This clarifies that the improved PET degradation we report is in\nterms of increased MHET release and PET weight loss, not the production of end-state monomers.\n• We added analyses of enzyme kinetics and factors underlying the performance gains. We measured in\nvitro kinetics (kcat and Km) by quantifying MHET formation over short intervals (initial-rate assay) for\nthree representative samples, AI9 (renamed to AI-L2 in the revision), wild-type, and codon-optimized. As\nsummarized in Supplementary Table S6, AI9 (renamed to AI-L2 in the revision) exhibited a pronounced\nkinetic advantage, driven by a substantial increase in kcat while maintaining a favorable Km.\nConsequently, its catalytic efficiency (kcat/Km) clearly exceeded that of the wild-type enzyme and the\ncodon‑optimized control tested. Beyond enzyme expression and kinetics factors, multiple DNA-level\nfeatures contribute in a non-linear manner. As shown in Figure 4b and 4c, variants differ in structural\nconservation (TM-score, backbone RMSD, per-residue pLDDT), predicted structural stability, RNA\nsecondary-structure free energy (MFE), and GC content. The top performer AI9 (renamed to AI-L2 in the\nrevision) achieves an optimal balance of characteristics: the highest predicted structural stability, a TM\nscore of 0.98 for the modelled 3D structure, a GC content of 37.0%, and a favorable mRNA secondary\nstructure (predicted folding free energy ΔG ≈ −281 kcal/mol for the full-length mRNA). These jointly\noptimized properties provide a mechanistic rationale for enhanced activity.\nComment 8 (Dataset and code clarity – test set independence and potential data leakage): Reviewer #1\nnoted that our manuscript states datasets and code are available, but asked for clarification on whether the\ntest set is completely independent of the training data. They were concerned that any overlap (data leakage)\ncould inflate performance metrics.\nResponse: We apologize for not making this point clearer initially. We have updated the Methods section to\nexplicitly describe how we constructed training and test splits to avoid any data leakage. We obtained\northologous gene pairs from the OMA database. We then partitioned ortholog groups such that no ortholog\npair in the test set shares a protein with any pair in the training set. This means if a particular gene (or any\ngenes from its species pair) was used for training, none of its orthologs or that gene itself appear in testing.\nComment 9 (Eukaryotic genes and long sequences – model limitations and potential use of structure\nguidance): Reviewer #1 asked whether OrthologTransformer can handle eukaryotic genes (for example,\ngenes with introns) and how it performs with longer genes (>2 kb) given Transformer context length limits.\nThey also wondered if incorporating AlphaFold-predicted structures could guide non-synonymous mutation\nselection.\nResponse: We address these points as follows in the revised manuscript:\n• Eukaryotic genes (intron-containing): We clarify that our current implementation of OrthologTransformer\nis designed for continuous coding sequences (CDS). To apply it to a eukaryotic gene with introns, one\nwould use the spliced mRNA/cDNA sequence (exons only) as input and output. The model itself has no\nfundamental limitation in handling eukaryotic coding sequences. However, we have added a note in the\noptimization with separate considerations for promoter, UTRs, and possibly alternative splicing variants.\nWe indicate that this is an area for future work, but that in principle the model could be extended to\neukaryotes by including intron-less ortholog pairs in training and perhaps adding tokens for other relevant\nfeatures if needed. This acknowledges the limitation while expressing that it’s not an insurmountable one,\njust outside our current scope.\n• Long genes (>2 kb) and Transformer context: In our Methods, we now specify the maximum sequence\nlength the model was trained on. We used positional encodings up to a certain length (in our case, 700\ntokens, which corresponds to ~2.1 kbp when using codon tokens, since each token is 3 bases). Thus,\ntypical genes (which are often 1 kbp or less in bacteria) are well within the context window. We do note\nthat extremely long genes (much beyond 3-4 kbp) could approach the limit of the model’s context. In\npractice, we did not encounter significant degradation in performance for the longest genes in our test. We\nhave added this information to reassure that the model can handle genes in the few-kilobase range. So, we\nacknowledge the context limit but indicate it did not affect our results and is a known consideration for\nfuture scaling.\n• AlphaFold-guided mutation selection: We find this to be an intriguing suggestion and we say so in the\nDiscussion. In the current work, OrthologTransformer’s mutation decisions are entirely data-driven\n(learning from sequences alone). The reviewer’s idea of incorporating structural predictions (e.g., using\nAlphaFold to evaluate a candidate mutation’s structural impact during generation) could indeed further\nensure functional conservation. We have added a forward-looking statement: “At the modeling level, we\nwill explore explicit evolutionary‑distance conditioning (e.g., a distance token or embedding) to\ncalibrate edit propensities across divergence scales, and we aim to integrate structure prediction or functional\nscoring (e.g., AlphaFold‑based rescoring) into the design loop to guide non‑synonymous edits –\nparticularly for more distantly related conversions where functional preservation is most challenging.”\nComment 10 (Balancing sequence adaptation with functional conservation – what features are\nprioritized by the model during generation?): Reviewer #1 asked how the model balances making changes\nfor host adaptation with preserving protein function. In other words, what features (e.g., conserved domains,\nstructural stability) does the model prioritize to ensure functionality is maintained while the sequence is\naltered?\nResponse: This question goes to the heart of how OrthologTransformer operates, and we have expanded our\nDiscussion to give insight into this “black box”. While the model does not explicitly encode notions like\n“conserved domain” or “structural stability” as separate variables, these aspects are implicitly learned through\nthe training on ortholog pairs. We now explain that:\n• Nature-defined trade-off (ortholog supervision, not manual prioritization). OrthologTransformer is trained\nin a supervised seq2seq manner on naturally occurring ortholog pairs that already embody the balance\nbetween host adaptation and function preservation. Consequently, the model does not decide this balance\nvia an explicit objective or hand-tuned weights; rather, it reproduces the patterns present in orthologs\n8\n(synonymous changes, conservative amino-acid substitutions, occasional indels). This point is now stated\nexplicitly in the revised text.\n• Connection to evolutionary selection (added analysis). Motivated by Reviewer #5’s suggestion, we\nconducted a site-wise dN/dS analysis on homologs to contextualize where OrthologTransformer proposes\namino-acid substitutions. Consistent with an ortholog-driven learning signal, model-edited sites were\nenriched in positions showing higher dN (less purifying constraint), whereas highly conserved (low-dN)\nsites were less changed. We have incorporated these results and their interpretation into the revised\nDiscussion to clarify that the model’s edits track evolutionary permissiveness rather than overriding\nfunctional constraints.\nIn summary, OrthologTransformer learns the adaptation–conservation balance from nature: its outputs\npreserve functional cores while adjusting DNA‑level features for the target host in a manner consistent with\nhow orthologs differ across species. The added dN/dS analysis makes this mechanism explicit in the revision.\nComment (Remarks on code availability): Reviewer #1 asked to specify key aspects of the computational\nenvironment, particularly the minimum GPU requirements (e.g., CUDA version, minimum GPU memory,\nGPUmodel compatibility). They also noted that given the model's size, users without sufficient GPU resources\nmay encounter difficulties reproducing the results..\nResponse: We appreciate the reviewer’s comment regarding the documentation and computational\nenvironment requirements. To address this concern, we have implemented the following updates:\n• Documentation of Environment: We have revised the documentation to clearly specify the required\ncomputational environment, including the necessary CUDA version, minimum GPU memory, and\ncompatible GPU models. This information ensures that users understand the hardware and software\nrequirements for running our model.\n• Pretrained Model Availability: The pretrained model checkpoint is now publicly available. This allows\nusers to access the model weights directly, which is particularly helpful for those who cannot undertake\nfull model training due to resource constraints.\n• As we mentioned above, our model’s size is moderate and smaller than large protein generative models\nsuch as ProGen and ESM, making training feasible on available hardware. In inference,\nOrthologTransformer can generate a redesigned gene in a matter of seconds on a single GPU, which is a\none-time cost per gene.\nReviewer #2\nComment 1 (Limited model innovation – use of standard Transformer architecture): Reviewer #2 felt\nthat the methodological innovation of our model appears limited, noting that we largely adopted a standard\nsequence-to-sequence Transformer with only minor modifications (special tokens for species). They argued\nthat adding species tokens is not a substantial advancement over the original Transformer architecture from\nmachine translation.\nResponse: We respectfully clarify that while the core architecture is indeed a standard sequence-to-sequence\nTransformer, the novelty lies in (i) the new problem formulation – orthologous gene conversion at the DNA\nlevel framed as translation between species – and (ii) tailoring the model and data pipeline to that problem,\nrather than inventing a new neural block. In the revision, we explicitly articulate these contributions in the\nIntroduction and Methods, and link them to empirical gains.\n• Problem formulation (conceptual advance). We cast cross-species gene redesign as a seq2seq “translation”\ntask on coding DNA, training on many ortholog pairs so the model learns evolution-consistent edits (from\nsynonymous swaps to amino-acid changes and indels) that preserve function in the target species. To our\nknowledge, applying a Transformer in this DNA-level ortholog conversion setting has not been previously\nexplored, and we now make this point explicit.\n• Species-conditioning (directionality and multi-species generalization). We prepend source/target species\ntokens to condition the model on the conversion direction, enabling one model to serve many species pairs\nwhile respecting their distinct codon biases and sequence preferences (Figure 1). This is essential to\ngeneralize across diverse evolutionary distances without training separate models per pair; we now\nemphasize this rationale in the text.\n• Two-stage learning (generality + precision). We combine many-to-many pretraining (to learn broad\ncross-species rules) with pair-specific fine-tuning (to capture nuances of a given source to target pair),\nas shown in Figure 1b. Figure 2 further illustrates accuracy improvements when moving from limited to\nbroad training regimes and when applying fine-tuning, motivating this design choice in the revised\nMethods.\n• Allowance of non-synonymous substitutions and indels (task-critical capability). Unlike synonym-only\ncodon optimizers, our autoregressive decoder is trained to output full target coding sequences – so length\nchanges and amino-acid changing edits arise naturally when supported by ortholog supervision. Figure 4a\ndocuments the resulting edit spectra per variant (insertions/deletions/synonymous/non-synonymous),\nillustrating that the model can introduce conservative amino-acid changes while preserving structure. We\nclarify this motivation and implementation in the Method and Results sections.\n• Alignment-aware supervision to expose biologically relevant edits. To explicitly teach the model where\nevolution edits proteins, we prepared alignment-processed and unprocessed datasets so the model sees\nexplicit patterns of non-synonymous changes and indels at corresponding positions. In the conversion\nfrom I. sakaiensis to B. subtilis task, alignment-based training sometimes improved conversion accuracy\nrelative to alignment-free training (Figure 2, gray vs black bars), supporting this choice; we now state this\n10\nrationale explicitly.\n• Additionally, we guided the model with multi-objective evaluation (Monte Carlo Tree Search (MCTS)\nroutine that jointly optimizes GC content and mRNA secondary structure stability) during candidate\nselection. While these might not constitute modifications to the Transformer architecture itself, they are\nsignificant at the application level that were essential for success in this domain.\nIn sum, our contribution is not the invention of a new Transformer variant, but a new, biology-grounded\nmodeling task and a carefully adapted training/data strategy that equips a standard Transformer to perform\ncross-species DNA sequence conversion – a capability that codon-only methods cannot offer. We have revised\nthe Introduction, Methods, and Discussion sections to make these design motivations explicit and to connect\nthem to the empirical improvements reported in Results section.\nComment 2 (Insufficient methodological comparisons – need for rigorous benchmarking): Reviewer #2\nfelt that the manuscript lacked quantitative comparisons to traditional approaches. They requested more\nstatistically rigorous benchmarking against conventional methods, with clear evaluation metrics, to\nconvincingly demonstrate the model’s advancement.\nResponse: We agree that thorough benchmarking is critical. In our original submission we already reported\nquantitative comparisons on two sequence-identity metrics – codon-level identity and amino-acid–level\nsimilarity – across multiple representative species pairs spanning diverse genomic or physiological\ncharacteristics, using large gene sets per pair (around several thousands of orthologous genes per pair);\nthese results are summarized in Table 1 and the associated Results text. The benchmarking was performed on\nrigorously constructed training and test splits to avoid any data leakage. We partitioned ortholog groups such\nthat no ortholog pair in the test set shares a gene with any pair in the training set. In response to Reviewer #2\nrequest, we conducted statistically rigorous test against conventional methods with three clearly defined\nevaluation metrics: codon-level sequence identity, CAI proximity and GC alignment between source, target,\nand model-generated sequences.\nIn the revised manuscript, we have strengthened this benchmarking:\n We conduct statistical analyses to test differences across methods, and we now report statistical\nsignificance of improvements in the Results (in a new Table 2).\n We added CAI and GC content as formal evaluation metrics alongside codon-level identity (in Table 2).\nAs context, our existing Results already indicate that OrthologTransformer substantially outperforms\ntraditional synonym-only approaches: we observe ~1.7-fold average improvement (with statistical\nsignificance p < 10−5) in codon-sequence identity over frequency-based codon optimization across species\npairs, and ~1.8-fold (p < 10−5) in the subset comparable to CodonTransformer (Table 1 and 2). In the revision,\nall three metrics – codon-level identity, CAI proximity, and GC‑content alignment – show statistically\nsignificant gains for OrthologTransformer relative to the baselines (significance values reported in the Results\nand Table 2).\nReviewer #2 requested to include a detailed description of our tokenization strategy in the main text. They\nalso questioned our use of AlphaFold-predicted TM-score as the measure of structural similarity, suggesting\nthat using structural alignment (e.g., RMSD from experimentally determined structures) might be more\nappropriate.\nResponse: We have addressed both of these points by adding clarifications to the manuscript:\n• In the Methods (OrthologTransformer Architecture), we now explicitly describe that sequences are\ntokenized at the codon level (each three‑nucleotide codon is one token, with start/stop as special tokens).\nWe also detail the species‑conditioning tokens – e.g., a [Source_species] token at the encoder input and a\n[Target_species] token at the decoder input – to indicate the conversion direction. The vocabulary thus\ncomprises 64 sense codons + start/stop + N species tokens. Each token (codon or species) is mapped to a\nlearnable 512‑dimensional embedding, initialized and trained end‑to‑end with the model. This section\nclarifies precisely how the DNA sequences and species context are represented for the Transformer.\n• Use of TM-score vs. RMSD: For structural comparison we used TM‑score on AlphaFold‑predicted\nstructures, as stated in the Evaluation Metrics, because experimental structures are not available for all\ndesigned variants and TM‑score provides a length‑normalized, fold‑level similarity measure that is less\nsensitive to local deviations than RMSD – appropriate when comparing predicted models across\nsequences with small edits. In response to the reviewer’s suggestion, the revised manuscript also reports\nRMSD and per‑residue pLDDT to complement TM‑score (in new Figure 4b and 4c). These additions\nprovide both a global (TM‑score) and local (RMSD/pLDDT) view of structural conservation.\nComment (Specific comments on the PETase case study): Reviewer #2 provided several specific comments\nrelated to the PETase experiment and its presentation.\nSpecific Comment 1 (Abstract claim – need quantitative validation of “surpassing codon-optimized\ncontrols): Reviewer #2 noted that in the Abstract we stated the redesigned PETase showed “robust activity\nsurpassing codon-optimized controls,” and they suggested we should quantify this improvement (e.g., how\nmany fold improvement) to substantiate the claim..\nResponse: We have revised the Abstract to include quantitative details of the improvement. It now reads (in\nsummary): “As proof of concept, an OrthologTransformer-redesigned PETase expressed in Bacillus subtilis\nshowed robust activity, producing approximately 10-fold more reaction product than the codon-optimized\nenzyme, and achieving higher expression levels, thereby demonstrating improved enzyme performance via\nAI-guided gene design.” By explicitly stating the fold-improvement (10× increase in MHET production in our\nPET degradation assay), we satisfy the reviewer’s request.\nSpecific Comment 2 (Validation in Gram-negative or eukaryotic hosts): Reviewer #2 The reviewer asked\nwhether we have tested the model’s performance in a Gram-negative bacterium or a eukaryotic system, given\ntheir distinct codon usage and regulatory elements. They noted our current experiments were only in B.\n12\nsubtilis (Gram-positive).\nResponse: We acknowledge this important point. As mentioned in Response to Reviewer #1, Comment 1\nabove, we have taken steps to demonstrate the model’s applicability to other host types. In the revised\nmanuscript, we report new experimental results for a Gram-negative scenario. We applied\nOrthologTransformer to redesign a gene for E. coli expression (Gram-negative host) and compared it with a\ntraditional codon-optimized version. Preliminary results indicate that the OrthologTransformer-designed gene\nyields higher protein expression and PET degradation activity than the codon-optimized gene. For a\neukaryotic system, we have not yet completed a wet-lab validation due to the complexity of eukaryotic\nexpression systems. Designing for a eukaryotic host will likely require coupling CDS optimization with\nseparate consideration of promoters, UTRs, and potential alternative splicing. This remains an area for\nfuture work.\nSpecific Comment 3 (HPLC data for PETase assay – provide chromatograms for TPA/BHET and\ndiscuss MHET focus): Reviewer #2 pointed out that in our Methods we mentioned measuring multiple PET\ndegradation products (TPA, MHET, BHET) by HPLC (Supplementary Figure S5), but the materials provided\n(Figure 5 and Supplementary) seemed to show only MHET quantification and a pathway diagram. They\nrequested that we include the actual HPLC chromatograms, standard curves for product detection, and\nexplain why only MHET was shown as the readout. Essentially, they want assurance that the assay was\nsensitive to TPA/BHET and a rationale for focusing on MHET.\nResponse: We have updated the Supplementary information to include the requested data and have expanded\nthe discussion of the degradation assay in the main text, as also mentioned in Response to Reviewer #1,\nComment 7. Specifically, we added Supplementary Table S5, which shows example HPLC chromatograms for\nMHET, TPA, and BHET standards and for sample supernatants from our PETase experiments. These additions\ndemonstrate the sensitivity of our HPLC method for all three products.\nIn the Results section, we now clarify why MHET was the primary quantified product. In our experimental\nsetup (PET film degradation by PETase in B. subtilis culture supernatants), we consistently detected MHET\naccumulation, but TPA and BHET remained below detectable level. We interpret this as follows: MHET is the\nmono-hydrolysis product of PET (one ester bond cleaved), whereas TPA would result from complete\ndegradation (both ester bonds cleaved). B. subtilis does not naturally have a MHET hydrolase to further break\ndown MHET into TPA, so MHET accumulates. BHET is an intermediate di-ester that is less commonly\nproduced except under certain conditions. We explain that, given the absence of additional enzymes, the\nreaction stalls at MHET, making it the most relevant readout of PETase activity in our system. We have added\na note that measuring MHET is sufficient to gauge PETase performance, since any TPA produced would\nindicate an even further extent of reaction. Importantly, we now explicitly state that no significant TPA was\ndetected in our assay, confirming that our redesigned PETase, like wild-type PETase, primarily generates\nMHET.\nWe also mention in the Discussion that future work could incorporate a downstream MHETase to fully\nconvert MHET to TPA, allowing a closed-loop plastic degradation system.\nReviewer #5\nMajor Comments:\nComment 1 (Evolutionary distance not explicitly modeled – potential limitation): Reviewer #5 pointed\nout that our model uses a large dataset of orthologs from species with varying evolutionary distances, but we\ndid not explicitly include evolutionary distance as an input or feature. They ask if not modeling evolutionary\ndistance could be a limitation, and suggest that including such information might help the model learn how\northologs evolve across different divergence levels.\nResponse: This is a thoughtful observation. In our current model, evolutionary distance is indeed not\nexplicitly given to the model (we do not, for example, feed a numeric distance or category). Instead, the model\nmust infer the “distance” indirectly from the sequences themselves and the species tokens. We agree that\nexplicitly incorporating a measure of divergence could potentially enhance performance by contextually\nbiasing the model’s propensity to make changes (e.g., a large evolutionary distance might warrant more\nnon-synonymous substitutions). We have added a sentence in the Discussion as future work: “At the modeling\nlevel, we will explore explicit evolutionary‑distance conditioning (e.g., a distance token or embedding) to\ncalibrate edit propensities across divergence scales, and we aim to integrate structure prediction or functional\nscoring (e.g., AlphaFold‑based rescoring) into the design loop to guide non‑synonymous edits –\nparticularly for more distantly related conversions where functional preservation is most challenging.”\nComment 2 (Most successful variants had few non-synonymous mutations – discuss implications):\nReviewer #5 observed that some of the best-performing redesigned PETase variants (notably AI9 and AI12)\ncontained few or no non-synonymous mutations, even though our approach emphasizes the ability to\nintroduce such mutations. They note that AI9, which had only a single amino acid change, was the top\nperformer (3× more MHET production, most PET surface erosion). This raises the question: does the success\nof low non-synonymous-change variants undermine the value of modeling non-synonymous/INDEL changes?\nThey ask us to discuss why the most powerful redesigned proteins in this case were those with minimal amino\nacid changes, and what that means for our approach that also explores non-synonymous changes.\nResponse: The reviewer brings up an interesting paradox: if the best designs ended up being close to the\noriginal protein sequence, one might wonder if our allowance for non-synonymous mutations was necessary.\nWe address this by analyzing why these variants succeeded and how we interpret the outcome. We appreciate\nthe observation and agree that the PETase case highlights an important point: enhanced activity in our system\ndid not require extensive amino‑acid changes. The result for AI9 (renamed AI‑L2 in the revision) does not\nundermine the value of enabling non‑synonymous substitutions and indels. Trained on natural ortholog pairs,\nOrthologTransformer learns when evolution typically leaves protein cores unchanged and when amino‑acid\nadaptations occur. As noted in the Discussion, for the S290N substitution introduced in AI‑L2, two prior\nstudies on protein engineering reported a change at the same residue (S290P) that enhanced PETase stability\nand activity, supporting the view that the model introduces amino‑acid/indel edits where beneficial. Moreover,\nacross species pairs, our benchmarking shows that sequences generated by OrthologTransformer (which may\n14\ninclude non‑synonymous substitutions) more closely resemble target‑species orthologs and improve\nfunctional similarity compared with synonym‑only approaches; this includes gains in codon‑sequence identity\nand BLOSUM‑based amino‑acid similarity, indicating conservative, biologically informed non‑synonymous\nsubstitutions when advantageous.\nComment 3 (Clarity on selection of the 12 final designs (AI1–AI12)): Reviewer #5 felt it was not very clear\nhow the 12 PETase variants were selected. They ask us to clarify what parameters or combinations were\nconsidered and whether there was a larger initial pool of candidates. If so, why were certain designs chosen\nor others excluded from experimental validation?\nResponse: We agree that the selection pipeline for the 12 PETase variants was not fully explained. In the\nrevision, we have clarified the process in the Methods/Results and aligned the presentation with a similar\nrequest in Reviewer #1, Comment 6. Concretely:\n• Moved and expanded table: The former Supplementary Table S2 (design “recipes” for AI1–AI12\n(renamed to AI-S1–AI-L5 in the revision)) has been moved into the main text (as a new Table 3) and\nexpanded with a brief step-by-step narrative of how the models/designs were generated under\nsystematically varied conditions – training source (Orthofinder vs OMA), breadth of training species (23,\n54, or 2,138), alignment processing (±), pair-specific fine-tuning (±), MCTS (Monte Carlo Tree Search\nOptimization) (±) to balance GC and RNA MFE – and that AI-S1–AI-L5 were selected as representative\noutputs spanning these axes.\nComment 4 (Figure 1A clarity – input/output depiction): Reviewer #5 suggested that Figure 1A (the\nschematic of the OrthologTransformer model) could be graphically improved for clarity. They found it not\nimmediately clear what the input and output of the model are from that figure, and they recommended\nelaborating the figure caption.\nResponse: We have revised Figure 1a in line with this feedback. The new version of the figure has a more\nstraightforward layout: we explicitly label the input sequence (with a notation like “Gene from Species A”)\nand the output sequence (“Orthologous gene for Species B”).\nThe figure caption has been expanded accordingly. It now clearly states: “Figure 1a: Input: a coding DNA\nsequence from Species A (source), prepended with a source_species token 𝑠 , is encoded. Output: the\n𝑡𝑔𝑡\ndecoder, conditioned by a target_species token 𝑠 , generates an orthologous coding sequence for Species B,\n𝑡𝑔𝑡\npermitting synonymous and conservative non‑synonymous substitutions and indels where supported by\northolog supervision. The model features a 20-layer encoder-decoder structure, with each layer equipped with\nAdd & Normalization layers and Multi-head Attention mechanisms. Species tokens (𝑠 , 𝑠 ) are prepended\n𝑠𝑟𝑐 𝑡𝑔𝑡\nto the input sequence, enabling species-specific sequence conversion… etc.” This caption now serves to\nreinforce what the input/output are and how the model operates on them.\nComment 5 (Embedding details not entirely clear): Reviewer #5 noted that while we mention using\ncodon-level tokens and an embedding dimension of 512 (lines 341-343 in the manuscript), we did not explain\n15\nhow the embedding was initialized or derived (one-hot vs learned, etc.). They asked for an explicit description\nof the embedding process.\nResponse: In the revised Methods (Model Architecture section), we have added a sentence to clarify the\nembedding. We now state: “Tokenization and embeddings: Input/output sequences are tokenized at codon\nresolution (three nucleotides per token). The vocabulary comprises 64 sense codons + start/stop + N species\ntokens, and each token is mapped to a 512‑dimensional learnable embedding. Embeddings are randomly\ninitialized and optimized end‑to‑end during training. The context length is 700 tokens, which corresponds to\n~2.1 kbp when using codon tokens, since each token is 3 bases.”\nComment 6 (Methods details missing – structure prediction and RNA folding): Reviewer #5 asked how\nthe protein structures were predicted for Figure 4a, noting that this was missing from the Methods. They also\nnoticed that we mention mRNA folding kinetics and the ViennaRNA package in our code, but did not describe\nthis in the Methods.\nResponse: We have corrected these omissions by updating the Methods section accordingly:\n• Protein structure prediction: We now explicitly state that we used AlphaFold3 (ver.3.0.1) to predict the 3D\nstructures of the redesigned PETase variants, as well as the wild-type and codon-optimized PETase, for\nstructural comparisons. We outline that we used the AlphaFold3 default pipeline with no templates (since\nwe wanted a prediction purely from sequence) and ran it for each sequence. We mention that we took the\ntop-ranked predicted structure for each and calculated TM-scores between each variant’s structure and the\nwild-type PETase structure.\n• mRNA structure and ViennaRNA: We now mention in the Methods that we evaluated mRNA secondary\nstructure for each designed sequence using the ViennaRNA package ver.2.6.4 (RNAfold). We specify that\nwe computed the minimum free energy (MFE) structure for the full mRNA sequence and recorded the\nMFE value as a measure of mRNA folding strength. This was done to compare whether our designs had\npotentially more favorable (less stable, thus more easily unwound) mRNA structures than the originals.\nMinor Comment:\nComment 7 (Relation to dN/dS (selection analysis)): Reviewer #5 made a comment out of curiosity (and\nnoted it’s not necessarily for the manuscript) about a methodological similarity: our approach reminds them\nof how one detects sites under positive/negative selection using dN/dS ratios. They wondered if the codon sites\nchanged by our model correlate with sites that would be identified as under selection.\nResponse: We find this insight very interesting and have decided to address it, even if not required, as it adds\nan extra perspective to our work. We have included this dN/dS analysis in the Results (and provided details in\nnew Table 4). In the Discussion (likely where we talk about how the model chooses mutations, related to\nReviewer#1’s comment 3 and 10), we include a short remark connecting to dN/dS analysis:\nWe say that indeed, OrthologTransformer’s behavior is analogous in some ways to a dN/dS-informed\napproach – positions that the model changes tend to be those that evolution (across species) has changed (i.e.,\nhigher dN), whereas conserved positions (low dN) are left unchanged by the model. We actually performed a\n16\nsite-wise dN/dS analysis on alignments of homologous sequences using the FEL method. For each site, the\nnonsynonymous (dN) and synonymous (dS) substitution rates were estimated; sites showing significant\nevidence of purifying selection (p < 0.01 and dN/dS < 1) were classified as “purifying”, whereas sites with\ndN/dS ≥ 1 or without significant purifying signal were classified as “non-purifying”. We found that\nnon-purifying sites exhibited significantly higher substitution rates than purifying sites (Fisher’s exact test, p <\n10−14). These findings indicate that OrthologTransformer preferentially introduces substitutions at\nevolutionarily variable sites, while changes at highly conserved positions are comparatively suppressed. We\nmention this observation to the reviewer to satisfy their curiosity.\nComment (Remarks on code availability): Reviewer #5 made the following remarks: We were able to install\nand run OrthologTransformer inside a Linux Docker container on a MacBook. The documentation on GitHub\nis great, but we noticed a few things you could improve: (a) make a release version of the code so that your\nanalyses are reproducible even when changing something in the future; (b) matplotlib is missing from\nrequirements.txt; and (c) inform the user that the results output folder should exist; we got an error at first\nbecause of that.\nResponse: Thank you very much for your valuable suggestions. We have addressed all three points as\nfollows:\na) Release version for reproducibility:\nTo ensure reproducibility of the analyses, we have created a fixed release version (v1.0.0) on GitHub. This\nallows users to reproduce the results even if the codebase evolves in the future.\nb) Added matplotlib to requirements.txt:\nAs suggested, we have added matplotlib==3.9.1 to the requirements.txt.\nc) Improved handling and documentation of the results output folder:\nWe have modified the code so that the results output folder is automatically created if it does not already exist.\nAdditionally, we updated the README to clearly inform users about this behavior to avoid confusion or\nerrors during the first run.\nhas been significantly improved in terms of clarity, completeness of results, and adherence to the journal’s\nrequirements.\nOnce again, we thank the reviewers and the editor for their valuable feedback, which has helped us improve\nour work. We hope that our revisions satisfactorily address all comments and that the updated manuscript is\nnow suitable for publication in Nature Communications.\nBest regards,\nYasubumi Sakakibara\nWe appreciate the reviewers’ further comments on our manuscript entitled “Cross-species gene redesign\nleveraging ortholog information and generative modeling.” We have carefully revised the manuscript to\naddress every point raised. Below, we provide a point-by-point response to each comment. All corresponding\nchanges in the manuscript are highlighted in red in the revised Word document.\nResponse to Reviewer #1\nComment 1: Reviewer #1 was concerned that while the added E. coli experiment improves the manuscript,\nboth validation systems (B. subtilis and E. coli) are prokaryotic, and both tests target a single enzyme\n(PETase). The authors claim that OrthologTransformer is broadly applicable to “any gene and any species,”\nyet this remains untested.\nResponse: We fully recognise that Comment 1 exposed an important issue with the way we described the\nscope of OrthologTransformer. We do acknowledge that the phrase “any gene and any species” was an\noverstatement, which may have led to unrealistic expectations, especially concerning eukaryotic applications.\nPlease note that in the manuscript itself, we did not claim universal applicability; instead, we focus on\nbacterial orthologs. We therefore no longer use that phrase and similar universal language in the response\nletter and clarified in the Discussion that our current study focuses on bacterial orthologs, and that the model\nwas trained and validated on prokaryotic genes. In addition, we note that our updated computational\nbenchmarking – the other main component of our work – already demonstrates superior performance, with\nclear improvements over existing codon‑optimization and deep learning method. In this second-round revision,\nwe conducted additional, more extensive computational benchmarking, and the results are presented in a new\nFigure 3. This new computational benchmarking result extends the performance comparison to the maximum\npossible scale, covering all 45 bacterial species included in both the OMA database (2,138 species) and the set\nof species that the existing deep-learning method CodonTransformer can handle (164 species), as well as 450\nsource-to-target species combinations (10 source species per target species).\nComment 2: Reviewer #1 was concerned that the authors report mRNA levels, protein secretion, and PET\ndegradation activity, but these datasets are not quantitatively correlated. It remains unclear whether the\nenhanced activity of AI-L2 results mainly from higher expression, improved folding, or kinetic enhancement.\nResponse: We agree that multiple interdependent factors underlie the enhanced activity of the AI-designed\nenzyme (AI-L2). In our previous revision, we emphasized that no single molecular property fully explains the\nperformance gain. Rather, several jointly optimized features – including codon usage adaptation, improved\nmRNA stability/structure, and possibly beneficial amino acid substitutions, which can affect mRNA\nexpression, protein folding, and enzyme kinetics – are likely to act together in a synergistic, non-linear manner\nto yield the higher PETase activity in B. subtilis. We have added a brief explanation to the Result highlighting\nthat PETase performance in vivo is governed by multiple interacting, non-linear factors – including mRNA\n1\nabundance/stability/structure, translation efficiency/codon adaptation, and beneficial amino acid substitutions,\nwhich in turn can affect mRNA expression, protein folding, and enzyme kinetics – rather than by any single\ndeterminant.\nComment 3: Reviewer #1 pointed out that AI-L2 (formerly AI9) has a higher kcat/Km but provide little\nstructural rationale for this effect. Without mapping mutation sites or discussing potential impacts on the\nactive site or enzyme stability, the functional link remains speculative.\nResponse: We did include an analysis of a key mutation in the previous manuscript, and we apologize if it\nwas not sufficiently highlighted. In the Discussion, we explicitly discussed the substitution S290N in AI-L2 as\na potential structural/functional contributor. To quote the relevant sentence from the previous manuscript:\n“Regarding the mutation S290N introduced in AI-L2, two earlier studies have reported a substitution at the\nsame residue position, specifically S290P. Serine (S) and asparagine (N) are both neutral amino acids with\npolar side chains, whereas proline (P) is characterized by a nonpolar side chain. Thus, the AI-L2 coding\nsequence included only a handful of amino acid substitutions that the model likely inferred from related\nBacillus enzyme orthologs.” This discussion (now clearly visible in the revised Discussion section) provides a\nstructural rationale: in addition to prior studies showing that introducing mutations at the serine residue at\nposition 290 contributes to enhanced PETase enzymatic activity and thermostability, the S290N change is a\nsafe substitution relative to the wild-type serine (maintaining polarity), unlike the previously attempted S290P\nsubstitution which introduces a nonpolar residue.\nComment 4: Reviewer #1 noted that the paper introduces “species-conditioning tokens,” “two-stage\ntraining,” and “MCTS optimization,” but the technical specifics are only briefly described. For example, the\nnumber and encoding method of species tokens, the weighting in the MCTS reward function, and\nhyperparameter details are incomplete.\nResponse: We have augmented the Methods and Supplementary Information to address these technical points\nin detail. First, the concept of species-conditioning tokens was already explained in the Materials and Methods\n(“Tokenization and embeddings” subsection) in the previous manuscript. We clarified there that the model’s\nvocabulary includes N special tokens corresponding to each training species, which are prepended to the input\nand output sequences to condition the Transformer on source and target species. (For clarity, we have now\nexplicitly noted that N is the number of distinct species used in training and test.) Second, we have added a\nmore technical description of our two-stage training strategy in the Methods. This new text outlines how we\nfirst pretrain the model on a broad many-to-many orthologous gene pairs (any combination of source–target\nspecies) to capture general patterns, and then fine-tune it on a specific source–target species pair to\nincorporate pair-specific nuances. Third, we expanded the Methods (and Supplementary Methods) to describe\nthe MCTS optimization procedure. We detail how our Monte Carlo Tree Search algorithm was used to refine\nthe generated sequences, including the exact multi-objective reward function that balances GC content and\nmRNA secondary structure stability (melting energy). The weighting scheme used in the MCTS reward\nfunction and the corresponding hyperparameter details are now provided in the Supplementary Methods.\nallow readers to understand and, if desired, reproduce our approach.\nComment 5: Reviewer #1 was concerned that although the authors claim that training and test sets are\nstrictly separated, the manuscript does not examine whether certain bacterial clades dominate the training\ndata. Uneven representation could bias the model toward overrepresented taxa.\nResponse: We apologize for the confusion and have clarified the partitioning scheme in the Methods. Our\ndata partition was done at the gene (ortholog group) level, not at the species level. In practice, this means that\nif a particular gene (and its orthologs) was used in training, none of those genes or their orthologous\ncounterparts appear in the test set. There is no overlap of orthologous groups between training and testing. We\nexplicitly stated in the previous manuscript that no ortholog pair in the test set shares any gene with the\ntraining set. This strategy ensures that no systematic bias in the representation of bacterial clades is introduced\nin the training data by our partitioning strategy.\nComment 6: Reviewer #1 noted that the dN/dS analysis shows that model-introduced mutations tend to occur\nin variable regions, which supports biological plausibility, but it does not prove that the model “learns\nevolutionary rules.”\nResponse: We appreciate this concern and have revised the wording to ensure objectivity. We would like to\nclarify that our manuscript never claimed that the model “learns evolutionary rules” (indeed, that phrase does\nnot appear in our paper or prior response). To prevent any misunderstanding, we have rephrased some related\nexpressions in the revised manuscript. For example, instead of saying the model captures “general principles\nof protein adaptation,” we now state that it captures patterns of gene adaptation observed in nature. Similarly,\nwe replaced the term “evolution’s imprint” with a more neutral description of observed ortholog variation.\nThese edits avoid overreaching statements and simply report that the model reproduces the kinds of\nsubstitutions seen in naturally evolved orthologs (as supported by our dN/dS analysis).\nResponse to Reviewer #2\nComment: The author has revised the manuscript as requested and has addressed the reviewers' comments\nthoroughly. I do not have additional comments for it.\nResponse: We sincerely appreciate Reviewer #2’s careful evaluation and thoughtful feedback.\nResponse to Reviewer #5\nComment: Thank you very much for this comprehensive review and for addressing all of our questions and\nconcerns (even beyond what we expected by also integrating the quite interesting dN/dS part). We are satisfied\nwith the current state of the manuscript and look forward to its publication.\nResponse: We sincerely thank Reviewer #5 for their careful evaluation and positive recommendation.\nrevisions satisfactorily address all comments and that the updated manuscript is now suitable for publication\nin Nature Communications.\nBest regards,\nYasubumi Sakakibara\nWe appreciate the reviewers’ further comments on our manuscript entitled “Cross-species gene redesign\nleveraging ortholog information and generative modeling.” Below, we provide a point-by-point response to\neach comment.\nResponse to Reviewer #1\nComment: (Remarks to the Author) The revised manuscript is much improved, and the authors have\nadequately addressed the concerns raised in the previous round. The scope of the study is now clearly defined,\nthe methodological descriptions are more transparent, and the interpretations are appropriately cautious. I\nhave no further scientific concerns at this stage. The experimental and computational results support the\ncurrent claims, and the manuscript is suitable for publication.\nResponse: We sincerely thank Reviewer #1 for their careful evaluation and positive recommendation.\nWe hope that the updated manuscript is now suitable for publication in Nature Communications.\nBest regards,\nYasubumi Sakakibara\n1"}]

Decision letter