glmm unbalanced nested random factors imagine you have the following data structure
                                     data
                                       |
                                       |
                   env1    env2    env3    env4      env5
                    |       |       |        |         |
    replicates   1,2,3    4,5,6    7,8       9         10

At each replicate you apply your nominal treatment A & B and measure your response for each replicate at under seqeuntial A and B treatments. So each replicate has a data point for A and B treatments.
Using lme4 a GLMM of this data might look like:
glmer(response ~ treatment + (1|env/replicate), family=poisson, data=data)

I would like to know what are the effects (repercussions for my model) of having only one replicate per environment in some of the cases?
 A: tl;dr I think a Poisson GLMM should work fine with this setup and only one replicate for some env levels.  Some of the individual latent variables will be confounded with each other, but the random-effect variance estimates should be OK (??). (Sorry, don't know if this counts as "answer drawing from credible and/or official sources" ...)
The relevant statistical model for the response of the $i$th replicate from the $j$th environment is
$$
Y_{ij} \sim \text{Poisson}(\mu_{ij}) \\
\mu_i \sim ({\mathbf X} \beta)_{ij} + \epsilon_i + \epsilon_j \\
\epsilon_i \sim \text{Normal}(0,\sigma^2_{\text{rep}}) \\
\epsilon_j \sim \text{Normal}(0,\sigma^2_{\text{env}}) \\
$$
(I could have written this out as ${\mathbf X} \beta + {\mathbf Z}b$, but for this simple case it's easier to write out the random effects explicitly.) Since there is a unique value of $ij$ for each observation, the replicate-level random effect corresponds to modeling overdispersion (i.e. a logNormal-Poisson model).
As discussed above, $\epsilon_i$ and $\epsilon_j$ will be confounded/perfectly correlated/jointly unidentifiable for any $j$ where there is only one replicate, but I don't think this will interfere with the fitting of the model. Since this is a technical/algorithmic point, I don't know of a way of checking other than by simulation.
If there were too large a proportion of cases with a very small number of replicates, that could make it practically difficult to jointly estimate $\sigma^2_{\text{rep}}$ and $\sigma^2_{\text{env}}$.
Here's a single simulation example (requires the development version of lme4) that says that at least in one reasonably best-case scenario, we can recover reasonable parameter estimates:
library(lme4)
packageVersion("lme4")  ## 1.1.2
set.seed(101)
nlevs <- 20
nreps <- rpois(nlevs,lambda=3)+1
sum(nreps==1)  ## this gives 2 cases with only 1 rep per treatment
d <- data.frame(env=factor(rep(seq(nlevs),nreps)),
               rep=factor(unlist(sapply(nreps,seq))))
d$trt <- factor(sample(c("control","trt"),size=nrow(d),replace=TRUE))
    d$y <- simulate(~trt+(1|env/rep),family=poisson,
        newdata=d,
        newparam=list(beta=c(1,3),theta=c(1,1)))[[1]]
(g1 <- glmer(y~trt+(1|env/rep),family=poisson,data=d))
## Random effects:
##  Groups  Name        Std.Dev.
##  rep:env (Intercept) 1.101   
##  env     (Intercept) 1.137   
## Number of obs: 79, groups: rep:env, 79; env, 20
## Fixed Effects:
## (Intercept)       trttrt  
##       1.137        2.690  

The Wald confidence intervals on the treatment effect do include the true value (3) in this case ...
