In this paper, we first address the space-time decay properties for higher order derivatives of strong solutions to the Boussinesq system in the usual Sobolev space. The decay rates obtained here are optimal. The proof is based on a parabolic interpolation inequality, bootstrap argument and some weighted estimates. Secondly, we present a new solution integration formula for the Boussinesq system, which will be employed to establish the existence of strong solutions in scaling invariant function spaces. We further investigate the asymptotic profiles and decay properties of these strong solutions. Our results recover and extend the important results in Brandolese and Schonbek (Tran. A. M.S. Vol 364, No.10, 2012, 5057-5090).