A total of eight synthetic wheats, produced from different cross combinations between tetraploid wheats and Aegilops species containing D-genome [ Table 1 ; (Mirzaghaderi et al., 2020)] were chosen. For each genotype, 12 seeds were planted in three pots containing loam soil (four seeds per pot). To ensure vernalization and more even flowering, the seedlings were first placed in a growth chamber at 4°C for 1 month, and then kept in a greenhouse with photoperiod of 16 hours light: 8 hours dark at 22°C. All pots were watered once a week and kept under the same conditions until the initiation of flowering and sampling.

Seeds were sampled during grain filling at 21 days post-anthesis when the kernels were at milk to soft dough stage. Seeds were immediately frozen in liquid nitrogen, and stored at -80°C. RNA isolation was performed using the single step extraction method by acid guanidinium thiocyanate–phenol–chloroform developed by (Chomczynski and Sacchi, 2006). Briefly, RNA was extracted from at least five immature seeds of each single plant by adding 1 ml of denaturing solution (4 M guanidinium thiocyanate, 25 mM sodium citrate, pH 7.0, 0.5% (wt/vol) N-laurosylsarcosine (Sarkosyl) and 0.1 M 2-mercaptoethanol) into the crushed grains in liquid nitrogen, and followed by phase separation in an acidic pH to ensure the separation of RNA from DNA and proteins. After precipitation of extracted RNA in isopropanol, RNA was washed using 75% ethanol. The RNA pellet was air-dried and then dissolved in 50 µl of DEPC-treated water and stored at -80°C until cDNA synthesis. RNA extraction using the protocol mentioned above at acidic pH usually results in pure RNA. To test whether our extracted RNA samples are free of DNA, we performed an additional PCR on each RNA sample using alpha-gliadin primers. No band was detected after running the PCR products on agarose gel ( Supplementary Figure S1A ), so no DNase treatment was applied. The quality and quantity of extracted RNA were measured using agarose gel electrophoresis and NanoDrop spectrophotometry. For all samples, the concentration of RNA was equalized. Three biological RNA samples were then pooled for each genotype, and 1 µg of the mixed RNA was used for cDNA synthesis. cDNA was prepared using the TaKaRa kit (Otsu, Shiga, Japan) according to the vendor’s protocol.

Alpha-gliadin sequences were amplified from cDNA samples in two steps. The first amplification was performed using gene specific primers of alpha-gliadin; AlphaF (5′-atgaaracmtttcycatc-3′; for the MKTF[LP]I- motif) and AlphaR (5′- ctgctgctgtgaaattrgwt-3′; for the PISQQQ-motif) ( Figure 1 ; Supplementary Figure S1C ). Both forward and reverse primers were chosen from (Salentijn et al., 2013), which specifically amplify the repetitive domain of a typical alpha gliadin gene coding the major CD epitopes including the p31-p43 peptide and the toxic 33-mer fragment containing DQ2.5-Glia-α1, DQ2.5-Glia-α2, and DQ2.5-Glia-α3 peptides ( Figure 1 , underlined in blue). To lower the amplification bias, amplification of alpha-gliadin amplicon was performed on two different PCR thermocyclers. PCR reaction was prepared in 10 µl including 25 ng template cDNA, 5 pmol of each primer, 2.5 mM of each dNTP, 2.5 mM MgCl₂ and 0.5 U Taq polymerase. The PCR cycling parameters included an initial step at 94°C for 5 min, followed by 30 cycles of (30 seconds at 94°C, 15 seconds at 50°C, and 30 seconds at 72°C), and a final extension at 72°C for 5 min. Then, for each sample, the PCR products (two technical replicates) were mixed thoroughly and diluted 1 to 10 using nuclease-free water. 1 µl of diluted PCR product from the first round of PCR was used as a template for the second round of amplification. The second round of PCR was performed in 10 µl reactions using fusion primers. Fusion primers included the sequences needed for Illumina sequencing, and a gene specific part. The forward fusion primer also contained an 8 bp ID sequence (barcode) specific for each genotype ( Supplementary Figure S1C ) to enable multiplexing. The second-round PCR condition was similar to the first PCR amplification. The PCR products from eight different genotypes were pooled, and 80 µl of the pooled PCR product was sent for paired-end Illumina amplicon sequencing (Novogene (UK) Company). The resulting reads belonging to each genotype could later be identified via the ID oligo sequence of the forward primer.

An overall pipeline for pre-processing and analyzing of alpha-gliadin amplicon reads is shown in Figure 2 . The paired-end Illumina sequencing of the pooled alpha-gliadin amplicons resulted in a total of 6,541,604 paired-end reads (250 bp length per read). The alpha-gliadin reads were first checked for their quality using FastQC (Andrews, 2010). Then, trimming and base quality control of the reads were performed using trimmomatic tool (Bolger et al., 2014). The forward and reverse reads were concatenated/assembled using FLASH tool (Magoč and Salzberg, 2011) in Linux. The concatenated sequences belonging to each genotype were separated using the genotype-specific oligonucleotide barcodes embedded in the forward primers ( Supplementary Figure S1C ). Next, for each genotype, the sequences were further trimmed to remove the genotype identifier barcode, signal peptide, and the forward and reverse primers ( Figure 1 ). These preprocessed codding sequences (CDSs) were translated to proteins using TBtools (Chen et al., 2023). Protein sequences were then clustered using CD-HIT at 100% identity (Li and Godzik, 2006) in Linux, and the representative sequences with more than 20 members were chosen for alpha-gliadin immunogenic peptides analysis. The frequency of CD epitopes was calculated based on the abundance of the amplicons containing them. Briefly, different alpha-gliadin epitopes, including the canonical epitopes and variants with one or two mismatches were searched only in the representative clusters with more than 20 members. For each single representative protein, the number of identified epitopes was multiplied by the number of its cluster members to determine the abundance of epitopes across the clusters. To estimate the abundance of a specific epitope across the lines, the frequency of that epitope was summed for all clusters. Subsequently, the number of each calculated epitope was normalized based on read depth per line, and represented as relative frequency/expression rate and their log 2 values ( Supplementary Table 1 ) were visualized by a heatmap in R software using the pheatmap package (Kolde, 2019). We also coded the alpha-gliadin epitopes based on their previously reported toxicity information into three groups (T: toxic, RT: decreased toxicity, *T: highly likely toxic) and assigned scores (3 for T, 2 for RT, and 1 for highly likely toxic) to estimate the overall toxicity load of each line ( Figure 3 ; Supplementary Table 2 ). To achieve this, the number of peptides for each group was multiplied by its corresponding score. Then, the line with the highest toxicity score was assumed 100% toxic, and the toxicity percentages of other lines were calculated relative to this line ( Supplementary Table 2 ).

To construct the phylogenetic tree, first, the ten most prevalent clusters of alpha-gliadin proteins were chosen in each genotype ( Supplementary Table 3 ). Then, the consensus protein sequences were created for each line separately using SeaView tool version 4 (Gouy et al., 2010). Logos for the alignment of the consensus proteins were generated using online WebLogo application (Crooks et al., 2004). A phylogenetic tree was also constructed using the consensus sequences in MEGA11 software (Tamura et al., 2021) with neighbor-joining method and bootstrapping of 1000 replicates.

To analyze the gliadin banding patterns of synthetic lines, we performed Acid-PAGE. First, gliadin proteins were extracted using the method described by Khan et al. (1985) with slight modifications. 200 µl of 70% ethanol was added to 0.1 g of fresh milled flour in a 1.5 ml tube. The samples were then vortexed every 10 min for 1 hour at room temperature and centrifuged at 13,000 rpm for 15 minutes. For each genotype, 20 µl of supernatant was transferred to a new 1.5 ml tube and heated at 70°C in a water-bath to form a protein pellet. Then, 20 µl of 5% acetic acid and 20 µl of sample solution (4.5 M urea solution prepared in 5% acetic acid) were added to the pellets and vortexed for 10 seconds. The prepared protein samples were stored at -20°C until use. Acid-PAGE was performed according to Watry et al. (2020). Proteins were heated at 70°C for 10 minutes before running the gel. In each well, 6 µl of proteins were loaded onto a 10% acrylamide gel and run for 90 min at 400 V. The gel image was captured using a smartphone camera of 108 MP (Xiaomi Poco X4).

Eight different synthetic wheat lines resulted from different cross combinations ( Table 1 ), were examined for their alpha-gliadin protein composition and CD epitope content. For this, RNA was extracted from immature seeds and converted into cDNA. The first variable repetitive domain of alpha-gliadins ( Figure 1 , underlined in blue) harboring the three major CD epitopes DQ2.5-glia-α1, DQ2.5-glia-α2 and DQ2.5-glia-α3 ( Figure 1 , highlighted in magenta, cyan and red, respectively), was amplified from each genotype, and the pooled PCR products were sequenced using Illumina platform. Because we used a specific barcode sequence in the used primers for amplification of each genotype ( Supplementary Figure S1C ), it was possible to pool the PCR products before amplicon sequencing. In this way, the sequences specific to each line could later be identified and separated during data processing. A summary of data pre-processing and analyzing the alpha-gliadin epitope composition across the lines is presented in Figure 2 .

Overall, 6,541,604 paired reads (an average of 600,000-800,000 reads per plant, 250 bp in size) were generated from the alpha-gliadin amplicons across all samples. After trimming, base quality control and concatenation/assembling of forward and reveres reads, collectively 2,005,705 CDSs (an average of 200,000-300,000 sequences per plant, 192-432 bp in size) were generated, from which 111,336 sequences contained stop codon(s) and were discarded from the analysis. The remaining sequences (1,894,369 CDSs) were used for analyzing and comparing the alpha gliadin CD epitopes among the genotypes. CDSs belonging to each line were separated using the specific barcode sequence ( Supplementary Figure S1C ). The sequences were then translated to proteins in TBtools, and the protein sequences of each line were clustered using CD-HIT tool with a threshold of 100% similarity. As a result, a total of 744,981 clusters were produced across all plants, which were later minimized by selecting only those clusters with at least 20 members. In this way, 700-1000 clusters were remained for each plant ( Supplementary Table 4 ).

For analyzing the expression pattern of alpha-gliadin epitopes among different synthetic wheat lines, first, a list of all previously predicted/identified alpha-gliadin epitopes, the canonical forms, all those variants with one or two mismatches and their toxicity information (if available) were retrieved from the literature (Shan et al., 2002; Ozuna et al., 2015; Wang et al., 2017; Altenbach et al., 2020; Sollid et al., 2020; Schaart et al., 2021; Marín-Sanz et al., 2023). These epitopes (121 in total) were searched in the alpha-gliadin protein clusters in each line. Overall, 54 epitopes were found in the analyzed lines ( Table 2 ): Eleven epitopes were DQ2.5-Glia-α1 and its variants. Eighteen were DQ2.5-Glia-α2 and its variants, among which variants 15, 16 and 17 containing only one mismatch, were reported for the first time. Sixteen were DQ2.5-Glia-α3 and its variants, of which, a new variant with one mismatch (variant 14) was identified. Hence, four new variants were identified which have been highlighted in Table 2 . Regarding the P31-43 peptide, beside the canonical epitope, six variants were detected. All the canonical epitopes (bolded in Table 2 ) were present across all the synthetic lines, though with different frequencies. However, the distribution and frequency of the variants varied significantly among the lines ( Table 2 ).

Specifically, the major and most immunogenic CD epitopes DQ2.5-Glia-α1a (PFPQPQLPY), DQ2.5-Glia-α1b (PYPQPQLPY), DQ2.5-Glia-α2 (PQPQLPYPQ), and DQ2.5-Glia-α3 (FRPQQPYPQ) which are located in the 33-mer region of Ae. tauschii alpha-gliadin (GenBank: ADD17012.1) ( Figure 1 ), were expressed across all synthetic wheat lines, albeit with different frequencies ( Figure 3 , Table 2 ; Supplementary Table 1 ). Collectively, the expression profile of DQ2.5-Glia-α1a and DQ2.5-Glia-α3 were higher than that of DQ2.5-Glia-α2 and DQ2.5-Glia-α1b. Among these canonical epitopes, the expression pattern of DQ2.5-Glia-α1b was the lowest especially in Lines 4, 7 and 8. The expression of each alpha-gliadin epitopes also differed across different synthetic lines. DQ2.5-Glia-α1 had the highest expression frequency in line 5 and 6, while the lowest frequency was observed in Ae. crassa derived lines 3, 4, and 7. DQ2.5-Glia-α2 was expressed with higher intensity in line 6 and with lower intensity in line 7 and 8. DQ2.5-Glia-α3 differed from other alpha-gliadins in its expression pattern among the different lines. All lines significantly expressed DQ2.5-Glia-α3 epitope with high intensity, except for lines 3 and 7 ( Figure 3 ; Supplementary Table 1 ). Our data also showed that a variant of alpha-gliadin epitope, DQ2.5-Glia-α3 (var1, FPPQQPYPQ) showed the highest expression rate across all lines. The newly identified variants of DQ2.5-Glia-α2 (var15, var16, and var17) were rarely expressed, and were only found in some synthetic lines. In contrast, the new variant of DQ2.5-Glia-α2 (var14) was found in all lines, showing higher expression levels in lines 4 and 7 ( Figure 3 ; Supplementary Table 1 ). The p31–43 peptide, which induces the innate immune response in patients susceptible to CD, and its variants (var1 and var4) had higher expression compared to other variants. It had a higher frequency in line 8 and a lower frequency in line 3. In addition to analyzing each of these CD epitopes separately, we also searched for the intact 33-mer region across the synthetic lines. The intact form of the 33-mer region (LQLQPFPQPQLPYPQPQLPYPQPQLPYPQPQPF; AFX69628.1) was found only in three lines (Lines 2, 3, and 5), showing the highest expression rate in line 2 and a lower expression rate in line 5 ( Table 1 ; Supplementary Table 5 ).

The wheat lines were also compared for the frequency of the toxic alpha-gliadin epitopes and those with the lower and higher toxicity were identified. For this, the normalized frequency of toxic (T), decreased toxicity (RT) and highly likely to be toxic (*T) epitopes were multiplied by 3, 2 and 1 scores relative to their toxic intensity, respectively. Then, the toxicity of the lines was expressed as a percentage relative to the line with the highest toxicity score as discussed in the methods section ( Figure 3 ; Supplementary Table 2 ). This generally showed that the major and toxic alpha-gliadin peptides were present in all lines, however, different synthetic lines exhibited dramatic variability in their toxicity ( Supplementary Table 2 ). Based on the frequency of toxic alpha-gliadins across the lines and our assumption (considering the line with the highest frequency of toxic CD epitopes as 100% toxic), the highest toxicity percentages were found for line 6 (100% with 435,534 T, RT, and *T peptides), followed by line 5 (~80% with 350,618 T, RT, and *T peptides). These two lines have received the D-genome from Ae. tauschii ‘299’, but their maternal parents belong to different emmer wheats ( Figure 3 , Table 1 ). Lines 1 and 2 were not significantly different for their CD toxicity (67% vs 69%, respectively). Both lines have different maternal parents, but they have been crossed with Ae. tauschii, and experienced the same backcross event with the hexaploid wheat T. aestivum ‘Pishgam’. Line 8 with the D-genome from Ae. ventricosa possessed a toxicity percentage of 49% and was the fourth line regarding CD toxicity. Lines 3, line 4 and 7, with the toxicity percentages of 37%, 42% and 32%, respectively, had the lowest contents of toxic peptides ( Supplementary Table 2 ). Interestingly, these three lines received their D genome from Ae. crassa and their AB subgenomes from the same tetraploid genotype of T. durum, except for line 3 ( Figure 3 , Table 1 ).

We also found that the frequency of epitopes with reduced toxicity varied among different synthetic lines, and some epitopes were not detected in certain lines ( Table 2 ; Supplementary Table 1 ). Epitopes DQ2.5-Glia-α1 (var3) PFSQPQLPY, DQ2.5-Glia-α2 (var8) PQPQLPYLQ, and DQ2.5-Glia-α3 (var7) FRPQQSYPQ were found in all lines. Peptide PFSQPQLPY with the higher frequencies (each containing 2,850 and 2,785 peptides, correspondingly) was found in the synthetic lines 2 and 1. While the later peptides, PQPQLPYLQ and FRPQQSYPQ had a higher proportion in the synthetic line 6. Peptide FRQQQPYPQ was a variant of DQ2.5-Glia-α3 (var11) and was only present in line 8. Other peptides such as FRPQKPYPQ, FPPQQPYTQ, PFPQLQQPY and FQPQQPYPQ), which were mainly variants of DQ2.5-Glia-alpha3, expressed only in synthetic line 8 ( Table 2 ; Supplementary Table 1 ).

The phylogenetic relationship of the synthetic lines was constructed based on the consensus sequences of the most common clusters for each line. The alignment graph of these consensus protein sequences indicated that the peptides p31-43 (PGQQQPFPPQQPY) and one of the major CD epitopes DQ2.5-Glia-α3 (FRPQQPYPQ) were conserved across all genotypes ( Figure 4A ). Based on the logo of the consensus proteins, the amino acid residues were almost equally distributed among all synthetic lines ( Figure 4B ). The phylogenetic tree of the synthetic lines also rather matched the toxicity scores of the lines. In this way, the most toxic lines (line 5 and 6) as we previously identified based on their frequency of toxic alpha-gliadin epitopes ( Figure 3 , Line toxicity), both were grouped together, and had closer phylogenetic relationship ( Figure 4C ).

We also assessed the gliadin proteins of our lines using Acid-PAGE and compared their banding patterns with those of Chinese spring wheat (T. aestivum). Though the banding patterns of ω-, β-, and γ-gliadins were not our concern in this study, they showed high variations across the genotypes ( Supplementary Figure S2 ). The banding patterns of alpha-gliadins also differed among the lines and had one or two bands missing compared to Chinese spring wheat ( Supplementary Figure S2 , red arrows). Collectively, the lines with higher abundance of toxic CD epitopes (lines 5 and 6) and the line with relatively decreased toxicity (line 8) had only one missing band, while the lines with a lower toxicity (lines 3, 4 and 7) had two bands missing. These findings generally agreed with the results of sequencing data, indicating that lines with higher amounts of toxic CD epitopes had more alpha-gliadin bands compared to those with lower toxicity.