Multi-step time series forecasting (TSF) is a crucial element to support tactical decisions (e.g., designing production or marketing plans several months in advance). While most TSF research addresses only single-point prediction, prediction intervals (PIs) are useful to reduce uncertainty related to important decision making variables. In this paper, we explore a large set of neural network methods for multi-step TSF and that directly optimize PIs. This includes multi-step adaptations of recently proposed PI methods, such as lower-upper bound estimation (LUBET), its ensemble extension (LUBEXT), a multi-objective evolutionary algorithm LUBE (MLUBET) and a two-phase learning multi-objective evolutionary algorithm (M2LUBET). We also explore two new ensemble variants for the evolutionary approaches based on two PI coverage-width split methods (radial slices and clustering), leading to the MLUBEXT, M2LUBEXT, MLUBEXT2 and M2LUBEXT2 methods. A robust comparison was held by considering the rolling window procedure, nine time series from several real-world domains and with different characteristics, two PI quality measures (coverage error and width) and the Wilcoxon statistic. Overall, the best results were achieved by the M2LUBET neuroevolution method, which requires a reasonable computational effort for time series with a few hundreds of observations.