  We present a variational study of pseudo-spin $1/2$ Bose gases in a harmonic
trap with weak 3D spin-orbit coupling of $\bmsigma\cdot\mathbf{p}$ type. This
spin-orbit coupling mixes states with different parities, which inspires us to
approximate the single particle state with the eigenstates of the total angular
momentum, i.e. superposition of harmonic $s$-wave and $p$-wave states. As the
time reversal symmetry is protected by two-body interaction, we set the
variational order parameter as the combination of two mutually time reversal
symmetric eigenstates of the total angular momentum. The variational results
essentially reproduce the 3D skyrmion-like ground state recently identified by
Kawakami {\it et al.}. We show that these skyrmion-like ground states emerging
in this model are primarily caused by $p$ wave spatial mode involving in the
variational order parameter that drives two spin components spatially
separated. We find the ground state of this system falls into two phases with
different density distribution symmetries depending on the relative magnitude
of intraspecies and interspecies interaction: Phase I has parity symmetric and
axisymmetric density distributions, while Phase II is featured with special
joint symmetries of discrete rotational and time reversal symmetry. With the
increasing interaction strength the transition occurs between two phases with
distinct density distributions, while the topological 3D skyrmion-like spin
texture is symmetry protected.
