CLY23 Phenotype Analysis: PCA and Heritability

Load Libraries

Purpose: Load required packages for data manipulation, visualization, and multivariate analysis.

Approach: Load tidyverse ecosystem and PCA visualization packages.

Expected outcome: All dependencies available for downstream analysis.

library(dplyr)
library(ggplot2)
library(lubridate)
library(MASS)
library(FactoMineR)
library(factoextra)
library(ggbeeswarm)
library(ggpubr)

Import LI-COR Data

Purpose: Import photosynthesis measurements from local CSV.

Approach: Read PhiPS2 data from exported CSV, parse sample IDs to extract field coordinates.

Expected outcome: Clean dataframe with CLY23_D4 plot IDs and PhiPS2 values.

# Read from local CSV (exported from Google Sheets)
licor_raw <- read.csv("licor_data.csv", skip = 1)

# Parse sample IDs
licor <- stringr::str_split(licor_raw$BZea, "_", simplify = TRUE) %>% 
  as.data.frame()
colnames(licor) <- c("field_year", "section", "CLY23_D4")
licor$CLY23_D4 <- as.numeric(licor$CLY23_D4)
licor$PhiPS2 <- as.numeric(licor_raw$PhiPS2)

# Diagnostics
{
cat("LI-COR records:", nrow(licor), "\n")
cat("PhiPS2 range:", range(licor$PhiPS2, na.rm = TRUE), "\n")
}

## LI-COR records: 2561 
## PhiPS2 range: -0.144505 0.975577

Import and Process CLY23 Field Data

Purpose: Import main phenotype data and calculate derived traits.

Approach: Read field book CSV, calculate flowering time metrics (DTS, DTA, ASI), parse donor identifiers, merge with LI-COR data.

Expected outcome: Complete CLY23 dataframe with all phenotypes and donor metadata.

planting_date <- mdy("4/4/2023")

CLY23 <- read.csv("CLY23_D4_FieldBook.csv") %>%
  rename(silking = "DTS", anthesis = "DTA", pedigree = "Female_genotype") %>%
  mutate(
    DTS = (mdy(silking) - planting_date) %>% as.integer(),
    DTA = (mdy(anthesis) - planting_date) %>% as.integer(),
    dSPAD = SPAD2 - SPAD1,
    donor_pop = gsub("-.*", "", pedigree, perl = TRUE),
    donor_line = gsub("_P.*?$", "", pedigree, perl = TRUE)
  ) %>%
  left_join(licor %>% dplyr::select(CLY23_D4, PhiPS2))

# Clean donor_line strings
CLY23$donor_line <- gsub("^.*?-", "", CLY23$donor_line, perl = TRUE)
CLY23$donor_line <- gsub("-+", "_", CLY23$donor_line, perl = TRUE)
CLY23$donor_line <- gsub("\\(.*", "", CLY23$donor_line, perl = TRUE)
CLY23$donor_line <- gsub("-", "_", CLY23$donor_line, perl = TRUE)
CLY23$donor_line <- gsub(" +$", "", CLY23$donor_line, perl = TRUE)

# Calculate derived traits
CLY23$ASI <- CLY23$DTS - CLY23$DTA
CLY23$LL <- CLY23$SL + CLY23$BL
CLY23$StPu <- as.factor(CLY23$StPu) %>% as.numeric() - 1
CLY23$StPi <- as.factor(CLY23$StPi) %>% as.numeric() - 1
CLY23$donor_line <- factor(CLY23$donor_line)

# Diagnostics
{
cat("CLY23 records:", nrow(CLY23), "\n")
cat("Donor populations:", length(unique(CLY23$donor_pop)), "\n")
cat("Donor lines:", nlevels(CLY23$donor_line), "\n")
cat("Available traits:", paste(colnames(CLY23), collapse = ", "), "\n")
}

## CLY23 records: 4446 
## Donor populations: 10 
## Donor lines: 83 
## Available traits: CLY23_D4, Rep, pedigree, Species, Notes, silking, anthesis, PH, EH, EN, BW, BL, SL, SPAD1, SPAD1_Date, SPAD2, SPAD2_Date, Leaf_Auringle, ST, StPu, StPi, Kinki, Prolif, NBR, LI_COR, LI_COR_Date, SCiO, SCiO_Date, DTS, DTA, dSPAD, donor_pop, donor_line, PhiPS2, ASI, LL

Define Trait Categories

Purpose: Organize phenotypic variables into functional categories for analysis and visualization.

Approach: Create vectors of trait names and corresponding category labels.

Expected outcome: Trait and category vectors for downstream analysis.

vars <- c("DTA", "DTS", "ASI", "PH", "EH", "EN", "BW", "BL", "SL", "LL", 
          "NBR", "SPAD1", "SPAD2", "dSPAD", "PhiPS2", "StPu", "StPi")

var_type <- c("FT", "FT", "FT", "Plant Size", "Ear", "Ear", "Plant Size", 
              "Plant Size", "Plant Size", "Plant Size", "Plant Size", 
              "Photosynthesis", "Photosynthesis", "Photosynthesis", "Photosynthesis", 
              "Highland", "Highland")

var_pal <- c("darkgoldenrod1", "darkred", "darkmagenta", "darkgreen", "burlywood")

# Diagnostics
{
cat("Traits defined:", length(vars), "\n")
cat("Categories:", paste(unique(var_type), collapse = ", "), "\n")
}

## Traits defined: 17 
## Categories: FT, Plant Size, Ear, Photosynthesis, Highland

Prepare PCA Data

Purpose: Aggregate phenotype data to genotype means for PCA.

Approach: Filter out controls, group by pedigree, calculate trait means across replicates.

Expected outcome: Genotype-level means suitable for multivariate analysis.

# Genotype categories for plotting
cats <- CLY23 %>%
  dplyr::filter(donor_pop != "Purple Check") %>%
  droplevels() %>%
  dplyr::select(donor_pop, donor_line, pedigree, Rep) %>%
  group_by(donor_pop, donor_line, pedigree) %>%
  summarize(n = n(), .groups = "drop")

# Genotype means for PCA
to_pca <- CLY23 %>%
  dplyr::select(donor_pop, donor_line, pedigree, Rep, all_of(vars)) %>%
  droplevels() %>%
  dplyr::filter(donor_pop != "Purple Check") %>%
  group_by(pedigree) %>%
  summarise(across(all_of(vars), ~ mean(.x, na.rm = TRUE)), .groups = "drop")

# Diagnostics
{
cat("Genotypes for PCA:", nrow(to_pca), "\n")
cat("Traits with complete data:", sum(colSums(is.na(to_pca[,-1])) == 0), "/", length(vars), "\n")
}

## Genotypes for PCA: 1329 
## Traits with complete data: 2 / 17

Principal Component Analysis

Purpose: Reduce dimensionality and identify major axes of phenotypic variation.

Approach: Run PCA on genotype means, visualize individual and variable contributions.

Expected outcome: PCA plots showing genotype clustering by donor population and trait loadings.

pca <- PCA(to_pca[, -1], graph = FALSE)

# Diagnostics
{
cat("=== PCA Variance Explained ===\n")
cat("PC1:", round(pca$eig[1, 2], 1), "%\n")
cat("PC2:", round(pca$eig[2, 2], 1), "%\n")
cat("Cumulative (PC1-2):", round(pca$eig[2, 3], 1), "%\n")
}

## === PCA Variance Explained ===
## PC1: 21.6 %
## PC2: 16.1 %
## Cumulative (PC1-2): 37.6 %

fviz_pca_ind(
  pca, 
  col.ind = cats$donor_pop, 
  geom = "point", 
  pointsize = 2,
  addEllipses = TRUE
) +
  guides(
    col = guide_legend(title = "Donor", override.aes = aes(label = "")),
    shape = "none", 
    fill = "none"
  ) +
  theme_classic2(base_size = 15)

PCA of genotypes colored by donor population

fviz_pca_var(pca, col.var = as.factor(var_type)) +
  scale_color_manual(values = var_pal) +
  guides(
    col = guide_legend(title = "Type", override.aes = aes(label = ""))
  ) +
  theme_classic2(base_size = 15)

PCA variable loadings by trait category

Trait Correlations

Purpose: Examine pairwise correlations among phenotypic traits.

Approach: Compute correlation matrix using pairwise complete observations.

Expected outcome: Correlation matrix revealing trait relationships.

cor_matrix <- cor(to_pca[, -1], use = "pairwise.complete.obs")
round(cor_matrix, 2)

##          DTA   DTS   ASI    PH    EH    EN    BW    BL    SL    LL   NBR SPAD1
## DTA     1.00  0.91  0.08  0.18  0.11 -0.18  0.10  0.15  0.17  0.17  0.06 -0.26
## DTS     0.91  1.00  0.44  0.16  0.07 -0.28  0.07  0.13  0.17  0.15  0.07 -0.27
## ASI     0.08  0.44  1.00 -0.03 -0.12 -0.30 -0.04 -0.03  0.06 -0.02 -0.02 -0.15
## PH      0.18  0.16 -0.03  1.00  0.74  0.02  0.28  0.42  0.54  0.48  0.34  0.17
## EH      0.11  0.07 -0.12  0.74  1.00  0.18  0.11  0.15  0.37  0.20  0.22  0.17
## EN     -0.18 -0.28 -0.30  0.02  0.18  1.00 -0.08 -0.10 -0.07 -0.11 -0.04  0.15
## BW      0.10  0.07 -0.04  0.28  0.11 -0.08  1.00  0.62  0.18  0.61  0.29  0.27
## BL      0.15  0.13 -0.03  0.42  0.15 -0.10  0.62  1.00  0.27  0.99  0.36  0.23
## SL      0.17  0.17  0.06  0.54  0.37 -0.07  0.18  0.27  1.00  0.42  0.18  0.05
## LL      0.17  0.15 -0.02  0.48  0.20 -0.11  0.61  0.99  0.42  1.00  0.37  0.22
## NBR     0.06  0.07 -0.02  0.34  0.22 -0.04  0.29  0.36  0.18  0.37  1.00  0.24
## SPAD1  -0.26 -0.27 -0.15  0.17  0.17  0.15  0.27  0.23  0.05  0.22  0.24  1.00
## SPAD2  -0.21 -0.27 -0.17  0.16  0.11  0.11  0.10  0.11  0.05  0.11  0.10  0.40
## dSPAD  -0.15 -0.20 -0.12  0.07  0.02  0.03 -0.01  0.00  0.02  0.00  0.00 -0.07
## PhiPS2 -0.06 -0.08 -0.08  0.12  0.07 -0.01  0.06  0.09  0.10  0.10  0.09  0.17
## StPu   -0.05 -0.05  0.00 -0.03  0.00 -0.02 -0.01 -0.03 -0.02 -0.03  0.06 -0.01
## StPi   -0.14 -0.14 -0.04  0.02  0.06  0.16 -0.08 -0.09 -0.03 -0.09  0.00 -0.06
##        SPAD2 dSPAD PhiPS2  StPu  StPi
## DTA    -0.21 -0.15  -0.06 -0.05 -0.14
## DTS    -0.27 -0.20  -0.08 -0.05 -0.14
## ASI    -0.17 -0.12  -0.08  0.00 -0.04
## PH      0.16  0.07   0.12 -0.03  0.02
## EH      0.11  0.02   0.07  0.00  0.06
## EN      0.11  0.03  -0.01 -0.02  0.16
## BW      0.10 -0.01   0.06 -0.01 -0.08
## BL      0.11  0.00   0.09 -0.03 -0.09
## SL      0.05  0.02   0.10 -0.02 -0.03
## LL      0.11  0.00   0.10 -0.03 -0.09
## NBR     0.10  0.00   0.09  0.06  0.00
## SPAD1   0.40 -0.07   0.17 -0.01 -0.06
## SPAD2   1.00  0.80   0.55 -0.01 -0.02
## dSPAD   0.80  1.00   0.46 -0.02  0.01
## PhiPS2  0.55  0.46   1.00 -0.02 -0.02
## StPu   -0.01 -0.02  -0.02  1.00  0.34
## StPi   -0.02  0.01  -0.02  0.34  1.00

Broad-Sense Heritability by Trait

Purpose: Estimate H² for each phenotypic trait.

Approach: For each trait, subset to genotypes with all 3 replicates, fit one-way ANOVA, partition variance.

Expected outcome: Heritability estimates for all measured traits.

her <- data.frame()

for (trait in vars) {
  # Remove NAs for this trait
  without_NA <- CLY23[!is.na(CLY23[, trait]), ]
  
  # Keep only genotypes with all 3 reps
  with_reps <- as.data.frame(table(line = without_NA$pedigree)) %>%
    filter(Freq == 3) %>%
    pull(line)
  
  to_H2 <- without_NA %>%
    filter(pedigree %in% with_reps)
  
  # Fit ANOVA model
  formula <- as.formula(paste0(trait, " ~ pedigree"))
  mod <- lm(data = to_H2, formula)
  mod_aov <- anova(mod)
  
  # Calculate H2
  her <- rbind(her, data.frame(
    trait = trait,
    H2 = mod_aov$`Sum Sq`[1] / sum(mod_aov$`Sum Sq`),
    n_geno = length(with_reps)
  ))
}

her$type <- var_type

# Diagnostics
{
cat("=== Heritability Estimates ===\n")
print(her %>% arrange(desc(H2)))
}

## === Heritability Estimates ===
##     trait        H2 n_geno           type
## 1      SL 0.7842433     12     Plant Size
## 2      PH 0.7212333   1010     Plant Size
## 3     DTA 0.7084150   1280             FT
## 4     DTS 0.6671166   1277             FT
## 5      EH 0.6432390   1008            Ear
## 6      LL 0.6175722     11     Plant Size
## 7     NBR 0.6163410     12     Plant Size
## 8      BL 0.6141236     11     Plant Size
## 9      EN 0.5899337   1008            Ear
## 10     BW 0.5273850     12     Plant Size
## 11  dSPAD 0.5096379     11 Photosynthesis
## 12  SPAD2 0.5087019    220 Photosynthesis
## 13    ASI 0.5063024   1274             FT
## 14   StPi 0.4515571   1307       Highland
## 15 PhiPS2 0.4506074    217 Photosynthesis
## 16  SPAD1 0.3891292     13 Photosynthesis
## 17   StPu 0.3884873   1307       Highland

her %>%
  mutate(trait = forcats::fct_reorder(trait, H2, .desc = TRUE)) %>%
  ggplot(aes(x = trait, y = H2, fill = type)) +
  geom_col() +
  ylab("Broad-Sense Heritability (H²)") +
  xlab("Trait") +
  scale_fill_manual(values = var_pal) +
  theme_classic2(base_size = 15) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "top"
  )

Broad-sense heritability by trait

Import CLY22 Brace Root Data

Purpose: Import previous year’s brace root data for cross-year comparison.

Approach: Read CLY22 field book, filter to J2Teo project, parse donor identifiers.

Expected outcome: Clean CLY22 dataframe for GxE analysis.

CLY22 <- read.csv("CLY-22-RRA-Field-Book - B6Field.csv") %>%
  dplyr::filter(Project == "J2Teo") %>%
  dplyr::mutate(
    donor_pop = gsub("-.*", "", Description, perl = TRUE),
    donor_line = gsub("_P.*?$", "", Description, perl = TRUE)
  )

CLY22$brace_root <- as.numeric(CLY22$brace_root)

# Clean donor_line strings
CLY22$donor_line <- gsub("^.*?-", "", CLY22$donor_line, perl = TRUE)
CLY22$donor_line <- gsub("-+", "_", CLY22$donor_line, perl = TRUE)
CLY22$donor_line <- gsub("\\(.*", "", CLY22$donor_line, perl = TRUE)
CLY22$donor_line <- gsub("-", "_", CLY22$donor_line, perl = TRUE)
CLY22$donor_line <- gsub(" +$", "", CLY22$donor_line, perl = TRUE)

# Diagnostics
{
cat("CLY22 J2Teo records:", nrow(CLY22), "\n")
cat("Brace root range:", range(CLY22$brace_root, na.rm = TRUE), "\n")
}

## CLY22 J2Teo records: 2408 
## Brace root range: 0 6

Combine Brace Root Data Across Years

Purpose: Merge brace root measurements from CLY22 and CLY23 for GxE analysis.

Approach: Standardize column names, stack datasets with field year indicator.

Expected outcome: Combined dataframe for cross-environment analysis.

BR <- rbind(
  CLY22 %>%
    dplyr::select(donor_line, donor_pop, NBR = brace_root) %>%
    dplyr::filter(donor_pop != "Purple Check") %>%
    mutate(field = "CLY22"),
  
  CLY23 %>%
    dplyr::filter(Rep == 3) %>%
    dplyr::select(donor_line, donor_pop, NBR) %>%
    dplyr::filter(donor_pop != "Purple Check") %>%
    mutate(field = "CLY23")
)

field_mean <- BR %>%
  group_by(field) %>%
  summarise(br_mean = mean(NBR, na.rm = TRUE), .groups = "drop")

# Diagnostics
{
cat("=== Combined Brace Root Data ===\n")
print(table(BR$field, BR$donor_pop))
cat("\nField means:\n")
print(field_mean)
}

## === Combined Brace Root Data ===
##        
##         B73 Bals Chal Dura Hueh Mesa Nabo Zdip Zlux
##   CLY22  17  570  195  359   84  338  212  424  197
##   CLY23 130  363  120  137   55  177  136  244  114
## 
## Field means:
## # A tibble: 2 × 2
##   field br_mean
##   <chr>   <dbl>
## 1 CLY22    1.88
## 2 CLY23    1.17

Brace Root Variation by Population and Year

Purpose: Visualize brace root number variation across donor populations and field years.

Approach: Faceted beeswarm plots with population means and field-level reference lines.

Expected outcome: Visual comparison of NBR distributions across populations and environments.

BR %>%
  ggplot(aes(x = donor_pop, y = NBR, group = donor_pop, col = donor_pop)) +
  geom_quasirandom() +
  geom_hline(
    data = field_mean,
    aes(yintercept = br_mean),
    colour = "black", linetype = 2
  ) +
  stat_summary(
    fun.data = mean_cl_normal,
    geom = "pointrange",
    col = "black"
  ) +
  facet_wrap(~field, ncol = 1) +
  scale_y_continuous(breaks = 0:6) +
  xlab("Donor Population") +
  ylab("Number of Brace Roots") +
  theme_classic2(base_size = 15) +
  theme(
    legend.position = "none",
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  )

Brace root number by donor population and field year

Brace Root Heritability by Year

Purpose: Estimate heritability of brace root number within each field year.

Approach: Fit one-way ANOVA (NBR ~ donor_line) separately for each year.

Expected outcome: Year-specific H² estimates for brace root number.

# CLY22 heritability
br22_aov <- lm(formula = NBR ~ donor_line, data = BR[BR$field == "CLY22", ])
s22 <- summary(br22_aov)

# CLY23 heritability
br23_aov <- lm(formula = NBR ~ donor_line, data = BR[BR$field == "CLY23", ])
s23 <- summary(br23_aov)

# Diagnostics
{
cat("=== Brace Root H² by Year ===\n")
cat("CLY22 H² (R²):", round(s22$r.squared, 3), "\n")
cat("CLY23 H² (R²):", round(s23$r.squared, 3), "\n")
}

## === Brace Root H² by Year ===
## CLY22 H² (R²): 0.257 
## CLY23 H² (R²): 0.064

Genotype × Environment Analysis

Purpose: Partition variance into genotype, environment, and GxE components for brace root.

Approach: Fit two-way ANOVA with interaction (NBR ~ donor_line * field), extract variance components.

Expected outcome: Proportion of variance explained by genotype, field year, and their interaction.

# Full GxE model
br_gxe_mod <- lm(formula = NBR ~ donor_line * field, data = BR)
gxe_aov <- anova(br_gxe_mod)

# Variance components
var_components <- data.frame(
  source = c("Genotype", "Field", "G×E", "Residual"),
  sum_sq = gxe_aov$`Sum Sq`,
  prop_var = gxe_aov$`Sum Sq` / sum(gxe_aov$`Sum Sq`)
)

# Summary table
gxe_summary <- data.frame(
  field = c("CLY22", "CLY23", "Combined"),
  H2_fam = c(
    s22$r.squared, 
    s23$r.squared, 
    gxe_aov$`Sum Sq`[1] / sum(gxe_aov$`Sum Sq`)
  ),
  field_year = c(NA, NA, gxe_aov$`Sum Sq`[2] / sum(gxe_aov$`Sum Sq`)),
  GxFY = c(NA, NA, gxe_aov$`Sum Sq`[3] / sum(gxe_aov$`Sum Sq`))
)

# Diagnostics
{
cat("=== Variance Components ===\n")
print(var_components)
cat("\n=== GxE Summary ===\n")
print(gxe_summary)
}

## === Variance Components ===
##     source    sum_sq   prop_var
## 1 Genotype  209.6703 0.08738041
## 2    Field  464.9371 0.19376325
## 3      G×E  150.4479 0.06269938
## 4 Residual 1574.4560 0.65615696
## 
## === GxE Summary ===
##      field     H2_fam field_year       GxFY
## 1    CLY22 0.25666810         NA         NA
## 2    CLY23 0.06435725         NA         NA
## 3 Combined 0.08738041  0.1937632 0.06269938

Population-Level GxE Analysis

Purpose: Estimate heritability and GxE within each donor population.

Approach: Split data by population, fit GxE models where both years are represented.

Expected outcome: Population-specific variance partitioning.

BR$field_donor <- interaction(factor(BR$field), factor(BR$donor_pop))

# Diagnostic: show donor_line counts per field × donor_pop combination
line_counts <- BR %>%
  group_by(field, donor_pop) %>%
  summarise(
    n_obs = n(),
    n_lines = n_distinct(donor_line),
    .groups = "drop"
  ) %>%
  arrange(n_lines)

# Diagnostics
{
cat("=== Donor Lines per Field × Population ===\n")
cat("(Combinations with n_lines < 2 will be skipped)\n\n")
print(as.data.frame(line_counts))
cat("\nSkipped combinations:", sum(line_counts$n_lines < 2), "\n")
}

## === Donor Lines per Field × Population ===
## (Combinations with n_lines < 2 will be skipped)
## 
##    field donor_pop n_obs n_lines
## 1  CLY22       B73    17       1
## 2  CLY23       B73   130       1
## 3  CLY22      Hueh    84       2
## 4  CLY22      Nabo   212       2
## 5  CLY23      Hueh    55       2
## 6  CLY23      Nabo   136       2
## 7  CLY22      Dura   359       3
## 8  CLY23      Dura   137       3
## 9  CLY22      Zdip   424       4
## 10 CLY23      Zdip   244       4
## 11 CLY22      Zlux   197       5
## 12 CLY23      Zlux   114       5
## 13 CLY22      Chal   195      13
## 14 CLY23      Chal   120      13
## 15 CLY22      Mesa   338      17
## 16 CLY23      Mesa   177      17
## 17 CLY22      Bals   570      35
## 18 CLY23      Bals   363      35
## 
## Skipped combinations: 2

# H2 by population within each field
br_list_field <- split(BR, f = BR$field_donor)

H2_pop <- lapply(names(br_list_field), FUN = function(x) {
  df <- br_list_field[[x]]
  n_lines <- length(unique(df$donor_line))
  if (nrow(df) < 3 || n_lines < 2) return(NULL)
  
  df$donor_line <- factor(df$donor_line)
  lm_aov <- lm(data = df, NBR ~ donor_line) %>% anova()
  
  data.frame(
    field = df$field[1],
    donor_pop = df$donor_pop[1],
    H2_fam = lm_aov$`Sum Sq`[1] / sum(lm_aov$`Sum Sq`)
  )
}) %>% 
  dplyr::bind_rows()

# GxE by population (across fields)
br_list_pop <- split(BR, f = BR$donor_pop)

# Diagnostic: show which populations have data in both fields
pop_field_summary <- BR %>%
  group_by(donor_pop) %>%
  summarise(
    n_fields = n_distinct(field),
    n_lines = n_distinct(donor_line),
    .groups = "drop"
  )

# Diagnostics
{
cat("\n=== Population Summary for GxE Analysis ===\n")
cat("(Populations with n_fields < 2 or n_lines < 2 will be skipped)\n\n")
print(as.data.frame(pop_field_summary))
}

## 
## === Population Summary for GxE Analysis ===
## (Populations with n_fields < 2 or n_lines < 2 will be skipped)
## 
##   donor_pop n_fields n_lines
## 1       B73        2       1
## 2      Bals        2      35
## 3      Chal        2      13
## 4      Dura        2       3
## 5      Hueh        2       2
## 6      Mesa        2      17
## 7      Nabo        2       2
## 8      Zdip        2       4
## 9      Zlux        2       5

H2_pop_gxe <- lapply(names(br_list_pop), FUN = function(x) {
  df <- br_list_pop[[x]]
  field_count <- table(df$donor_line, df$field) %>% as.data.frame()
  n_field <- length(unique(field_count[, 2]))
  n_lines <- length(unique(df$donor_line))
  
  if (n_field < 2 || n_lines < 2) return(NULL)
  
  df$donor_line <- factor(df$donor_line)
  lm_aov <- lm(data = df, NBR ~ donor_line * field) %>% anova()
  
  data.frame(
    donor_pop = df$donor_pop[1],
    H2_fam = lm_aov$`Sum Sq`[1] / sum(lm_aov$`Sum Sq`),
    GxF = lm_aov$`Sum Sq`[3] / sum(lm_aov$`Sum Sq`),
    field = lm_aov$`Sum Sq`[2] / sum(lm_aov$`Sum Sq`),
    residuals = lm_aov$`Sum Sq`[4] / sum(lm_aov$`Sum Sq`)
  )
}) %>% 
  dplyr::bind_rows()

# Diagnostics
{
cat("=== H² by Population and Field ===\n")
print(H2_pop)
cat("\n=== Population-Level GxE ===\n")
print(H2_pop_gxe)
}

## === H² by Population and Field ===
##    field donor_pop      H2_fam
## 1  CLY22      Bals 0.184208997
## 2  CLY23      Bals 0.091356100
## 3  CLY22      Chal 0.142280183
## 4  CLY23      Chal 0.134074433
## 5  CLY22      Dura 0.069233218
## 6  CLY23      Dura 0.011233303
## 7  CLY22      Hueh 0.014584654
## 8  CLY23      Hueh 0.003906250
## 9  CLY22      Mesa 0.366538899
## 10 CLY23      Mesa 0.131949990
## 11 CLY22      Nabo 0.005970679
## 12 CLY23      Nabo 0.005851729
## 13 CLY22      Zdip 0.120948596
## 14 CLY23      Zdip 0.012394894
## 15 CLY22      Zlux 0.201291366
## 16 CLY23      Zlux 0.006585232
## 
## === Population-Level GxE ===
##   donor_pop       H2_fam         GxF     field residuals
## 1      Bals 0.0741531279 0.046077328 0.1749622 0.7048073
## 2      Chal 0.0488590559 0.051865210 0.2526499 0.6466259
## 3      Dura 0.0214560436 0.015308706 0.2347739 0.7284613
## 4      Hueh 0.0053188460 0.004473805 0.4374807 0.5527267
## 5      Mesa 0.1575494230 0.127015582 0.0461216 0.6693134
## 6      Nabo 0.0008535911 0.003171211 0.4405808 0.5553944
## 7      Zdip 0.0275514370 0.039477679 0.1627404 0.7702305
## 8      Zlux 0.0594412897 0.039878067 0.2614525 0.6392282

Session Info

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.2
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggpubr_0.6.2     ggbeeswarm_0.7.3 factoextra_1.0.7 FactoMineR_2.13 
## [5] MASS_7.3-65      lubridate_1.9.4  ggplot2_4.0.1    dplyr_1.1.4     
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     vipor_0.4.7          farver_2.1.2        
##  [4] S7_0.2.1             fastmap_1.2.0        digest_0.6.39       
##  [7] rpart_4.1.24         timechange_0.3.0     estimability_1.5.1  
## [10] lifecycle_1.0.4      cluster_2.1.8.1      multcompView_0.1-10 
## [13] magrittr_2.0.4       compiler_4.5.2       rlang_1.1.6         
## [16] Hmisc_5.2-4          sass_0.4.10          tools_4.5.2         
## [19] utf8_1.2.6           yaml_2.3.12          data.table_1.18.0   
## [22] knitr_1.51           ggsignif_0.6.4       labeling_0.4.3      
## [25] htmlwidgets_1.6.4    scatterplot3d_0.3-44 RColorBrewer_1.1-3  
## [28] abind_1.4-8          withr_3.0.2          foreign_0.8-90      
## [31] purrr_1.2.0          nnet_7.3-20          grid_4.5.2          
## [34] xtable_1.8-4         colorspace_2.1-2     emmeans_2.0.1       
## [37] scales_1.4.0         dichromat_2.0-0.1    flashClust_1.01-2   
## [40] cli_3.6.5            mvtnorm_1.3-3        rmarkdown_2.30      
## [43] generics_0.1.4       otel_0.2.0           rstudioapi_0.17.1   
## [46] cachem_1.1.0         stringr_1.6.0        base64enc_0.1-3     
## [49] vctrs_0.6.5          jsonlite_2.0.0       carData_3.0-5       
## [52] car_3.1-3            rstatix_0.7.3        ggrepel_0.9.6       
## [55] Formula_1.2-5        htmlTable_2.4.3      beeswarm_0.4.0      
## [58] tidyr_1.3.2          jquerylib_0.1.4      glue_1.8.0          
## [61] DT_0.34.0            stringi_1.8.7        gtable_0.3.6        
## [64] tibble_3.3.0         pillar_1.11.1        htmltools_0.5.9     
## [67] R6_2.6.1             evaluate_1.0.5       lattice_0.22-7      
## [70] backports_1.5.0      leaps_3.2            broom_1.0.11        
## [73] bslib_0.9.0          Rcpp_1.1.0           coda_0.19-4.1       
## [76] gridExtra_2.3        checkmate_2.3.3      xfun_0.55           
## [79] forcats_1.0.1        pkgconfig_2.0.3