Strictly speaking, I think you can sometimes have a p-value with a categorical variable with more than 2-levels, but you can't have a fold-change value.

So, I like your ANOVA example, but I'm don't think it is correct to claim that no method for calculating RNA-Seq p-values can handle a categorical value with more than 2 levels. For example, I know DESeq2 / edgeR / limma-voom can all handle continuous variable analysis (even though I was previously confused about this): so, I'm 100% certain that you could make a comparison if you converted your categories into a numeric score.

Use likelihood ratio tests in edgeR.