Numerical Integration involves evaluating the integrand not only at the limits of integration but at points inbetween them as well. So if there's anywhere in the domain of the integration where a numerical evaluation could produce a slightly complex value due to floating point inaccuracy, you will see this behavior.

NIntegrate is a complex beast. It is both a numerical algorithm and a symbolic one. Unless you ask it not to, NIntegrate will symbolically analyze your input and find equivalent Integrals to give a quicker or better answer. There's no guarantee that the new integral won't introduce a trivial imaginary value.

Normally, you use the Chop function along with the result of NIntegrate to remove such small imaginary values.

Alternatively, you can pick a numerical method for NIntegrate to use with Method option that prevents it from doing anything which might introduce a complex value. This generally isn't worth the headache. I'd use Chop.

In my case the integrand cannot be complex in the limits on integration.

When i'm changing the integrand to Im[integrand]:

NIntegrate[Im[MyD[0.001, Ez, Ep, 10, 2, WLU, WRD, mLUPi, m1i, mRDOWNi]], {Ez, -((mRDOWNi/mLUPi)WRD + e0.001), -((mRDOWNi/mLUPi)(WRD + e0.001))}, {Ep, (-e0.001 - Ez), (mRDOWNi/(mLUPi - mRDOWNi))(Ez + WRD + e0.001)}, MinRecursion -> 90, MaxRecursion -> 100*]

the output is:

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

0.