Horizontal line arrays are often employed in underwater environments to estimate the direction of arrival (DOA) of a weak signal. Conventional beamforming (CB) is robust but has wide beamwidths and high-level sidelobes. High-resolution methods such as minimum-variance distortionless response (MVDR) and subspace-based MUSIC algorithm, produce low sidelobe levels and narrow beamwidths, but are sensitive to signal mismatch and require many snapshots and the knowledge of number of sources. In addition, heteroscedastic noise where the variance varies across observations and sensors due to nonstationary environments degrades the conventional methods significantly. This paper studies DOA in heteroscedastic noise (HN) environment, where the variance of noise is varied across the snapshots and the antennas. By treating the DOAs as random variables and the nuisance parameters of the noise variance different across the snapshots and the antennas, multi-snapshot variational line spectral estimation (MVALSE) dealing with heteroscedastic noise (MVHN) is proposed, which automatically estimates the noise variance, nuisance parameters of the prior distribution, number of sources, and provides the uncertain degrees of DOA estimates. When the noise variance only varies across the snapshots or the antennas, the variants of MVHN, i.e., MVHN-S and MVHN-A, can be naturally developed. Finally, substantial numerical experiments are conducted to illustrate the proposed algorithms' performance, including a real data set in DOA application.