Better late than never… Just had an “off course!” revelation reading that thread (Low effective sample size for random effects only - #5 by mike-lawrence) about my model with three factors (\mu = \alpha_{site} + \beta_{species} + \gamma_{plot}): each combination of levels has only one observation! Hence, we have to make a choice and take one factor out. In my example case, we can assume plot (which can be seen as intra-site spatial variability) is the least important factor.
The issues with the two factors model (\mu = \alpha_{site} + \beta_{species}), which has 4 observations per combination of levels, remain to be solved: high Rhats (>1.15), low ESS (< 35 for \alpha and \beta), high correlation between sampled parameters, poor trace plots… with the default 4 chains with 1000 iterations.