We measure the halo bispectrum covariance in a large set of N-body simulations and compare it with theoretical expectations. We find a large correlation among (even mildly) squeezed halo bispectrum configurations. A similarly large correlation can be found between squeezed triangles and the long-wavelength halo power spectrum. This shows that the diagonal Gaussian contribution fails to describe, even approximately, the full covariance in these cases. We compare our numerical estimate with a model that includes, in addition to the Gaussian one, only the non-Gaussian terms that are large for squeezed configurations. We find that accounting for these large terms in the modeling greatly improves the agreement of the full covariance with simulations. We apply these results to a simple Fisher matrix forecast, and find that constraints on primordial non-Gaussianity are degraded by a factor of ∼2 when a non-Gaussian covariance is assumed instead of the diagonal, Gaussian approximation.