Background

This exploratory data analysis examines differential gene expression in Arabidopsis thaliana abscission zone mutants. Abscission is the controlled detachment of plant organs such as flowers and fruits, occurring at specialized abscission zones (AZs). The development and activation of these zones are regulated by both hormonal signals, primarily auxin and ethylene, and key transcription factors.

Following fertilization, AZs undergo differentiation into two distinct cell layers: a separation layer that facilitates organ detachment through the secretion of hydrolytic enzymes, and a lignified layer that defines the break point. Understanding how hormonal signaling interacts with genetic regulation in this process is essential for identifying mechanisms controlling abscission timing.

The dataset used in this analysis was generated by a previous MSc student in our lab, Risham Osahan. Abscission zone tissues were collected from several Arabidopsis mutants, including bop1 bop2, ath1 knat2 knat6 (tripletale), bp, and hae hsl2. RNA sequencing was performed on these samples, and raw sequencing data (FASTQ files) were processed using DESeq2 to identify differentially expressed genes (done by Alicia H. from the Rowland lab).

The dataset consists of genes filtered based on standard thresholds (|log2 fold change| ≥ 1 and adjusted p-value ≤ 0.05). In this analysis, I focus on comparing the bp (early abscission) and ath1 knat2 knat6 (blocked abscission) mutants to determine whether they exhibit antagonistic gene expression patterns. In particular, I aim to identify candidate genes associated with auxin and ethylene signaling that may contribute to differences in abscission phenotypes.

Research Questions

To guide this analysis, the following questions are addressed:

1. How do gene expression patterns differ between bp and ath1 knat2 knat6 mutants relative to wild-type?
2. Among shared differentially expressed genes, do bp and ath1 knat2 knat6 mutants exhibit similar or opposing patterns of regulation?
3. What biological processes are enriched among genes showing opposing regulation between the two mutants?
4. Do hormone-related genes show patterns of regulation that are consistent with the contrasting abscission phenotypes?

Hygiene

Initial inspection of the dataset is performed to assess its structure, consistency, and completeness.

Viewing the first few rows of the bp dataset:

# Viewing the first few rows
head(bp_raw)

# A tibble: 6 × 13
  `Gene Name...1` `WT-1` `WT-2` `WT-3` `bp-1` `bp-2` `bp-3` baseMean `Log2 (FC)`
  <chr>            <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>    <dbl>       <dbl>
1 AT5G15360            0      0  0     1616.  1622.  1749.    332.         13.2 
2 AT5G41660            0      0  0      809.   810.   856.    165.         12.2 
3 AT5G17100            0      0  0      500.   472.   453.     95.1        11.4 
4 AT3G43340            0      0  0.809  136.   147.   159.     29.7         8.78
5 AT5G18880            0      0  0       57.7   45.5   88.8    12.9         8.53
6 AT4G00130            0      0  0       50.4   41.0   36.4     8.64        7.95
# ℹ 4 more variables: `P-adj.` <dbl>, `Up/Down` <chr>, `Gene Name...12` <chr>,
#   Description <chr>

Viewing the first few rows of the ath1 knat2 knat6 dataset:

head(tripletale_raw)

# A tibble: 6 × 13
  `Gene Name...1` `WT-1` `WT-2` `WT-3` `ath1knat2/6-1` `ath1knat2/6-2`
  <chr>            <dbl>  <dbl>  <dbl>           <dbl>           <dbl>
1 AT1G64790         0      0     0               78.4           47.5  
2 AT4G35690         0      3.62  1.62             8.97         256.   
3 AT1G62510         2.26   0     2.43            15.7          163.   
4 AT5G24200         0      0     1.62             1.12          71.3  
5 AT2G31930         1.13   0     0.809            7.84          56.1  
6 AT1G20860         0      0     0                6.72           0.951
# ℹ 7 more variables: `ath1knat2/6-3` <dbl>, baseMean <dbl>, `Log2(FC)` <dbl>,
#   `P-adj.` <dbl>, `Up/Down` <chr>, `Gene Name...12` <chr>, Description <chr>

Checking the dimensions of the bp dataset:

# Checking dimensions
dim(bp_raw)

[1] 1711   13

Checking the dimensions of the ath1 knat2 knat6 dataset:

dim(tripletale_raw)

[1] 1152   13

Checking the column names of the bp dataset:

# Checking column names
names(bp_raw)

 [1] "Gene Name...1"  "WT-1"           "WT-2"           "WT-3"          
 [5] "bp-1"           "bp-2"           "bp-3"           "baseMean"      
 [9] "Log2 (FC)"      "P-adj."         "Up/Down"        "Gene Name...12"
[13] "Description"   

Checking the column names of the ath1 knat2 knat6 dataset:

names(tripletale_raw)

 [1] "Gene Name...1"  "WT-1"           "WT-2"           "WT-3"          
 [5] "ath1knat2/6-1"  "ath1knat2/6-2"  "ath1knat2/6-3"  "baseMean"      
 [9] "Log2(FC)"       "P-adj."         "Up/Down"        "Gene Name...12"
[13] "Description"   

Checking for any missing values in the bp dataset:

# Checking for missing values
colSums(is.na(bp_raw))

 Gene Name...1           WT-1           WT-2           WT-3           bp-1 
             0              0              0              0              0 
          bp-2           bp-3       baseMean      Log2 (FC)         P-adj. 
             0              0              0              0              0 
       Up/Down Gene Name...12    Description 
             0            753            121 

Checking for any missing values in the ath1 knat2 knat6 dataset:

colSums(is.na(tripletale_raw))

 Gene Name...1           WT-1           WT-2           WT-3  ath1knat2/6-1 
             0              0              0              0              0 
 ath1knat2/6-2  ath1knat2/6-3       baseMean       Log2(FC)         P-adj. 
             0              0              0              0              0 
       Up/Down Gene Name...12    Description 
             0            474             94 

The inspection showed that both datasets contain 13 columns, including expression values, fold change, and gene annotations. The Excel file includes two columns labeled “Gene Name”, which are read in as “Gene Name…1” and “Gene Name…12”.

A small inconsistency was observed in the fold change column names: the bp dataset uses “Log2 (FC)”, whereas the tripletale dataset uses “Log2(FC)”. These differences should be standardized to ensure consistency between datasets.

Missing values are present in the gene annotation columns (gene name and description), which is expected since not all TAIR gene IDs have complete annotations. Importantly, there are no missing values in the main analysis columns (gene ID, expression values, log2 fold change, and adjusted p-values), indicating that the dataset is suitable for downstream analysis.

Tidying & Wrangling

Columns should be renamed so that both datasets use consistent variable names. This ensures that the two mutants can be compared directly and simplifies merging and downstream analysis

bp <- bp_raw %>%
  dplyr::rename(gene_id = `Gene Name...1`,
    wt_1 = `WT-1`,
    wt_2 = `WT-2`,
    wt_3 = `WT-3`,
    bp_1 = `bp-1`,
    bp_2 = `bp-2`,
    bp_3 = `bp-3`,
    log2fc = `Log2 (FC)`,
    padj = `P-adj.`,
    direction = `Up/Down`,
    gene_name = `Gene Name...12`,
    description = Description) %>%
  dplyr::mutate(gene_id = trimws(gene_id))

tripletale <- tripletale_raw %>%
  dplyr::rename(gene_id = `Gene Name...1`,
    wt_1 = `WT-1`,
    wt_2 = `WT-2`,
    wt_3 = `WT-3`,
    tripletale_1 = `ath1knat2/6-1`,
    tripletale_2 = `ath1knat2/6-2`,
    tripletale_3 = `ath1knat2/6-3`,
    log2fc = `Log2(FC)`,
    padj = `P-adj.`,
    direction = `Up/Down`,
    gene_name = `Gene Name...12`,
    description = Description) %>%
  dplyr::mutate(gene_id = trimws(gene_id))

After renaming, the datasets were checked for duplicated gene IDs and missing values to ensure data integrity before proceeding.

No duplicated gene IDs were detected for either mutant:

sum(duplicated(bp$gene_id))

[1] 0

sum(duplicated(tripletale$gene_id))

[1] 0

Missing values reassessed:

bp:

# Check missing values after tidying
colSums(is.na(bp))

    gene_id        wt_1        wt_2        wt_3        bp_1        bp_2 
          0           0           0           0           0           0 
       bp_3    baseMean      log2fc        padj   direction   gene_name 
          0           0           0           0           0         753 
description 
        121 

ath1 knat2 knat6:

colSums(is.na(tripletale))

     gene_id         wt_1         wt_2         wt_3 tripletale_1 tripletale_2 
           0            0            0            0            0            0 
tripletale_3     baseMean       log2fc         padj    direction    gene_name 
           0            0            0            0            0          474 
 description 
          94 

Missing values are limited to annotation columns and do not affect the main analysis variables. Taken together with the absence of duplicated gene IDs, this confirms that the datasets are suitable for downstream analysis.

# Wrangling  ----------------------------------------------

# I am now creating variables that will help classify genes as significantly upregulated, downregulated, or not significant. These variables will be used for plotting and for comparing shared genes between mutants.

deg_cutoff <- 0.05
lfc_cutoff <- 1

bp <- bp %>%
  dplyr::mutate(significance = dplyr::case_when(
      padj <= deg_cutoff & log2fc >= lfc_cutoff ~ "Upregulated",
      padj <= deg_cutoff & log2fc <= -lfc_cutoff ~ "Downregulated",
      TRUE ~ "Not significant"),
    neg_log10_padj = -log10(padj + 1e-10))

tripletale <- tripletale %>%
  dplyr::mutate(significance = dplyr::case_when(
      padj <= deg_cutoff & log2fc >= lfc_cutoff ~ "Upregulated",
      padj <= deg_cutoff & log2fc <= -lfc_cutoff ~ "Downregulated",
      TRUE ~ "Not significant"),
    neg_log10_padj = -log10(padj + 1e-10))

# Keeping only significantly differentially expressed genes
bp_sig <- bp %>%
  dplyr::filter(significance != "Not significant")

tripletale_sig <- tripletale %>%
  dplyr::filter(significance != "Not significant")

Genes were classified as significantly upregulated, downregulated, or not significant using combined thresholds of statistical significance (adjusted p-value ≤ 0.05) and biological relevance (|log2 fold change| ≥ 1). These categories will be used for visualization and for comparing shared genes between mutants. Adjusted p-values were transformed using −log10 to improve visualization in volcano plots.

Descriptive Statistics

# Descriptive statistics --------------------------------------------------

# Before moving into the main research questions, I want to summarize how many genes fall into each significance category for each mutant. This provides a simple overview of the dataset and also serves as a table for the report.

summary_table <- dplyr::bind_rows(bp %>%
    dplyr::count(significance) %>%
    dplyr::mutate(Mutant = "bp"),
  tripletale %>%
    dplyr::count(significance) %>%
    dplyr::mutate(Mutant = "ath1 knat2 knat6")) %>%
  dplyr::rename(Regulation = significance,
    `Number of Genes` = n) %>%
  tidyr::pivot_wider(names_from = Mutant,
    values_from = `Number of Genes`)

# Print clean table
summary_table %>%
  knitr::kable(caption = "Number of Upregulated and Downregulated Genes per Mutant",
    align = "c") %>%
  kableExtra::kable_styling(full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover"))

Both mutants show a wide range of gene expression changes. While the median log2 fold change is positive in both cases, bp exhibits more extreme values and a slightly negative mean, suggesting stronger downregulation compared to ath1 knat2 knat6. This suggests that bp may involve stronger repression or activation of regulatory pathways, consistent with its early abscission phenotype.

Q1. Exploring global expression patterns

To begin, overall gene expression patterns can be examined to assess how strongly each mutant differs from wild-type.

bp_volcano <- bp %>% 
  ggplot(aes(log2fc, neg_log10_padj, color = significance)) +
    geom_point(alpha = 0.6) +
  labs(title = "bp vs WT",
    x = "Log2 Fold Change (bp vs WT)",
    y = "-log10 Adjusted p-value") +
    theme_classic()
bp_volcano 

tripletale_volcano <- tripletale %>% 
  ggplot(aes(log2fc, neg_log10_padj, color = significance)) +
    geom_point(alpha = 0.6) +
  labs(title = "ath1 knat2 knat6 vs WT",
       x = "Log2 Fold Change (ath1knat2knat6 vs WT)",
       y = "-log10 Adjusted p-value") +
    theme_classic()
tripletale_volcano

Volcano plots are used to visualize both the magnitude of expression changes (log2 fold change) and their statistical significance. The x-axis represents log2 fold change, where genes on the right are upregulated and those on the left are downregulated relative to wild-type. The y-axis shows the negative log10 adjusted p-value, with higher values indicating greater statistical significance.
In both the bp and ath1 knat2 knat6 mutants, a large number of genes are significantly differentially expressed, indicating widespread transcriptional changes. Both mutants show a mixture of upregulated and downregulated genes; however, the bp mutant exhibits more extreme fold change values, particularly among downregulated genes, suggesting stronger repression of certain pathways. In contrast, ath1 knat2 knat6 shows more moderate expression changes, with a slight tendency toward upregulation.
Together, these results suggest that while both mutants broadly alter gene expression, bp induces stronger transcriptional shifts, whereas ath1 knat2 knat6 results in more moderate but still widespread changes.

Having established these global differences, the next step is to determine whether the same genes are affected in both mutants and how their regulation compares.

Q2. Shared genes and opposite regulation

To compare gene expression patterns between mutants, genes that are significantly differentially expressed in both bp and ath1 knat2 knat6 are identified. These shared genes are then classified based on whether they are regulated in the same direction or in opposite directions between the two mutants.

# Q2: Shared genes and opposite regulation --------------------------------

shared_genes <- dplyr::inner_join(
  bp_sig %>% dplyr::select(gene_id, log2fc_bp = log2fc),
  tripletale_sig %>% dplyr::select(gene_id, log2fc_tt = log2fc),
  by = "gene_id") %>%
  dplyr::mutate(regulation = dplyr::case_when(
      log2fc_bp > 0 & log2fc_tt > 0 ~ "Same",
      log2fc_bp < 0 & log2fc_tt < 0 ~ "Same",
      TRUE ~ "Opposite"))

# Venn diagram showing DEG overlap between mutants
print(ggVennDiagram(list("bp" = bp_sig$gene_id,
      "tt" = tripletale_sig$gene_id)) +
    scale_fill_gradient(low = "lightblue", high = "lavender") +
    ggtitle("Overlap of differentially expressed genes")+
      theme(legend.position = "none"))

The Venn diagram shows the overlap of significantly differentially expressed genes between the bp and ath1 knat2 knat6 (tt = tripletale) mutants. A total of 404 genes (16%) are shared between both mutants, indicating a subset of genes that are commonly affected despite their contrasting abscission phenotypes.

In contrast, 1,307 genes (53%) are unique to bp, while 748 genes (30%) are unique to ath1 knat2 knat6. This suggests that although there is some overlap in gene expression changes, each mutant also exhibits a large number of distinct transcriptional responses.

The greater number of unique genes in bp is consistent with its broader and more extreme changes in gene expression observed in the volcano plot.

# Scatterplot showing whether shared genes are regulated similarly or oppositely

sharedgenes_scatter <- shared_genes %>%   
 ggplot(aes(log2fc_bp, log2fc_tt, color = regulation)) +
    geom_point(size = 2) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_vline(xintercept = 0, linetype = "dashed") +
    scale_color_manual(values = c("Same" = "lightblue", "Opposite" = "coral")) +
  labs(title = "Comparison of shared gene expression between mutants",
    x = "Expression change in bp (log2 fold change)",
    y = "Expression change in ath1 knat2 knat6 (log2 fold change)")+
  geom_point(size = 2, alpha = 0.7)+
    theme_classic()

sharedgenes_scatter 

The scatterplot further compares expression changes for shared genes. Genes located in the top-right and bottom-left quadrants are regulated in the same direction, whereas genes in opposite quadrants show opposing regulation. Most shared genes follow similar expression patterns, although a smaller subset displays opposite regulation.

shared_genes %>%
  dplyr::count(regulation)

# A tibble: 2 × 2
  regulation     n
  <chr>      <int>
1 Opposite      46
2 Same         358

In total, 46 genes are oppositely regulated, while the majority show consistent regulation across mutants. These oppositely regulated genes are of particular interest, as they may contribute to the contrasting abscission phenotypes.

Based on this observation, the next step focuses specifically on these oppositely regulated genes to determine whether they are associated with particular biological functions.

Q3: Biological processes enriched among opposite genes

To investigate the functional significance of oppositely regulated genes, Gene Ontology (GO) enrichment analysis is performed. These genes are divided into two groups based on their direction of regulation: genes upregulated in bp and downregulated in ath1 knat2 knat6, and those showing the opposite pattern.

# Q3: Biological processes enriched among opposite genes ------------------

opposite_genes <- shared_genes %>%
  dplyr::filter(regulation == "Opposite")

bp_up_tt_down <- opposite_genes %>%
  dplyr::filter(log2fc_bp > 0 & log2fc_tt < 0)

bp_down_tt_up <- opposite_genes %>%
  dplyr::filter(log2fc_bp < 0 & log2fc_tt > 0)

gene_universe <- unique(c(bp$gene_id, tripletale$gene_id))
#Function for running the GO enrichment analysis
run_go <- function(df, title_text) {

  if (nrow(df) == 0) {
    message("No genes in ", title_text)
    return(NULL)
  }

  ego <- clusterProfiler::enrichGO(
    gene = unique(na.omit(df$gene_id)),
    universe = gene_universe,
    OrgDb = org.At.tair.db,
    keyType = "TAIR",
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.2
  )

  if (!is.null(ego) && nrow(as.data.frame(ego)) > 0) {
    print(
      clusterProfiler::dotplot(ego) +
        ggtitle(paste("GO enrichment:", title_text))
    )
    go_table <- as.data.frame(ego) %>%
      dplyr::select(
        Description,
        GeneRatio,
        Count,
        p.adjust
      ) %>%
      dplyr::arrange(p.adjust) %>%
      head(10) %>%
      dplyr::mutate(
        p.adjust = signif(p.adjust, 3)  # cleaner p-values
      ) %>%
      dplyr::rename(
        `GO Term` = Description,
        `Gene Ratio` = GeneRatio,
        `Gene Count` = Count,
        `Adjusted p-value` = p.adjust
      )
    go_table %>%
      knitr::kable(
        caption = paste("Top GO terms:", title_text),
        align = "c"
      ) %>%
      kableExtra::kable_styling(
        full_width = FALSE,
        position = "center",
        bootstrap_options = c("striped", "hover")
      )

  } else {
    message("No GO enrichment for ", title_text)
  }
}
run_go(bp_up_tt_down, "bp up / tripletale down")
run_go(bp_down_tt_up, "bp down / tripletale up")

# Inspecting genes in the smaller category (bp up / tripletale down). This is to identify candidate genes manually because GO enrichment did not work. 

#| echo: true
#| results: asis
#| code-fold: false

# Create a clean table of candidate genes (bp up / tripletale down)
bp_up_tt_down %>%
  dplyr::left_join(
    bp %>% dplyr::select(gene_id, gene_name, description),
    by = "gene_id"
  ) %>%
  dplyr::mutate(
    gene_label = dplyr::coalesce(gene_name, gene_id),  # use gene name if available
    description = stringr::str_trunc(description, 60)  # shorten long text
  ) %>%
  dplyr::select(
    `Gene` = gene_label,
    `Description` = description,
    `Log2FC (bp)` = log2fc_bp,
    `Log2FC (tripletale)` = log2fc_tt
  ) %>%
  dplyr::arrange(desc(abs(`Log2FC (bp)`))) %>%   # sort by strongest change
  head(10) %>%                                   # show top 10 genes
  knitr::kable(
    caption = "Candidate genes (bp up / tripletale down)",
    align = "c",
    digits = 2
  ) %>%
  kableExtra::kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover")
  )

A total of 46 oppositely regulated genes are identified, with a smaller subset in the bp up / ath1 knat2 knat6 down category and a larger subset in the reverse category.

GO enrichment is not detected for the smaller group, likely due to limited statistical power. To further explore this gene set, individual genes are examined manually to identify potential functional roles. A subset of these candidate genes is shown above. Although this group is too small to detect statistically significant enrichment, the presence of annotated genes with diverse functions suggests that it may still contribute to phenotypic differences between mutants.

In contrast, the larger group (bp down / ath1 knat2 knat6 up) shows significant enrichment for several metabolic and biosynthetic processes, including sulfur compound biosynthesis, glucosinolate metabolism, and amino acid biosynthesis. The enriched GO terms are summarized in both the dotplot and table above.

These results suggest that metabolic activity is relatively elevated in ath1 knat2 knat6 but reduced in bp. This functional divergence provides insight into how differences in gene regulation may contribute to their contrasting abscission phenotypes.

Given that abscission is strongly influenced by hormonal signaling, the next step examines whether hormone-related genes show similar or divergent patterns of regulation.

Q4: Hormone-related genes and phenotype interpretation

To determine whether hormone-related pathways differ between mutants, shared genes are screened for associations with auxin or ethylene signaling using gene names and functional descriptions. This approach is exploratory and intended to identify candidate genes rather than perform formal pathway analysis.

# Q4: Hormone-related genes and phenotype interpretation ------------------

shared_annotated <- dplyr::left_join(
  shared_genes,
  bp %>% dplyr::select(gene_id, gene_name, description),
  by = "gene_id")

shared_hormone <- shared_annotated %>%
  dplyr::mutate(text = tolower(paste(
      dplyr::coalesce(gene_name, ""),
      dplyr::coalesce(description, ""))),
    hormone = dplyr::case_when(
      stringr::str_detect(text, "auxin|iaa|arf|pin|yuc") ~ "Auxin",
      stringr::str_detect(text, "ethylene|ein|etr|erf|acs") ~ "Ethylene",
      TRUE ~ NA_character_ )) %>%
  dplyr::filter(!is.na(hormone))

hormone_summary <- shared_hormone %>%
  dplyr::count(hormone, regulation)

hormone_summary

# A tibble: 4 × 3
  hormone  regulation     n
  <chr>    <chr>      <int>
1 Auxin    Opposite       1
2 Auxin    Same          34
3 Ethylene Opposite      29
4 Ethylene Same         216

Most hormone-related genes are regulated in the same direction in both mutants. Specifically, 97% of auxin-related genes and 88% of ethylene-related genes show consistent expression patterns. This suggests that a large portion of the core hormone-related transcriptional response is shared between bp and ath1 knat2 knat6.

However, a smaller subset of genes displays opposite regulation, and this pattern is more pronounced among ethylene-related genes. Only one auxin-associated gene is oppositely regulated, compared with 29 ethylene-associated genes. This indicates that differences between the mutants are more strongly reflected in ethylene-responsive gene expression than in auxin-related gene expression.

hormone_summary_prop <- hormone_summary %>%
  dplyr::group_by(hormone) %>%
  dplyr::mutate(prop = n / sum(n)) %>%
  dplyr::ungroup()

hormoneprop <- hormone_summary_prop %>% 
  ggplot(aes(x = hormone, y = prop, fill = regulation)) +
    geom_col() +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)),
      position = position_stack(vjust = 0.5),
      color = "black") +
    labs(title = "Proportion of regulation patterns in hormone-related genes",
      x = "Hormone pathway",
      y = "Proportion") +
  scale_fill_manual(values = c(
  "Same" = "lightblue",
  "Opposite" = "pink"
)) +
    theme_classic()
hormoneprop

opp_hormone <- shared_hormone %>%
  dplyr::filter(regulation == "Opposite") %>%
  dplyr::mutate(gene_label = ifelse(is.na(gene_name), gene_id, gene_name))

if (nrow(opp_hormone) > 0) {
  print(
    ggplot(opp_hormone, aes(log2fc_bp, gene_label)) +
      geom_point(aes(color = hormone), size = 3) +
      geom_point(aes(x = log2fc_tt, color = hormone), size = 3) +
      geom_segment(
        aes(x = log2fc_bp, xend = log2fc_tt, y = gene_label, yend = gene_label),
        color = "grey") +
      scale_color_manual(values = c("Auxin" = "lightblue", "Ethylene" = "violet")) +
      labs(title = "Oppositely regulated hormone-related genes",
        x = "Log2 Fold Change",
        y = "Gene" ) +
      theme_classic() )
}

This plot highlights individual hormone-related genes that are oppositely regulated between mutants. Many ethylene-associated genes show contrasting expression patterns between bp and ath1 knat2 knat6, with several exhibiting higher expression in ath1 knat2 knat6 and lower expression in bp. In contrast, fewer auxin-related genes display this pattern, suggesting that differences in ethylene-responsive pathways may contribute more strongly to the contrasting abscission phenotypes.

At the gene level, these differences are reflected in both the direction and magnitude of expression changes, identifying candidate genes that may contribute to the contrasting abscission phenotypes.