A three-dimensional turbulence simulation of the SOL and edge regions of a toroidally limited tokamak is carried out. The simulation couples self-consistently the drift-reduced two-fluid Braginskii equations to a kinetic equation for neutral atoms. A diagnostic neutral gas puff on the low-field side midplane is included and the impact of neutral density fluctuations on D alpha light emission investigated. We find that neutral density fluctuations affect the D alpha emission. In particular, at a radial distance from the gas puff smaller than the neutral mean free path, neutral density fluctuations are anti-correlated with plasma density, electron temperature, and D alpha fluctuations. It follows that the neutral fluctuations reduce the D alpha emission in most of the observed region and, therefore, have to be taken into account when interpreting the amplitude of the D alpha emission. On the other hand, higher order statistical moments (skewness, kurtosis) and turbulence characteristics (such as correlation length, or the autocorrelation time) are not significantly affected by the neutral fluctuations. At distances from the gas puff larger than the neutral mean free path, a non-local shadowing effect influences the neutral density fluctuations. There, the D alpha fluctuations are correlated with the neutral density fluctuations, and the high-order statistical moments and measurements of other turbulence properties are strongly affected by the neutral density fluctuations.