$$\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*}$$