How to design contrasts for a three way mixed effect model with interactions?

Question

I have an experiment which comprise a numerical dependent variable, say a feature such as growth_rate, and three independent factor variables describing where the samples were collected, i.e. locality, were the collected samples were growth, i.e. medium, and which taxonomical group they belong to, i.e. taxa. plus, there is variable that is adding some random noise which will be considered the random effect.

What I need to test is the combined effect of taxa and medium, interaction, while controlling for locality. Further complication is that I need to test all the possible combinations between taxa and medium.

To solve such problem I am thinking about a model like:

lme(growth_rate ~ medium*taxa + locality, block=random, data)

but then how to construct the matrix for the multcomp function glht?

I am reading the vignette for the multcomp package and I can understand the two way anova and how it works but I am unable to extend it to what i need. Specifically, i was looking at the Two-Way Anova part but I am still missing the how to add the locality variable.

I was also thinking about merging the medium and taxa variables into one, as mentioned in this thread but I am not sure how to manage the fact that i also have the locality variable

These are some data, to give an idea what i am dealing with. I randomised the whole table, so it's not the real data.

structure(list(locality = c("L1", "L1", "L1", "L1", "L1", "L1", 
"L1", "L1", "L1", "L1", "L1", "L1", "L1", "L1", "L1", "L1", "L1", 
"L1", "L1", "L1", "L1", "L1", "L1", "L1", "L1", "L2", "L2", "L2", 
"L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", 
"L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", "L2", 
"L2", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", 
"L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", "L3", 
"L3", "L3", "L3"), medium = c("M1", "M2", "M3", "M1", "M2", "M3", 
"M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", 
"M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", 
"M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", 
"M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", 
"M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", 
"M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", "M1", "M2", "M3", 
"M1", "M2", "M3"), random = c("rnd9", "rnd9", "rnd9", "rnd9", 
"rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", 
"rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", 
"rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd9", "rnd7", 
"rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd6", "rnd6", "rnd6", 
"rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", 
"rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", 
"rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", "rnd7", 
"rnd8", "rnd8", "rnd8", "rnd8", "rnd8", "rnd8", "rnd6", "rnd6", 
"rnd6", "rnd6", "rnd6", "rnd6", "rnd6", "rnd6", "rnd6"), taxa = c("g1", 
"g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", 
"g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", "g1", 
"g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g3", 
"g3", "g3", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", "g2", 
"g2", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", "g3", 
"g3", "g3", "g3", "g3", "g3", "g3", "g2", "g2", "g2", "g2", "g2", 
"g2", "g1", "g1", "g1", "g3", "g3", "g2", "g2"), growth_rate = c(7L, 
2L, 7L, 4L, 5L, 1L, 6L, 10L, 0L, 5L, 4L, 0L, 10L, 0L, 1L, 3L, 
8L, 8L, 0L, 0L, 0L, 5L, 0L, 6L, 5L, 3L, 10L, 1L, 7L, 5L, 0L, 
1L, 7L, 10L, 3L, 3L, 6L, 6L, 6L, 2L, 2L, 1L, 10L, 0L, 5L, 7L, 
1L, 2L, 8L, 5L, 9L, 1L, 4L, 10L, 0L, 4L, 3L, 3L, 5L, 7L, 3L, 
5L, 10L, 5L, 2L, 0L, 10L, 0L, 9L, 9L, 3L, 1L, 10L, 1L, 0L)), class = "data.frame", row.names = c(NA, 
-75L))