We consider a system of d non-linear stochastic fractional heat equations in spatial dimension 1 driven by multiplicative d-dimensional space-time white noise. We establish a sharp Gaussian-type upper bound on the two-point probability density function of (u(s, y), u(t, x)). From this result, we deduce optimal lower bounds on hitting probabilities of the process {u(t, x) : (t, x) is an element of [0, infinity[xR} in the non-Gaussian case, in terms of Newtonian capacity, which is as sharp as that in the Gaussian case. This also improves the result in Dalang et al. (2009) for systems of classical stochastic heat equations. We also establish upper bounds on hitting probabilities of the solution in terms of Hausdorff measure. (C) 2020 Elsevier B.V. All rights reserved.