Is the shape of the wavefunction of a particle in a free space always a 3D gaussian? Is the shape of the wavefunction of a particle in a free space at any particular time (thus a particle that is not subject to any external influence including the higgs field which is reponsible for its mass) always a 3D gaussian?
 A: Definitely not always a Gaussian.
All it means for a particle to be a "free particle" is to not be confined in a potential. (i.e., V(x) = 0 in the Schrodinger Equation.)
Pick an initial state that your particle is in, now solve the differential equation
$$\begin{align}
 i\hbar\frac{\partial}{\partial t}\psi = &H \psi \\
 = &(\frac{p^2}{2m} + V(x))\psi \\
 = &\frac{p^2}{2m} \psi \\
 = &-\frac{\hbar^2}{2m} \nabla^2 \psi \\
\end{align}
$$
So if you crack open a numerical partial differential equation solver, or solve it analytically, you have your solution for the behavior of $\psi(x,y,z, t)$ as a function of 3D as a function of time.
If your initial condition isn't a gaussian, then your solution won't be either!
A: The wavefunction is a solution of the relevant wave equation without a potential , Dirac for fermions, Klein Gordon for bosons, etc.
Wave equations' simplest solutions are plane waves, which are not good for describing a particle in free space, because , as it represents the probability of finding the particle, a plane wave cannot give a local probability , the particle could be anywhere in fourspace.
So for representing  non interacting free particles on uses the wave packet solutions,

For calculating interactions particles one uses the Feynman diagrams of quantum field theory, where it is not necessary to use such a complication in the the wavefunction. Related answer of mine here and here .
