We present skeleton studies of non-Gaussianity in the CMB temperature anisotropy observed in the WMAP5 data. The local skeleton is traced on the 2D sphere by cubic spline interpolation which leads to more accurate estimation of the intersection posit
ions between the skeleton and the secondary pixels than conventional linear interpolation. We demonstrate that the skeleton-based estimator of non-Gaussianity of the local type (f_NL) - the departure of the length distribution from the corresponding Gaussian expectation - yields an unbiased and sufficiently converged f_NL-likelihood. We analyse the skeleton statistics in the WMAP5 combined V- and W-band data outside the Galactic base-mask determined from the KQ75 sky-coverage. The results are consistent with Gaussian simulations of the the best-fitting cosmological model, but deviate from the previous results determined using the WMAP1 data. We show that it is unlikely that the improved skeleton tracing method, the omission of Q-band data, the modification of the foreground-template fitting method or the absence of 6 extended regions in the new mask contribute to such a deviation. However, the application of the Kp0 base-mask in data processing does improve the consistency with the WMAP1 results. The f_NL-likelihoods of the data are estimated at 9 different smoothing levels. It is unexpected that the best-fit values show positive correlation with the smoothing scales. Further investigation argues against a point-source or goodness-of-fit explanation but finds that about 30% of either Gaussian or f_NL samples having better goodness-of-fit than the WMAP5 show a similar correlation. We present the estimate f_NL=47.3+/-34.9 (1sigma error) determined from the first four smoothing angles and f_NL=76.8+/-43.1 for the combination of all nine. The former result may be overestimated at the 0.21sigma-level because of point sources.
