In this paper, we carry out the next-to-leading-order (NLO) studies on $Z to Upsilon(1S)+g+g$ via the color-singlet (CS) $bbar{b}$ state. We find the newly calculated NLO QCD corrections to this process can significantly influence its leading-order (LO) results, and greatly improve the dependence on the renormalization scale. By including the considerable feeddown contributions, the branching ratio $mathcal{B}_{Z to Upsilon(1S)+g+g}$ is predicted to be $(0.56 sim 0.95)times 10^{-6}$, which can reach up to $19% sim 31%$ of the LO predictions given by the CS dominant process $Z to Upsilon(1S)+b+bar{b}$. Moreover, $Z to Upsilon(1S)+g+g$ also seriously affect the CS predictions on the $Upsilon(1S)$ energy distributions, especially when $z$ is relatively small. In summary, for the inclusive $Upsilon(1S)$ productions in $Z$ decay, besides $Z to Upsilon(1S)+b+bar{b}$, the gluon radiation process $Z to Upsilon(1S)+g+g$ can provide indispensable contributions as well.