Non parametric 2 way anova In my dataset, there are two distinct species, each with three unique treatments. My objective is to investigate whether there are any significant differences among the various species and treatments. However, my data does not meet the normal distribution assumption for a two-way ANOVA. Is it possible to conduct a non-parametric two-way ANOVA instead? Additionally, can you advise me on how to create a ggplot2 visualization to represent the results?

The traditional nonparametric test for a two-way factorial design is the Scheirer–Ray–Hare test.  The only implementation I know of in R is in the rcompanion package, of which I am the author.  However, this test receives mixed reviews. In particular, in my experience, it's not always as likely to find a significant interaction as some alternative approaches.
A more contemporary approach is aligned ranks transformation anova. It has an excellent implementation in R in the ARTool package. It has functions, for example, for post-hoc tests and effect size.
You may want to consider other approaches, like ordinal regression or generalized linear model, depending on your particular data and situation.
In regard to presenting the results graphically, you might just look for "interaction plot" with the ggplot package.
