  We propose to use high-redshift long $\gamma$-ray bursts (GRBs) as
cosmological tools to constrain the amount of primordial non-Gaussianity in the
density field. By using numerical, N-body, hydrodynamic, chemistry simulations
of different cosmological volumes with various Gaussian and non-Gaussian
models, we self-consistently relate the cosmic star formation rate density to
the corresponding GRB rate. Assuming that GRBs are fair tracers of cosmic star
formation, we find that positive local non-Gaussianities, described in terms of
the non-linear parameter, \fnl, might boost significantly the GRB rate at high
redshift, $z \gg 6$. Deviations with respect to the Gaussian case account for a
few orders of magnitude if \fnl$\sim 1000$, one order of magnitude for
\fnl$\sim 100$, and a factor of $\sim 2$ for \fnl$\sim 50$. These differences
are found only at large redshift, while at later times the rates tend to
converge. Furthermore, a comparison between our predictions and the observed
GRB data at $z > 6$ allows to exclude large negative \fnl, consistently with
previous works. Future detections of any long GRB at extremely high redshift
($z\sim 15-20$) could favor non-Gaussian scenarios with positive \fnl. More
stringent constraints require much larger high-$z$ GRB complete samples,
currently not available in literature. By distinguishing the contributions to
the GRB rate from the metal-poor population III regime, and the metal-enriched
population II-I regime, we conclude that the latter is a more solid tracer of
the underlying matter distribution, while the former is strongly dominated by
feedback mechanisms from the first, massive, short-lived stars, rather than by
possible non-Gaussian fluctuations. This holds quite independently of the
assumed population III initial mass function.
