Why do we assume that Gauss' law of gravitation to remain same even when there are extra spatial dimensions? Every kid knows that the magnitude of the gravitational force between two point masses is inversely proportional to the square of the distance $r$ between them i.e., $F\sim r^{-2}$. But experimentally, due to finite precision of measurement, you'll never get an exact value -2 for the exponent of the distance. There are practical and aesthetic reason why we believe and preach that the exponent is -2 . In the scenario of extra spatial dimensions the force varies like $r^{-(2+N)}$. 
I presume that we derive $r^{-(2+N)}$ dependence of the gravitational force in presence of extra dimensions by assuming Gauss' law of gravitation holds even in its pristine form $\nabla\cdot\textbf{g}=-4\pi G\rho$ in presence of extra dimensions.
$\bullet$ If yes, why do we assume that the Gauss' law of gravitation to remain same even when there are extra spatial dimensions?
$\bullet$ If not, how do we derive $F\sim r^{-(2+N)}$ dependence?
 A: It's mainly because the "minimal" extension of GR to higher dimensions is simply obtained by writing down the Einstein-Hilbert action in $D$ dimensions and calling it a day:
$$
S_4 = \int \sqrt{-g} (R + \mathcal{L}_m) \, d^4x \quad \rightarrow \quad S_D = \int \sqrt{-g} (R + \mathcal{L}_m) \, d^Dx 
$$
($\mathcal{L}_m$ here is the matter Lagrangian).  If you do this, then the equations of motion are still $G_{ab} = 8 \pi T_{ab}$, even though we're working in $D$ dimensions.  Moreover, if you construct the equations of motion for "small" perturbations about a flat background, the derivation of the linearized Einstein equation $\Box \bar{h}_{\mu \nu} = 16 \pi T_{\mu \nu}$ proceeds exactly as in the 4-D case (where $\bar{h}_{\mu \nu}$ is the trace-reversed metric perturbation).  And if you further take the limit where time derivatives are negligible and the matter sources are moving slowly, you obtain the Poisson equation for $h_{00}$:
$$
\nabla^2 \bar{h}_{00} = -16 \pi \rho.
$$
We can then identify the Newtonian gravitational potential in this limit as $\phi = - \frac{1}{4} \bar{h}_{00}$, and $\mathbf{g} = - \nabla \phi$.  Gauss's Law of gravitation is then recovered.  None of this derivation relies on spacetime having four dimensions;  it carries over naturally to $D$ dimensions.
That said, there are some other extensions of general relativity that are allowed in dimensions higher than 4.  In particular, Lovelock's theorem says that in $D$ dimensions, there are a total of $\lfloor (D-1)/2 \rfloor$ terms you can put into your Lagrangian that lead to a version of Einstein's equation that is a second-order PDE in the metric.  (Higher-order PDEs often have problems due to so-called Ostrogradsky instabilities, which is why they're usually avoided.)  
In $D = 3, 4$, this means that there's only one possible Lagrangian you can write down, namely the Einstein-Hilbert action. But for $D > 4$, there could be new curvature combinations that would enter into the Lagrangian, each with its own coefficient.  The higher-dimensional Einstein-Hilbert action above is a special case of one of these;  you just need to set all of the coefficients of the new terms equal to zero.  Moreover, the linearized gravitational limit of Lovelock gravity would  still (I think) be Poisson's equation: all of the new curvature combinations are higher-order in the curvature, meaning they would contribute terms quadratic in $\bar{h}_{\mu \nu}$ to the equations of motion and would thus be negligible in the linearized limit.  
