We present a first QCD analysis of next-to-next-leading-order (NNLO) contributions of the spin-dependent parton distribution functions (PPDFs) in the nucleon and their uncertainties using the Jacobi polynomial approach. Having the NNLO contributions of the quark-quark and gluon-quark splitting functions in perturbative QCD [Nucl. Phys. B889, 351 (2014)], one can obtain the evolution of longitudinally polarized parton densities of hadrons up to NNLO accuracy of QCD. Very large sets of recent and up-to-date experimental data of spin structure functions of the proton g1p, neutron g1n, and deuteron g1d have been used in this analysis. The predictions for the NNLO calculations of the polarized parton distribution functions as well as the proton, neutron and deuteron polarized structure functions are compared with the corresponding results of the NLO approximation. We form a mutually consistent set of polarized PDFs due to the inclusion of the most available experimental data including the recently high-precision measurements from COMPASS16 experiments [Phys. Lett. B 753, 18 (2016)]. We have performed a careful estimation of the uncertainties using the most common and practical method, the Hessian method, for the polarized PDFs originating from the experimental errors. The proton, neutron and deuteron structure functions and also their first moments, Γp,n,d, are in good agreement with the experimental data at small and large momentum fractions of x. We will discuss how our knowledge of spin-dependence structure functions can improve at small and large values of x by the recent COMPASS16 measurements at CERN, the PHENIX and STAR measurements at RHIC, and at the future proposed colliders such as the Electron-Ion Collider.