A couple of points:

• It is better to fit the model using the adaptive Gaussian quadrature with a sufficient number of quadrature points, e.g., 10 or 15. This would provide a better approximation of the log-likelihood of the model. You could also give a try in the GLMMadaptive package that can also fit random slopes with the adaptive Gaussian quadrature.
• From your description, it seems that variable N.Level has a nonlinear effect on the log expected counts. You could account for that using splines or polynomials (the former are prefered). For example, you could first load the splines package, and the define your formula as Nodules ~ ns(N.Level, 3) * Species + (1 | Block).

A couple of points:

• It is better to fit the model using the adaptive Gaussian quadrature with a sufficient number of quadrature points, e.g., 10 or 15. This would provide a better approximation of the log-likelihood of the model. You could also give a try in the GLMMadaptive package that can also fit random slopes with the adaptive Gaussian quadrature.
• From your description, it seems that variable N.Level has a nonlinear effect on the log expected counts. You could account for that using splines or polynomials (the former are prefered). For example, you could first load the splines package, and then define your formula as Nodules ~ ns(N.Level, 3) * Species + (1 | Block).