Non-normality of residuals in a negative binomial GLMM I am testing for the effect of treatment and fish length on the school size of fish. Treatment is a categorical variable with two levels (e.g., treatment A & treatment B). Fish length is continuous. There are five species of fish in this dataset and school size varies among fish species. Fish species is therefore treated as a random intercept. There are 4,622 data used to fit this model.  
That is:
glmer.nb(School.Size~Treatment+Fish.Length+(1|Fish.Species),data=dat).
The diagnostics of this model do not indicate a normal distribution of residuals. I think this is because the school size of fish can be vastly greater than predicted values, but not vastly smaller. This dataset is also zero truncated (i.e., there cannot be observations of a school size of zero), which I think is the cause for a small curve without residuals in the bottom left corner of the plot. I also think that this data may be peculiar because it's not "true" count data. Distributions of fish schools are almost certainly going to be over-dispersed relative to typical count data, and observers had to estimate (i.e., round) counts when fish schools were large.
My question is, given the non-normality (but not heteroskedasticity) of the residuals, can I move forward with this model? Or, how can I improve this model? For reference, I have tried to fit a GLMM with a gamma distribution on log-transforming school sizes and the residuals are similarly non-normal.
Thank you very much for any suggestions.

Pearson residuals vs. log-transformed fitted values of negative binomial GLMM. Numbers indicate values of school size (the response variable) and text color indicates species (levels of the random effect). Transparency is set to 25% such that overlaid points appear darker. A loess smoother (not shown) fit to these points shows homogenous variance around a mean of zero.
 A: Reginald, you have probably moved on by now, but for people that stumble across this post I would like to note that the DHARMa package (available from CRAN, see here) that I have created solves this problem, i.e. it will allow you to test if the residuals are compatible with the assumptions of nb (or any other distributions for that matter). 
From the package description:

The DHARMa package uses a simulation-based approach to create readily
  interpretable scaled residuals from fitted generalized linear mixed
  models. Currently supported are all 'merMod' classes from 'lme4'
  ('lmerMod', 'glmerMod'), 'glm' (including 'negbin' from 'MASS', but
  excluding quasi-distributions) and 'lm' model classes. Alternatively,
  externally created simulations, e.g. posterior predictive simulations
  from Bayesian software such as 'JAGS', 'STAN', or 'BUGS' can be
  processed as well. The resulting residuals are standardized to values
  between 0 and 1 and can be interpreted as intuitively as residuals
  from a linear regression. The package also provides a number of plot
  and test functions for typical model mispecification problem, such as
  over/underdispersion, zero-inflation, and spatial / temporal
  autocorrelation.

