Why NIntegrate does not give correct answer for smaller values of input

The integral over negative z has a large positive value and the integral over the positive domain has a large nearly equal negative value. I managed to get high precision results by breaking the integral into two parts:

mu22[k_] :=
 NIntegrate[
   z*1/Sqrt[2 \[Pi]] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
    (Erf[z/Sqrt[2 (1 + (k*z)/3)]] + 
      Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]),
   {z, -3/k, 0}, WorkingPrecision -> 50] +
  NIntegrate[
   z*1/Sqrt[2 \[Pi]] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
    (Erf[z/Sqrt[2 (1 + (k*z)/3)]] + 
      Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]),
   {z, 0, Infinity}, WorkingPrecision -> 50]
Table[mu22[i/10], {i, 1, 9}]

{err, err, err, 0.6, 0.55939805, 0.557353,0.55498491455072,0.55231,0.54935762525154}

I'm not convinced this is a correct result, but it seems reasonable. The "exact" .6 is interesting..

As @george2079 says, the second term in the parenthesis integrates to zero:

Assuming[k > 0, 
  Integrate[z*1/Sqrt[2π] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
  (Exp[18/k^2]*Erf[-(z + 6/k)/Sqrt[2 (1 + (k*z)/3)]]), {z, -3/k, ∞}]]

0

So we only need to integrate the first term in the parenthesis, which is easy:

mu22[k_] := NIntegrate[z*1/Sqrt[2π] (3/(3 + z*k))^(3/2) Exp[(-3 z^2)/(2 (3 + k*z))]*
            (Erf[z/Sqrt[2 (1 + (k*z)/3)]]), {z, -3/k, ∞}]

Table[mu22[k], {k, 0.1, 1, 0.1}]

{0.563994, 0.563409, 0.562441, 0.561099, 0.559398, 0.557354, 0.554985, 0.552312, 0.549358, 0.546144}