Heavy Neutral Leptons (HNLs) are strongly motivated by theory due to their capability of simultaneously explaining the observed phenomena of dark matter, neutrino oscillations and the baryon asymmetry of the Universe. The existence of such particles would affect the expansion history of the Universe and the synthesis of primordial abundances of light elements. In this work we review, revise and extend the phenomenology of HNLs during the Big Bang Nucleosynthesis (BBN) epoch for masses up to 1 GeV. This is of great importance, as BBN is able to provide complementary bounds to those from upcoming and proposed laboratory experiments. To this end we have developed a high-precision Boltzmann code that simulates BBN in the presence of HNLs and takes into account all relevant HNL decay channels, as well as subsequent interactions of decay products (thermalization and decay showers), dilution due to QCD phase transition, active neutrino oscillations and corrections to the weak reaction rates. We present robust bounds on the lifetime and mixing angles of HNLs for masses $3\,\mathrm{MeV}\leq m_N \leq 1\,\mathrm{GeV}$ and show that BBN is able to constrain HNL lifetimes down to $0.03 - 0.05$ s, depending on the mixing pattern. Moreover, combining our results with current experimental searches, we can exclude HNLs that mix purely with electron neutrinos up to ${\sim}$450 MeV and those that mix purely with muon neutrinos up to ${\sim}$360 MeV, in both cases for lifetimes up to at least a few tens of seconds. Finally, we compare the BBN constraints with those obtained from Cosmic Microwave Background observations and explore how our results will be improved by a number of upcoming and proposed laboratory experiments.