Non-Gaussianity and nonlinearity have been shown to be ubiquitous characteristics of El Nino Southern Oscillation (ENSO) with implication on predictability, modelling, and assessment of extremes. These topics are investigated through the analysis of third-order statistics of El Nino 3.4 index in the period 1870-2018, namely bicovariance and bispectrum. Likewise, the spectral decomposition of variance, the bispectrum provides a spectral decomposition of skewness. Positive and negative bispectral contributions identify modes contributing respectively to La Ninas and El Ninos, mostly in the period range 2-6 years. The ENSO bispectrum also shows statistically significant features associated with nonlinearity. The analysis of bicovariance reveals a nonlinear correlation between the Boreal Spring and following Winter, coming from an asymmetry of the persistence of El Nino, contributing hence to a reduction of Spring Predictability Barrier. The positive skewness and main features of the ENSO bicovariance and bispectrum are shown to be well reproduced by fitting a bilinear stochastic model. This model shows improved forecasts, with respect to benchmark linear models, especially of the amplitude of extreme El Ninos. This study is relevant, particularly in a changing climate, to better characterize and predict ENSO extremes coming from non-Gaussianity and nonlinearity.