We compare the rapidity, $y$, and the beam energy, $\sqrt{s}$, behaviours of the cross section of the data for $D$ meson production in the forward direction that were measured by the LHCb collaboration. We describe the observed cross sections using NLO perturbative QCD, and choose the optimal factorization scale for the LO contribution which provides the resummation of the large double logarithms. We emphasize the inconsistency observed in the $y$ and $\sqrt{s}$ behaviours of the $D$ meson cross sections. The $y$ behaviour indicates a very {\it flat} $x$ dependence of the gluon PDF in the unexplored low $x$ region around $x\sim 10^{-5}$. However, to describe the $\sqrt{s}$ dependence of the data we need a steeper gluon PDF with decreasing $x$. Moreover, an even steeper behaviour is needed to provide an extrapolation which matches on to the well known gluons found in the global PDF analyses for $x\sim 10^{-3}$. The possible role of non-perturbative effects is briefly discussed.