I'm sure there are better ways to do it but maybe this is an idea to get started.

Divide your "Distance" variable in bins, each containing a reasonable number of datapoints. At a glance, bins from -2 to +2 in 0.5 step could do [i.e. seq(-2, 2, by= 0.5)].

Your data table now should have columns: "response", "bin", "treatment", "genotype".

Then fit the ANOVA model with interaction between bin, treatment, genotype:

aov1 <- aov(data= dat, response ~ as.factor(bin) * treatment * genotype)
summary.lm(aov1)

This should pick up that bins around 0 in treatment 1, genotype 4 are different from the baseline.

You can then check for comparisons between bins, treatments and genotypes with:

TukeyHSD(aov1)

I'm sure there are better ways to do it but maybe this is an idea to get started.

Divide your "Distance" variable in bins, each containing a reasonable number of datapoints. At a glance, bins from -2 to +2 in 0.5 step could do [i.e. seq(-2, 2, by= 0.5)].

Your data table now should have columns: "response", "bin", "treatment", "genotype".

Then fit the ANOVA model with interaction between bin, treatment, genotype:

aov1 <- aov(data= dat, response ~ as.factor(bin) * treatment * genotype)
summary.lm(aov1)

This should pick up that bins around 0 in treatment 1, genotype 4 are different from the baseline.

You can then check for comparisons between bins, treatments and genotypes with:

TukeyHSD(aov1)

EDIT after whuber's comment: A simple improvement to the above solution may be to use the absolute distance from 0 for the binning, e.g. use abs(seq(-2, 2, by= 0.5)), since bins to the left and right of 0 are assumed to be equivalent with respect to the response. This will half the number of bins and increase power. Bins could either be treated as a nominal variable or as an ordinal variable to reflect that there is an increasing trend moving towards 0.