$\begin{eqnarray*} n^{p+1} \int_{0}^{1} dx\, e^{-nx} \ln (1+x^p) &=& n^{p+1} \int_0^n \frac{dz}{n} e^{-z} \log\left(1 + \left(\frac{z}{n}\right)^p \right) \\ &=& n^p \int_0^n dz\, e^{-z} \sum_{k=1}^\infty \frac{(-1)^{k+1}}{k} \left(\frac{z}{n}\right)^{p k} \\ &=& \sum_{k=1}^\infty \frac{(-1)^{k+1}}{k} \frac{1}{n^{p (k-1)}} \int_0^n dz\, e^{-z} z^{p k} \\ &=& \sum_{k=1}^\infty \frac{(-1)^{k+1}}{k} \frac{1}{n^{p (k-1)}} \gamma(p k + 1,n) \\ &=& \sum_{k=1}^\infty \frac{(-1)^{k+1}}{k} \frac{1}{n^{p (k-1)}} \left[(p k)! + O\left(e^{-n}n^{p k}\right)\right] \\ &=& p! - \frac{(2p)!}{2 n^p} + O(n^{-2p}) \end{eqnarray*}$ Above we use the asymptotic expansion for the lower incomplete gamma function, $\begin{eqnarray*} \gamma(s,n) &=& \int_0^n dt\, e^{-t} t^{s-1} \\ &=& \Gamma(s) - \int_n^\infty dt\, e^{-t} t^{s-1} \\ &=& \Gamma(s) - e^{-n}n^{s-1} + O(e^{-n}n^{s-2}) \hspace{5ex} (\textrm{integrate by parts}). \end{eqnarray*}$