NIntegrate samples values at points in the domain. When the domain is so wide, it can miss features. The fix is to break the domain up so NIntegrate will integrate over each section separately.
In[8]:= NIntegrate[f[x], {x, -1000, 1000}, AccuracyGoal -> 8]
Out[8]= 9.52777*10^-136
In[9]:= NIntegrate[f[x], {x, -1000, -10, 10, 1000}, AccuracyGoal -> 8]
Out[9]= 5.54866