Why do we integrate square of wave function over N-1 electron coordinates?
$\newcommand{\el}{_\mathrm{e}}
\newcommand{\dif}{\mathrm{d}}
\newcommand{\braket}[2]{\langle{#1}\vert{#2}\rangle}$Consider an $n$-electron system in an arbitrary electronic state represented by the electronic wave function $\psi\el(\vec{q}_{1}, \vec{q}_{1}, \dotsc, \vec{q}_{n})$, where $\vec{q}_{i}$ stands for joint spin-spatial coordinate of $i$-th electron, $\vec{q}_{i} = \{ m_{si}, \vec{r}_{i} \}$.

Integrating the square of the absolute value of the electronic wave function over the spatial coordinates of all the electrons except for electron-one and summing up over the spin coordinates of all of them without any exception gives the probability density of finding electron-one with arbitrary spin at $\vec{r}_{1}$ while having the remaining $n - 1$ electrons with arbitrary spins at arbitrary positions,
\begin{equation*}
    \rho_{\mathrm{e-1}}(\vec{r}_{1})
    =
    \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int
    \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2}
    \dif \vec{r}_{2} \dif \vec{r}_{3} \dotsc \dif \vec{r}_{n} \, .
\end{equation*}
And yes, if we "integrate" over all coordinates in this expression we'll just get an inner product of the wave function with itself,
$$
    \braket{\psi\el}{\psi\el}
    =
    \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int
    \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2}
    \dif \vec{r}_{1} \dif \vec{r}_{2} \dif \vec{r}_{3} \dotsc \dif \vec{r}_{n} \, ,
$$
which is usually equals to 1 since the wave function is usually normalized to unity.

Electrons are, of course, indistinguishable, so that their labeling is arbitrary and can be changed any time.
Thus, the expression above should work equally well for calculating the probability density of finding electron-two, electron-three, etc., or, in other words, it gives the probability density of finding one particular electron, no matter which, with arbitrary spin at $\vec{r}_{1}$ while having the remaining $n - 1$ electrons with arbitrary spins at arbitrary positions.

The electron density, defined as the probability density of finding any of the $n$ electrons with arbitrary spin at $\vec{r}_{1}$, is evidently $n$ times as large as the probability of finding one particular electron,
\begin{equation} \label{eq:electron_density}
    \rho(\vec{r}_{1})
    =
    n \sum_{m_{s1}} \sum_{m_{s2}} \dotsi \sum_{m_{sn}} \iint \dotsi \int
    \left| ψ\el(\vec{q}_{1}, \vec{q}_{2}, \dotsc, \vec{q}_{n}) \right|^{2}
    \dif \vec{r}_{2} \vec{r}_{3} \dotsb \dif \vec{r}_{n} \, .
\end{equation}
Note here that, for a wave function that is normalized to unity, if we "integrate" over all coordinates in this expression we'll just get the number of electrons $n$, which is one of the well-known properties of the electron density,
\begin{equation}
    \int \rho(\vec{r}_{1}) \dif \vec{r}_{1} = n \, .
\end{equation}