To search for the nucleotide sequences of garlic γ-glutamyl transpeptidase transcripts in the National Center for Biotechnology Information (NCBI) GenBank (26), the following query was used: ‘gamma-glutamyl transpeptidase’ AND ‘Allium sativum’[porgn:__txid4682]. The three resulting nucleotide sequences were stored in FASTA format. To identify and download their corresponding peptide sequences in FASTA format, the link of each NCBI GenBank entry was followed to its corresponding NCBI Protein (https://www.ncbi.nlm.nih.gov/protein/) entry.

The garlic genome sequence in FASTA format, as well as its corresponding annotation in general feature format (GFF) format (27), were downloaded from NCBI Genome (https://www.ncbi.nlm.nih.gov/genome/). The GenBank accession of the genome assembly used in the analysis in the present study was GCA_014155895.2(25). Due to the size of the Allium sativum genome chromosomes, it was deemed necessary to fragment them, as basic local alignment search tool (BLAST)+ (28) that was to be used to search for sequences of interest in the garlic genome cannot handle sequences longer than 1 Gbp. For this reason, a script in PHP programming language was created to split the chromosome sequences into sequences with maximum size 1 Gbp. The preliminary identification of the genomic regions of the three known genes was based on the pairwise alignments between their nucleotide and protein sequences, and those of the genome, as obtained from BLAST+ searches: The boundaries of the AsGGT exons were identified using BLASTN search for the nucleotide sequences. Similarly, TBLASTN searches for the peptide sequences revealed the genomic coordinates of the coding sequence (CDS) of each exon for every AsGGT gene. Exonerate (29) was then used for the precise mapping of the nucleotide and protein sequences on the genome. Exonerate determined the exact boundaries for each gene, exon, CDS and 5' and 3' untranslated regions (UTRs), and produced a GFF file for each gene. Finally, a manual inspection of these boundaries was performed using Integrative Genomics Viewer (IGV, Version: 2.16.1) (30). Exonerate and manual inspection ensured that the splice sites of the introns belong to the GU-AG group (31).

To compare the manual annotation for the three genes with the annotations produced automatically (25), all relevant genomic features from the automatic GFF file were extracted to create the GFF specifically for these genes. The automatically and manually created exon-intron boundaries were visually inspected using IGV and the protein sequences were extracted from the GFF files and the garlic genome with another Gtf/Gff analysis toolkit (AGAT) (32). Global pairwise alignments (33) between the automatically and manually produced protein sequences were performed.

A multiple sequence alignment of the three protein sequences was performed using Clustal Omega (34). It was noted that apart from the first exon CDSs which did not present a high degree of conservation among the three peptides, the boundaries of the CDSs of all other exons matched perfectly. Based on this finding, HMM profiles for the multiple sequence alignments of the six conserved CDSs were built using HMMER (35). All characterized Allium sativum proteins were collected from NCBI Protein. This search yielded 804 peptide sequences which were stored as a single FASTA file. A search for the HMM of each exon against all characterized Allium sativum proteins was performed. This search did not yield any new peptide sequence beyond the three already known proteins; therefore, it was assumed that no other characterized protein was homologous to the known AsGGTs.

With the aid of a PHP script, the translation of the whole Allium sativum genome into the six open reading frames (ORFs) was performed as previously described (36). With this procedure, all possible peptides, as well as their corresponding GFFs were generated. To search for homologous sequences potentially encoding AsGGT enzymes in addition to those already characterized ones, a PHP script searched the HMM of each exon against all potentially genome-coded peptides. This produced a GFF file per chromosome or scaffold which contained all genomic regions that could code for AsGGT CDSs. To identify chromosomes or scaffolds coding AsGGT genes, a manual check of their corresponding GFF files was performed. The GFF files for chromosome 8 (CM031537.1), chromosome 5 (CM031533.1) and chromosome 4 (CM031532.1) contained the CDSs of the already characterized AsGGT1, AsGGT2 and AsGGT3 genes, respectively. The GFF file for chromosome 6 (CM031534.1) contained consecutive CDSs which corresponded to exons 2-7, suggesting the existence of a near complete gene structure. Mapping the three known protein sequences on chromosome 6 with Exonerate and consequently performing manual inspection with IGV, a GFF file containing genomic features (exons 2-7) of this candidate gene was created.

The protein sequence of the newly discovered AsGGT-like (AsGGT4) gene was extracted from its GFF and the Allium sativum genome sequence, using AGAT. A BLASTP search limited to Liliopsida (monocots) was performed in UniProt (37) to examine whether the peptide in question had already been identified and to identify homologous peptide sequences in related species of Allium sativum. The peptide had not been previously identified and the UniProt BLASTP search revealed that its amino terminal end was occasionally aligned with the signal peptide of homologous proteins, suggesting that the potential enzyme may contain a signal peptide at its amino terminal end. Thus, SignalP 6.0(38) and DeepTMHMM (39) were used to check for the existence of a eukaryotic signal peptide in the known AsGGT enzymes, as well as the potential one.

The confirmation of the existence of a signal peptide in the new enzyme allowed for the prediction of the length of the CDS region of its missing first exon: As signal peptides usually have a length of 22-25 amino acids and 19 amino acids are already found in the CDS region of exon 2 in this enzyme, it was expected that the CDS region of the first exon would encode ~4 amino acids, beginning with a methionine (AUG codon). In addition, the boundaries of the intron between exon 1 and exon 2 of the gene might be of GU-AG type, as explained above. Visually scanning the genomic area upstream exon 2 in IGV, a sequence matching the above criteria was discovered, thus completing the potential gene sequence. GFF features of the genomic region where AsGGT4 gene is located, were extracted from the automatically produced GFF file (25). A comparison between the automatically and manually produced annotations as they were depicted in their respective GFFs, was performed by visual inspection in IGV.

To ensure that AsGGT4 is actually transcribed, an alignment of all available Allium sativum short RNA-seq reads was performed on its potential transcript, using a PHP script. As a positive control, the same alignment was performed on the AsGGT3 transcript. AsGGT3 was selected as it was found that it was the most similar known AsGGT to AsGGT4. SRA (40) was searched for all Illumina RNA-Seq runs of Allium sativum samples. SRR IDs of these entries were used for the download of their corresponding fastq.gz files from the European Nucleotide Archive (ENA) (41). These files were then unzipped and converted into FASTA files using EMBOSS (42). Each FASTA file was converted into a BLAST database and the two transcript sequences were BLASTN-searched against the RNA-Seq read database. The outputs of each search were then manually examined to identify reads which could correspond to the AsGGT3 or AsGGT4 transcripts.

As no solved 3D structures for AsGGT proteins are available, AlphaFold (43) was used for structure predictions from primary peptide sequences. The predicted structures of the known AsGGTs were downloaded from the AlphaFold Protein Structure Database (44). A prediction of the tertiary structure of the AsGGT4 was performed within the ColabFold website (45), which uses a simplified version of AlphaFold v2.3.2. For AsGGT4 structure prediction, homology modelling was also performed using the SWISS-MODEL Workspace (46).

Since SWISS-MODEL used the solved structure of Bacillus licheniformis γ-glutamyl-transpeptidase [Protein Data Bank (PDB) ID: 4Y23] to predict the structure of the AsGGT4 peptide by homologous modelling, the bacterial peptide sequence was used as a distant evolutionary relative (outgroup) of Allium sativum AsGGT peptides. This sequence was downloaded from PDB (47) and was then searched using BLASTP in UniProt. Its corresponding sequence in UniProt ID was Q65KZ6_BACLD. MUSCLE (48) was used to generate a multiple sequence alignment and a phylogenetic tree of the peptide sequences of the four AsGGTs and their distant relative, Bacillus licheniformis γ-glutamyl transpeptidase. Visualization of the multiple sequence alignment was performed with JalView (49). The phylogenetic tree of the five peptides in Newick format (50) was visualized with Dendroscope (51).

The constructed GFF files for the AsGGT1, AsGGT2 and AsGGT3 genes were visualized (Fig. 2A-C). The analysis revealed that AsGGT1 is located on chromosome 8: 403,215,661-403,234,269, AsGGT2 on chromosome 5: 617,284,828-617,294,697 and AsGGT3 on chromosome 4: 182,866,703-182,871,895. All three genes consist of seven exons of comparable size. Although the size of each corresponding exon is comparable among the three genes, the size of the corresponding introns varies, resulting in a considerable difference in gene size: The AsGGT1 gene size is 18,609 bp, that of AsGGT2 is 9,870 bp and that of AsGGT3 is 5,193 bp.

Through a manual inspection of the genomic features of AsGGT3 on chromosome 4, it became apparent that there are four loci in the genomic sequence where single base deletions or insertions (indel events) occurred (Fig. 3). An addition of an adenine appears at position 182,867,809 in the second exon, an addition of a guanine appears at position 182,869,779 in the third exon, a deletion of an adenine appears at position 182,870,778 in the fourth exon and an addition of a guanine appears at position 182,870,929 in the seventh exon.

Search for the HMMs of the six coding sequences which correspond to already characterized AsGGT exons in the potential proteome, which resulted from the in silico translation of the genome into the six ORFs, in addition to identifying the coding regions of known genes, allowed for the discovery of novel genomic regions that could encode parts of γ-glutamyl transpeptidase enzymes.

By filtering the search results and keeping only those containing adjacent regions corresponding to all exons that comprise a potential AsGGT gene, a crude GFF was generated that described the coding regions of a potential fourth γ-glutamyl transpeptidase gene (AsGGT4). The complete GFF file for AsGGT4 was constructed by optimizing the primary GFF. AsGGT4 is located on chromosome 6: 92,586,715-92,589,081 and has a size of 2,336 bp. Similar to the three already characterized genes (AsGGT1, AsGGT2 and AsGGT3), the AsGGT4 (Figs. 2D and 4) gene consists of seven exons. Unlike the already characterized genes, AsGGT4 is transcribed in the reverse direction.

A BLASTP search in UniProt for homologous monocot protein sequences with the potential γ-glutamyl-transpeptidase sequence revealed a number of sequences belonging to different species (Fig. 5). The visual observation of the produced pairwise alignments revealed that a number of homologous peptides possessed a signal peptide. It was therefore examined whether the three known enzymes and the potential one, contain a signal peptide in their amino terminal sequence using SignalP and DeepTMHMM. SignalP did not predict a signal peptide for AsGGT1 or AsGGT2 peptides. DeepTMHMM predicted a single transmembrane helix in 32-47 and 30-46 amino acid regions of AsGGT1 and AsGGT2, respectively. It also predicted that their C-terminus is extracellular. SignalP revealed a potential signal peptide profile, mainly in the 23-46 amino acid region of AsGGT3 (Fig. 6A), beginning with a methionine at position 23. When the first 22 amino acids were removed from the sequence, SignalP revealed the presence of a signal peptide, where amino acids 1-7 form its amino terminal region, amino acids 8-18 are hydrophobic and finally amino acids 19-24 form its carboxy terminal end (Fig. 6B). DeepTMHMM also predicted a signal peptide for both the untruncated and truncated AsGGT3 peptide.

SignalP predicted a 19 amino acid long signal peptide in the amino terminal coding region of the second exon (Fig. 7A), implying that the coding sequence of the first exon of AsGGT4 would code for approximately another four amino acids. By manually inspecting the genomic region upstream of the second exon, a first exon which could code for four amino acids was proposed. SignalP predicted a signal peptide of 23 amino acids in the full length AsGGT4 peptide, where the first four amino acids form its amino terminal region, amino acids 5-18 are hydrophobic, while amino acids 19-23 form its carboxy terminus (Fig. 7B). DeepTMHMM also predicted a signal peptide.

The coding regions of AsGGT1 and AsGGT2 are identical to the coding regions obtained by automated genomic annotation; i.e., from the use of the GFF file that was provided by the group that performed the genome sequencing (Fig. 8A and B); by contrast, automated genomic annotation failed to correctly assign the exon boundaries and coding regions of AsGGT3 (Fig. 8C). The extracted protein sequence also differed from the one previously reported (24).

Pairwise alignments between the protein sequences extracted from manual and automatic annotation revealed that the peptide sequences of the AsGGT1 and AsGGT2 genes from the automatic annotation matched almost perfectly with the sequences obtained experimentally and with the annotation performed in the present study-small differences may be due to point mutations. By contrast, marked differences were observed between the automatic annotation of the peptides of the AsGGT3 and AsGGT4 genes and the annotation performed the present study. In AsGGT3, differences appear at the four points where indel events occur in the genomic region, as expected, but the exon boundaries are generally correctly described beyond these points. By contrast, the automatic annotation of AsGGT4 appears to have failed almost completely. Although the exon boundaries are correctly predicted (Fig. 8D), the ORF is incorrect in all exons, and the first exon is completely missing-the start of the code region is incorrectly located in the second exon. Thus, pairwise alignment does not detect any similarity-apart from a minute region that is probably accidental, between the two extracted peptide sequences (Fig. 9). By changing the phase from 2 to 1 in the automatically produced GFF, the resulting peptide sequence is the correct one.

A multiple sequence alignment (Fig. 10) of the four garlic AsGGT peptide sequences and the Bacillus licheniformis γ-glutamyl transpeptidase sequence was performed and the radial phylogram of the five peptides was constructed (Fig. 11). In the radial phylogram, the root of the AsGGTs subtree is the point at which the outgroup (the Bacillus licheniformis peptide) connects to the subtree of the AsGGT peptides. The root, which represents the common ancestor of the four AsGGTs, appears to have given rise to two ancestral peptides, which subsequently gave rise to the peptides AsGGT1 and AsGGT2, and the peptides AsGGT3 and AsGGT4.

The predicted structures of all AsGGTs are relatively similar and consist of one β-sandwich which is surrounded by clusters of interacting α-helices. The two cartoon-like visualizations of AsGGT4 produced by SWISS-MODEL (Fig. 12 and AlphaFold (Fig. 13) (the color spectrum in both visualizations reflects the confidence of prediction) were superimposed (Fig. 14) with PyMOL software (52). This revealed that predicted structures by both algorithms are very similar. Most importantly, the catalytic T372 remains at the same position in the active site in both predictions. The position of a flexible loop relative to the active site appears displaced between the two structures, with low prediction confidence in both structures. The main difference between the two predictions is that an unstructured sequence which corresponds to the predicted signal peptide, only appears in the AlphaFold structure, with low prediction confidence. As with AsGGT4, AlphaFold also fails to predict any secondary structure for the predicted signal peptide of AsGGT3. Notably, AlphaFold predicts a helical structure for the DeepTMHHM-predicted transmembrane helix of AsGGT1, but not for the same feature in AsGGT2.

A total of 366 paired-end and five single-end Illumina RNA-Seq runs were identified (until March 24, 2023) in garlic. In the vast majority of these, reads corresponding to AsGGT3 transcripts were identified. On the other hand, a small proportion of runs contained AsGGT4 reads. The majority of these were identified on the SRR13219906 run, which corresponded to a leaf tissue replicate of garlic grown under normal soil moisture conditions. Other runs containing reads which corresponded to AsGGT4 transcripts included those of the PRJNA566287 and PRJNA489986 BioProjects. Their runs derive from long-vernalization, short-vernalization and non-vernalized stored clove and leaf and clove samples, respectively.