Predicted protein 3D structures provide essential insights into the genetic architecture underlying phenotypic diversity in maize

Protein structure predictions for 26 diverse maize inbred lines

To investigate the structural diversity, we used AlphaFold2 to predict the structure of proteins for the maize inbred line B73 and the 25 other NAM founder lines. Whether AlphaFold2 can predict the effects of single amino-acid polymorphisms (SAPs) is controversial (Stein and Mchaourab 2024), so we first compared the protein sequence from the same core pangene group to characterize the sequence diversity of orthologous proteins among the 26 maize NAM founder lines (Hufford et al. 2021). We identified sequence variants in ∼70% of pairwise protein comparisons, whereas only ∼10% differed by a single point mutation (Fig. 1A), indicating the high level of sequence diversity within this panel.

Figure 1.

Protein 3D structure predictions for 26 maize NAM founder lines using AlphaFold2. (A) Protein sequence diversity in maize NAM founder lines. Each NAM founder line was compared with other NAM founder lines in pairs. (SAP) Single amino-acid polymorphism, (indel) insertion/deletion. (B) Percentage of predicted protein 3D structures at different confidence levels in B73 (left; total count: 72,539) and NAM founder lines (right; total count: 1,089,242); colors represent median confidence scores. Unresolved protein structures were not predicted; see Methods. (C) Distribution and percentage of confidence scores for all NAM founder lines residues. The histogram shows the density distribution of confidence scores, and the pie chart shows the percentage of confidence scores. (pLDDT) Predicted local distance difference test. (D) Distribution of RSA for residues in all NAM founder lines. (RSA) Relative solvent-accessible surface area. (E) Heatmap representing median confidence scores of protein structures as a function of protein length in all NAM founder lines.

Open in new tabDownload PowerPointLink to figure

We predicted the 3D structures of 68,262 proteins from B73 (v5.0; including canonical proteins from 37,724 genes) using AlphaFold2 (Jumper et al. 2021), as B73 has been widely used as the reference maize line for more than a decade (Schnable et al. 2009). For the other 25 founder lines, we merged their canonical protein sequences into a nonredundant collective protein set (excluding protein sequences identical to their respective B73 reference). Overall, we predicted the structures of 313,639 unique protein sequences, corresponding to 795,649 proteins for all 26 NAM founder lines (Supplemental Fig. S1A,B). The folded protein structures represent 94.1% of all proteins from B73 and ∼73% of proteins from the NAM founder lines (Fig. 1B; Supplemental Table S1). To assess the prediction confidence for each protein and each residue, we used the predicted local distance difference test score (pLDDT) (Jumper et al. 2021). At a per-protein level, 59.2% and 49.1% of all proteins showed high confidence (median pLDDT ≥ 70) in the B73 and the NAM founder lines, respectively (Fig. 1B; Supplemental Table S1). A more stringent cut-off (median pLDDT ≥ 90) predicted 33.6% (B73) and 29.6% (NAM founder lines) of folded proteins with very high confidence (Fig. 1B). At the per-residue level, 61.0% (B73) (Supplemental Fig. S2A) and 64.4% (NAM founder lines) (Fig. 1C) of the residues were predicted with high confidence (pLDDT ≥ 70) (Supplemental Table S2). The pLDDT scores for each residue followed a clear bimodal distribution (Fig. 1C; Supplemental Fig. S2A), similar to that of humans (Alderson et al. 2023).

RSA is a measure of the exposure of an amino acid residue to the protein's solvent (Ramsey et al. 2011; Savojardo et al. 2021), with values ranging from zero for a totally buried residue to one for a totally exposed residue. RSA followed a bimodal distribution, with low pLDDT scores mainly distributed in the high RSA region and high pLDDT scores mainly distributed in the low RSA region for both B73 (Supplemental Fig. S2B) and the NAM founder lines (Fig. 1D). Furthermore, RSA was significantly negatively correlated with pLDDT (Supplemental Fig. S3). Residues with high-confidence scores are potentially functional domains and are tightly folded (Akdel et al. 2022). Structural domains are fundamental units critical to biological function, typically ranging in length from 100 to 500 residues (Wheelan et al. 2000). Our results reveal that high-confidence levels are predominantly enriched in domain-like length regions in B73 (Supplemental Fig. S2C) and the NAM founder lines (Fig. 1E). The folded maize proteins structures with high sequence diversity represent valuable resources to investigate the evolution of protein structures and conduct quantitative genetics analysis.

Buried residues are under stronger purifying selection than exposed residues

Natural selection acting on protein 3D structures plays key roles in environmental adaptation and fitness with structural constraints representing one of the main determinants of the rates of evolution (Wolf et al. 2008). Therefore, we aimed to explore the molecular mechanisms underlying natural selection via predicted 3D structures. We identified 495,357 nonsynonymous and 490,815 synonymous single-nucleotide polymorphisms (SNPs) (Bukowski et al. 2018) from the maize 282 association panel (Flint-Garcia et al. 2005). High-confidence and buried regions (with high pLDDT and low RSA) held higher proportions of synonymous SNPs than nonsynonymous SNPs, indicating that these regions are likely subject to stronger purifying selection pressure (P-value < 2.2 × 10⁻¹⁶, Fisher's exact test) (Fig. 2A,B). To further investigate how natural selection acted on protein structure, we used Z. mays ssp. mexicana and Z. mays ssp. parviglumis as outgroups to infer ancestral states and observed that nonsynonymous SNPs in the lowest derived allele frequency bin held the lowest mean RSA value compared with the other bins (P-value < 0.05, one-way analysis of variance and Fisher's least significant difference test) (Fig. 2C), suggesting that purifying selection acts on buried residues.

Figure 2.

Protein structure predictions provide functional insights into natural selection in maize. (A) Density distribution of pLDDT scores between synonymous and nonsynonymous SNPs. (B) Density distribution of RSA between synonymous and nonsynonymous SNPs. (C) Distribution of RSA for different derived allele frequencies. (D) Relationship between GERP scores and buried (RSA < 0.5) or exposed residues. Median value is shown as a solid middle line within the box plot, which spans from the 25th to 75th percentiles. For definition of buried and exposed residues, see Supplemental Figure S4. (E) Correlation between RSA and GERP scores divided by pLDDT score bins (as indicated by the gray shading above the plots).

Open in new tabDownload PowerPointLink to figure

Genomic sites with higher genomic evolutionary rate profiling (GERP) scores are interpreted to be under stronger purifying selection (Davydov et al. 2010; Yang et al. 2017). Genomic sites at buried residues (RSA < 0.5) had greater GERP scores than those at exposed residues (RSA ≥ 0.5), suggesting that buried residues are under stronger purifying selection (P-value < 2.2 × 10⁻¹⁶, t-test) (Fig. 2D; Supplemental Fig. S4). RSA is negatively correlated with pLDDT (Supplemental Fig. S3), and pLDDT decreases substantially when the median alignment depth is less than about 30 sequences in the multiple sequence alignment (MSA) from the AlphaFold2 model (Jumper et al. 2021). A high depth of MSA is likely to reflect sequence conservation. To avoid this source of confounding (Chakravarty and Porter 2022), we analyzed the strength of purifying selection conditionally on pLDDT bins. In each bin, residues with low solvent accessibility held higher GERP scores compared with residues with high solvent accessibility (Fig. 2E), suggesting that they are subjected to stronger purifying selection than those with high solvent accessibility, even after accounting for the confounding effect of pLDDT. Moreover, πN/πS (the ratio of nonsynonymous to synonymous nucleotide site diversity) was lower in buried regions than in exposed regions, even when considering the confounding effect of pLDDT scores, confirming purifying selection in buried regions (Supplemental Fig. S5). Taken together, these observations indicate that predicted protein 3D structure by AlphaFold2 at population scale could reveal associations between structural features of protein and evolutionary constraint, providing key functional insights into natural selection.

Predicted structural variants provide novel insights into genetic effects on phenotypes

To conduct association analysis based on protein variation, we measured the structural similarity as well as sequence similarity for the canonical isoforms of each core pangene group. We used the principal components (PCs) of the structure- and sequence-based similarity matrices as independent variables to perform association analysis on 32 traits. To account for population structure and relatedness, we included a genome-wide relationship matrix (based on genomic SNPs) and a genome-wide IBD matrix (based on protein haplotypes) in a linear mixed model (Supplemental Fig. S6). Our PWAS focused on PCs with eigenvalues higher than a predefined threshold (1 × 10⁻⁵), which included most PCs with nonzero variance (Supplemental Fig. S7A). Lowering this predefined threshold resulted in near-identical results from structure-based tests (Supplemental Fig. S7B–E). It is noteworthy that this lower threshold also resulted in spurious associations with extremely high significance values (e.g., pangene 9228, with −log₁₀(P) > 200) (Supplemental Table S3). On average, the associations between predicted structure and phenotype were stronger than those between sequence and phenotype (Fig. 3A; Supplemental Fig. S8), suggesting higher statistical power through structure-based associations. We identified 108 significant proteins through both structure-based and sequence-based association analyses (96 and 84 by structure- and sequence-based tests, respectively). Of these, 72 proteins were shared between the two methods, whereas structure-based association identified 24 additional significant protein–trait associations unique to predicted structures (Fig. 3B; Supplemental Fig. S9), suggesting that structure-based association analysis is a valuable and informative complement to sequence-based association analysis. The variability of PCs was almost identical with different significant levels for all significant proteins (Supplemental Fig. S10). To test whether the relative advantage of structure-based PWAS was not owing to the use of specific sequence-based MSA software, we conducted the analysis using two other MSA tools: MUSCLE and T-Coffee (Notredame et al. 2000; Edgar 2004). The protein sequence similarity matrices generated from different MSA approaches were highly correlated with each other but were different from the structural similarity matrices (Supplemental Fig. S11). The observed improvement in statistical significance from predicted protein structures was consistent across MSA software (Supplemental Figs. S12, S13). Furthermore, we investigated whether the significant proteins identified via structure-based PWAS could be identified via standard SNP-based GWAS. The results suggested structured-based PWAS identified additional loci that were not identified by SNP-based GWAS (Fig. 3C,D; Supplemental Table S4).

Figure 3.

Comparison between association analyses. (A) Decay of significance over physical distance from significance peaks in PWAS. The peak is at the starting point of a significant gene. (B) Significant association signals from structure- and sequence-based PWAS. The dashed lines correspond to a 5% significance threshold with Bonferroni correction. (C) Significant loci from structure-based PWAS and SNP-based GWAS. The dashed lines correspond to a 5% significance threshold with Bonferroni correction (5.32 in green for structure-based PWAS, 7.34 in purple for SNP-based GWAS). In SNP-based GWAS, significant genes are selected based on the most significant SNP in coding regions. (D) Manhattan plot of structure-based PWAS and SNP-based GWAS. Traits shown are ear row number, node number above ear, and upper leaf angle. The dashed lines are identical to those shown in C.

Open in new tabDownload PowerPointLink to figure

To assess our ability to identify candidate causal genes via structure-based PWAS, we inspected associated protein structures encoded by known genes with phenotypes (Supplemental Fig. S14). The predicted protein structure of the transcription factor LIGULELESS1 (encoded by Zm00001eb067740) was significantly associated with upper leaf angle. LIGULELESS1 regulates ligule and leaf angle development, and the mutant phenotype of this gene has been genetically studied (Li et al. 2017; Mantilla-Perez and Salas Fernandez 2017). The predicted protein structure of TUBTF6 (encoded by Zm00001eb188370) was significantly associated with upper leaf angle. This gene was reported to be enriched in a leaf gene regulatory network (Bertolini et al. 2025). DWARF and IRREGULAR LEAF1 (DWIL1; encoded by Zm00001eb287100) was associated with node number above the ear and was previously reported to influence plant height, internode length, and potentially node number (Jiang et al. 2012). GNARLEY (encoded by Zm00001eb117820) was associated with node number above the ear. Its mutation alters the overall plant architecture (Foster et al. 1999), and its ortholog in rice, Oryza sativa homeobox 15 (OSH15), affects node development (Sato et al. 1999). Because of linkage disequilibrium (LD), the extra power gained from structured PWAS might not point to the causal genes. The predicted structure of GOLDEN2-LIKE 37 (GLK37; encoded by Zm00001eb073790) showed a stronger association with ear row number than its sequence, and KRN2 (Zm00001eb073740) is in high LD with GLK37. KRN2 was previously reported to regulate ear row number (Chen et al. 2022). Genetic research is needed to confirm and further understand the underlying molecular mechanisms of these significant protein–trait associations.

Protein 3D structure improves genomic prediction accuracy

We assessed the additional contributions of protein structure variation to the genomic variance of the 32 traits in the NAM population by performing variance partitioning using a genomic prediction model. Because the genotypic variants used for association are not fully independent of each other, we estimated their effect sizes (i.e., the amount of phenotypic variance explained) using a mixed model and included the whole-genome SNP-based genomic value as an independent variable. We constructed three proteome-wide similarity matrices using structure, sequence, and IBD from all pangene groups. We first estimated the variance explained by the sequence similarity matrix and then added the protein structural similarity matrix as an additional source of variance (Fig. 4A). The inclusion of the structural similarity matrix provided a better fit compared with using the sequence similarity matrix alone (likelihood ratio test, P-value < 0.05) (Fig. 4A). Structure-based PWP improved prediction accuracy by 3.8% on average compared with sequence-based PWP (P-value = 1.1 × 10⁻⁷, paired t-test) (Fig. 4B), in which node number below ear, fructose content, and glucose content yielded >10% improvement in prediction accuracy. Structure-based PWP incorporating sequence information improved prediction accuracy by 4.1% compared with sequence-based PWP on average (P-value = 2.7 × 10⁻⁷, paired t-test) (Fig. 4C), in which fumarate, node number below ear, protein, fructose content, and glucose content yielded >10% improvement in prediction accuracy. The increased prediction accuracy was independent of MSA methods used (Supplemental Figs. S15, S16). Structure-based similarities, alone or combined with IBD-based similarities, resulted in increased prediction accuracy compared with IBD-based similarity alone (1.9% and 2.5%, respectively) (Supplemental Fig. S17), whereas node number below ear led to the highest improvement in prediction accuracy (12.0% and 11.4%). These results indicate that protein structures explain genotypic variation not captured by protein sequence or IBD information, which translates into higher genomic prediction accuracy.

Figure 4.

Protein 3D structure improving prediction accuracy compared to protein sequence. (A) Variance partitioning in different genomic models for 32 traits. Sequences were aligned using MAFFT. Likelihood values are listed in Supplemental Table S5. (B) Comparison of prediction accuracy between sequence-based PWP and structure-based PWP. (C) Comparison of prediction accuracy between sequence-based PWP and sequence-based plus structure-based PWP. Dots represent the mean values of 25 leave-one-family-out repeats, and error bars indicate their standard errors.

Open in new tabDownload PowerPointLink to figure

Prediction of protein structure requires high-quality genome sequence and annotation

Prediction of protein structure heavily depends on high-quality protein sequences, as the 3D conformation of a protein is determined by the arrangement of its amino acids. Sequencing errors or gaps in the sequence can lead to inaccurate or misleading predicted structures. Therefore, ensuring the completeness and accuracy of the genomic sequence data is crucial for reliable structure prediction. Highly accurate gene annotation (the process of finding and designating the locations of individual genes) is also required to improve structure predictions (Salzberg 2019). With the advent of Oxford Nanopore Technology (ONT) and Pacific Biosciences (PacBio) HiFi sequencing technologies, more and more genomes are being assembled at the telomere-to-telomere (T2T) level (Xie et al. 2024). As a result, the acquisition of high-quality protein sequences has become increasingly feasible, enabling large-scale protein structure predictions and opening new avenues for the exploration of functional mechanisms.

Genetic effects are captured by high-confidence protein 3D structure

Proteome-wide protein 3D structure prediction accurately captures structural differences between genes in the reference genome of human and other model organisms (Tunyasuvunakool et al. 2021; Akdel et al. 2022). As protein sequence alleles are conserved across mammalian species, AlphaFold2 might not have the power to predict the impact of SAPs (Stein and Mchaourab 2024). However, AlphaFold2 may effectively predict allelic differences in protein structures in plant species, which are significantly more diverse than in mammal species (Buckler et al. 2006; Song et al. 2024). In this study, we extend proteome-wide protein 3D structure prediction to a population scale by taking advantage of the design of the maize NAM population. This allowed us to use predicted structure to analyze phenotypic variation in quantitative genetics studies, thereby enabling the decomposition of genotypic variability, which would otherwise not be possible.

Predicted protein structural similarity is more informative than protein sequence similarity

Structural similarity between proteins is a strong predictor of functional similarity (Erdin et al. 2011). Vastly different amino acid sequences can lead to similar structures, although even two closely related sequences may fold into distinctly different structures (Krissinel 2007). Our structure-based PWAS results suggest that 3D structures predicted from sequences are more useful than sequences alone for discovering candidate genes. Our method filtered some PCs with low variance, which introduced extra false positives in PWAS. We examined several candidate genes that contribute to phenotypic variation. Although we identified promising candidates, sequence differences in proteins do not necessarily affect their predicted 3D structure. For omic-wide associations, we applied the canonical Bonferroni multiple tests correction to adjust significance thresholds for SNPs and pangene groups separately. Owing to the substantially larger number of SNPs compared with pangene groups, the resulting significance thresholds (expressed as −log₁₀(P-value)) were higher for SNPs. Nevertheless, even applying the same threshold to SNP-based GWAS and PWAS, the observation that structure-based PWAS provides complementary information to SNP-based GWAS still holds (Fig. 3C). In our analysis, we only kept pangenes that were present in all individuals, although additional genes are dispensable or private to some founder lines. Our analysis focused on sequence variability and did not consider potential effects of gene absence. Moreover, the dosage of functional molecules in cells, which is also essential for trait regulation (Wingo et al. 2021), could not be captured by predicting protein structure. Including dosage regulation and protein complexes in the model might further improve our understanding of the genetic mechanisms underlying trait diversity. Predicted 3D structure may only describe one component of gene activity, independent from gene regulation. Therefore, our proposed structure-based PWAS approach should ultimately be combined with other approaches (e.g., proteome-wide dosage association studies) to obtain a detailed biological description of genotypic variability.

Predicted protein structural variants capture genotypic variability

There have been mixed results about the usefulness of AlphaFold2 for predicting the effect of single-residue variants on protein structure (Buel and Walters 2022; McBride et al. 2023). Our PWAS and PWP results clearly indicate that AlphaFold2 can capture genotypic variability for agronomic, morphological, and compositional traits in maize. Importantly, predicted protein structural variants in our study consisted of haplotype variants, not necessarily single-residue variants. Based on our results, we hypothesize that AlphaFold2 can accurately predict allelic variants for protein structure provided that there are large enough differences among protein alleles. This hypothesis calls for further research to quantify the structural differences between protein alleles and to test the accuracy of their prediction.

In our analyses, precision was prioritized over power, as our PWAS approach was essential to account for biases introduced by LD, by including both SNP-based relationship and haplotype-based IBD matrices (Supplemental Fig. S6). A previous study predicted mRNA abundance for each haplotype in the NAM population, similar to our prediction of protein structures (Giri et al. 2021). Importantly, in this study, permutations of predicted gene activity across the 26 haplotypes in the NAM population were predictive of phenotypes, as long as haplotypes were still correctly assigned to individual NAM lines. This result highlights the potential source of confounding owing to haplotype variation alone, as well as the necessity to account for haplotype IBD in association studies of gene activity imputed by haplotype.

Although our PWP results suggest that predicted 3D structure provides useful functional information about genetic effects, the gains in prediction accuracy realized in our study are certainly too low to justify routine applications of AlphaFold2 in genomic selection programs. However, the PHG used in this study may still be used for imputing haplotypes in any breeding population, and our structure-based similarity matrices may then be used without running AlphaFold2 again. Further research is needed to determine whether haplotype imputation by this PHG database (or future databases containing more genome assemblies) is accurate enough for effective incorporation of predicted protein structure information into genomic prediction models.

Protein sequence diversity investigation

To investigate sequence diversity, we extracted the pangenes from previously published pangene annotations (Hufford et al. 2021) with at least one predicted protein structure for each of the 26 founder lines (see below), resulting in 17,633 pangenes. MSA was performed using MAFFT to investigate the sequence diversity (Katoh and Standley 2013) for each pangene group. Protein sequences from the 26 NAM founder lines were iteratively used as the reference. The query sequences were compared with the reference sequence for each pangene of each iteration to count the number of SAPs and indels.

Protein 3D structure prediction

Protein structure prediction was performed using AlphaFold v2.0.0 with the preset reduced_dbs parameter as described in the AlphaFold database publication (Tunyasuvunakool et al. 2021). The genome annotations and pangene results were generated by Hufford et al. (2021). All B73 v5.0 protein sequences and canonical protein sequences from the other NAM founder lines were selected as input. Sequences with residue codes “*” (premature stop codon) or “X” (unknown residue) were also excluded. The 3D structures of all B73 protein sequences were predicted, and canonical protein sequences from the other NAM founder lines were extracted for 3D structure prediction. All the kept protein sequences were merged to check for duplicated protein sequences. By default, AlphaFold v2.0.0 outputs five structural models for each input sequence; the model with the highest pLDDT was selected for subsequent analyses. Because of the high computational cost, we only folded proteins with fewer than 1200 amino acids for the 25 NAM founder lines, except for B73. Because of high GPU memory requirements, some proteins failed to fold. Protein structures were defined as resolved if they were successfully predicted by AlphaFold2 and were otherwise considered to be unresolved. Protein folding was conducted using high-performance computing (HPC) computational nodes with 4 NVIDIA Quadro RTX 5000 GPU cards.

Calculation of RSA

The accessible surface area (ASA) for each amino acid residue in a protein structure was calculated using DSSP v2.3.0 (Touw et al. 2015). The maximum possible solvent accessible surface area (MaxASA) was obtained from a previous publication (Tien et al. 2013). The RSA of each amino acid residue was calculated using the formula RSA = ASA/MaxASA.

SNP annotation

The variant calling file for 282 association panels was downloaded from maize HapMap 3.2.1 (Bukowski et al. 2018) and uplifted to B73 v5.0 reference coordinates via CrossMap (Zhao et al. 2014). Allele frequencies were counted using VCFtools (Danecek et al. 2011). Based on genome annotation information in general feature format (GFF3), SNPs in protein-coding regions were extracted and classified as either synonymous SNPs or nonsynonymous SNPs. The SNP data set was annotated by GERP score from a previous publication (Kistler et al. 2018) uplifted to B73 v5.0 coordinates.

Calculation of nucleotide diversity (π)

The biallelic SNP file for each line in the 282 association panel was converted to a genomic sequence file using B73 v5.0 reference coordinates as the reference. Protein-coding sequences were extracted using GffRead (Pertea and Pertea 2020), and the same transcripts from all lines were merged into MSA files. Each codon was classified based on pLDDT bin and RSA bins. RSA value was divided into 10 bins ranging from zero to one at intervals of 0.1. πN and πS (π of nonsynonymous and synonymous mutations, respectively) were calculated by bppsuite (Guéguen et al. 2013), with MSA files used as input. The average πN and πS values for all the codons in each bin were used to perform statistical analysis and generate the plot (Supplemental Fig. S5).

Estimation of derived allele frequency

To infer the derived allele for each SNP in the maize population, the genomes of Z. mays ssp. mexicana (TIL18) and Z. mays ssp. parviglumis (TIL11) were used as outgroups. The sequences of two genomes are available at MaizeGDB (https://download.maizegdb.org/). The genomes of the two outgroups were aligned against the B73 v5.0 reference genome (Hufford et al. 2021) using AnchorWave (Song et al. 2022). The genome alignments were converted into a GVCF format file using the MAFToGVCF plugin of TASSEL (Bradbury et al. 2007). For each SNP in the maize population, if the reference allele matched either the TIL18 or TIL11 allele, the reference allele was categorized as an ancestral allele. Otherwise, if the alternative allele was identical to the TIL18 or TIL11 allele, the reference allele was defined as a derived allele. For a SNP, whose reference position was not aligned to either of those two outgroups and was therefore likely located in newly derived sequence fragments, the SNP was excluded from the analysis.

Variant calling for the NAM population and kinship matrix construction

SNP variant calling of the NAM founder lines was conducted using PHG v1 (Bradbury et al. 2022). A PHG database was constructed from the genomes of 26 NAM founder lines (Hufford et al. 2021) using the B73 v5.0 assembly as the reference. GBS reads from 4736 accessions were mapped to the pangenome, and SNP variants were imputed in the NAM RIL population. The PLINK (Purcell et al. 2007) program was used to retain only biallelic variants with the parameter “‐‐max-alleles 2,” resulting in 42,329,376 SNPs for 4736 inbred lines. The SNP-based kinship matrix G was computed using the Balding–Nichols (BN) methods in the EMMAX software (Kang et al. 2010).

Protein similarity matrix construction

The NAM population was used for PWAS and PWP. To predict protein structures at population scale, all metrics for the NAM founder lines were projected onto the NAM population with 25 families using the PHG (Supplemental Fig. S18; Bradbury et al. 2022), assuming no recombination crossover in the middle of coding genes. To measure the structural similarity within each core pangene group j, structure alignments were performed using US-align, and the structural similarity matrices were computed as TM scores with values in [0,1], using structure alignment results as input (Zhang et al. 2022). To construct sequence similarity matrices, protein sequences were aligned via one of three MSA approaches, that is, MAFFT (Katoh and Standley 2013), MUSCLE (Edgar 2004), or T-Coffee (Notredame et al. 2000). The similarity scores for each pair of sequences were calculated as numerical values in [0,1] as well. Each pair of amino acids was treated as similar to each other if they had a positive score from the BLOSUM62 matrix; otherwise, they were considered to be dissimilar to each other. For each pair of protein sequences, the similarity score was calculated as the number of similar amino acids divided by the average length of these two protein sequences. The IBD matrix for each pangene group was constructed using an identity matrix to NAM founder haplotypes.

For each type of similarity (structure-, sequence-, or IBD-based) and pangene group j, K_j was the 26 × 26 relationship matrix among the 26 founder haplotypes computed from similarity matrix S_j (https://github.com/shuaiwang2/Protein3D_QG/blob/main/PWAS/PWAS-data_processing.R). First, each relationship matrix was obtained by centering the corresponding similarity matrix: K_j = HS_jH′, where H is the centering matrix H = I₂₆ − J₂₆/26. The pangene groups were then filtered according to the following criteria: (1) no missing haplotypes among the 26 NAM founders, (2) successful protein structure prediction for all observed haplotypes, and (3) trace of similarity matrix (sum of diagonal elements after centering) greater than 1/25. This filtering step resulted in m = 11,927 retained pangene groups across similarity metrics (US-align, MAFFT, MUSCLE, T-Coffee, IBD). The distribution of pangene is almost even on the chromosomes, except for a missing portion of Chromosome 10 (Supplemental Fig. S19).

The input to PWAS consisted of one or more PCs per retained pangene group, obtained by eigen decomposition of K_j for each pangene group j. To retain eigenvalues that contribute substantially to the overall variance of the matrix, for each pangene group j, PCs were included in PWAS if their variance explained in K_j was greater than 1 × 10⁻⁵. If there were no PCs kept, the pangene group would not be used for the following analyses. The average number of structured-based PCs included (k) was seven, with a range from one to 25. The average number of sequence-based PCs included was six, with a range from one to 25. For each pangene group j and relationship matrix K_j, the matrix of PCs used in PWAS was ${\mathbf{X}}_j={\mathbf{Z}}_j{\mathbf{U}}_j{\mathbf{\Lambda }}_j^{1/2}$, where Z_j is the n × 26 design matrix assigning each individual to one of the 26 founder haplotypes at the pangene group j, U_j is the 26 × k matrix of retained eigenvectors of K_j, and Λ_j is the k × k diagonal matrix of the corresponding eigenvalues.

The input to PWP consisted of the sum of sequence-based, IBD-based, or structure-based similarity matrices across all pangene groups. To ensure mathematical validity, similarity matrices were adjusted to positive definite matrices using the nearPD function from the Matrix package in R (Higham 2002; R Core Team 2025). To normalize proteome-wide relationship matrices, the mean of their diagonal elements was set to one as follows:

PWAS

PWAS was performed using the following linear mixed model for each pangene group j:

A total of 11,927 pangenes were investigated in PWAS separately. A total of 32 traits were analyzed in PWAS. The phenotypes used were selected from previous publications (Buckler et al. 2009; Brown et al. 2011; Kump et al. 2011; Poland et al. 2011; Tian et al. 2011; Cook et al. 2012; Hung et al. 2012; Peiffer et al. 2013; Olukolu et al. 2014; Wallace et al. 2014; Benson et al. 2015; Foerster et al. 2015; Bian and Holland 2017; Diepenbrock et al. 2017), of which the most highly correlated traits were excluded (Supplemental Table S6).

SNP-based GWAS

GWAS was implemented by GEMMA (Zhou and Stephens 2012) using default parameters. A total of 19,458,794 SNPs were used for SNP-based GWAS. SNPs with a minor allele frequency greater than 0.05 and missing rates <20% were retained using PLINK (Purcell et al. 2007). We calculated the GWAS significant threshold using Bonferroni multiple-testing correction based on the number of effective SNPs. The effective SNP count was derived with an R² threshold of 0.99, a window size of 100,000, and a step size of 50,000, resulting into 1,090,989 SNPs. The GWAS used the same phenotypic traits as those used in the PWAS.

Variance partitioning

To assess the significance of structure-based PWP, variance partitioning was performed using the MMEst function implemented by the MM4LMM package (Laporte et al. 2022). Proteome-wide relationship matrices were used as input for variance partitioning. The “baseline” prediction model only included a SNP-based kinship matrix G. The “sequence” model included a G matrix and a proteome-wide relationship matrix P based on protein sequences. The “IBD” model included a G matrix and a P matrix based on haplotype IBD. The “structure” model included a G matrix and a P matrix based on predicted structures. The “sequence+structure” model included a G matrix and two P matrices based on sequences and predicted structures. The “IBD+structure” model included a G matrix and two P matrices based on haplotype IBD and predicted structures.

To test the statistical significance of additional random effects (e.g., structure-based genomic values), a likelihood ratio test was conducted to compare different genomic models using the following formula:

Proteome-wide prediction

To estimate the improvement of prediction accuracy using structure information, PWP was carried out using the greml function implemented in the qgg package (Rohde et al. 2020). The model applied was as follows:

In models with only one proteome-wide relationship matrix, p is a vector of “proteomic” values, such that $\mathbf{p}\sim \mathcal{N}(\mathbf{0},\mathbf{P}{\sigma }_p^2)$, where P is based on either structure, sequence, or IBD. In models with two proteome-wide relationship matrices, p = p₁ + p₂, where ${\mathbf{p}}_1\sim \mathcal{N}(0,{\mathbf{P}}_1{\sigma }_{p_1}^2)$, ${\mathbf{p}}_2\sim \mathcal{N}(\mathbf{0},{\mathbf{P}}_2{\sigma }_{p_2}^2)$, with P₁ and P₂ being two different proteome-wide relationship matrices.

Considering the population properties of the NAM population, PWP models were validated by leave-one-family-out cross-validation. Each NAM family was left out as a validation data set to compute the prediction accuracy separately, whereas the other 24 NAM families were used for training. Prediction accuracy was computed as $Cor(\mathbf{y},\hat{\mathbf{y}})$, where $\hat{\mathbf{y}}$ is the vector of predicted genotypic values: $\hat{\mathbf{y}}=\hat{\mathbf{g}}+\hat{\mathbf{p}}$. This process was repeated for 25 iterations.

Statistical analysis

All statistical analyses used for this paper were performed in R (R Core Team 2025) and are indicated in the main text and figure legends.