(1) involves a multivariate Gaussian with potentially off-diagonal covariance matrix terms. You can see this from the previous section of Bishop, which gives the form of $q(w)$ as
$$
q(w) = \mathcal{N}(w \mid w_{\text{MAP}}, S_N),
$$
where $S_N$ can generally have off-diagonal terms. This is replaced by an integral over a single variable (as opposed to $D$ variables) with a univariate Gaussian. You see that even then, the integral isn't analytically tractable, which the author fixes by approximating the sigmoid function with a probit, $\Phi(\lambda a)$ and makes use of the identity,
$$
\int \Phi(\lambda a) \mathcal{N}(a \mid \mu \sigma^2) da = \Phi\left(\frac{\mu}{(\lambda^{-2} + \sigma^2)^{1/2}}\right),
$$
which wouldn't have been possible had he stuck with the original multivariate Gaussian with an unknown covariance. (I'm not sure if an analogous identity exists for a Gaussian with arbitrary covariance, though if you are able to reduce the expression to one that makes use of a simpler identity, then why not?)
As for the projection comment. Consider the following integral,
$$
S = \int_{-\infty}^\infty \int_{-\infty}^\infty f(a x + by) g(x,y) dx dy
$$
where $a$ and $b$ are constants. The integral is over the entire two-dimensional plane. Imagine, on the plane, the directional vector $(a,b)'$. You can define a line spanning the space along this direction. The inner product of any vector with $(a,b)'$ is a projection of that vector on this line. For example, suppose we have a vector $(c,d)'$ that's orthogonal to $(a,b)'$, so that $(c,d) (a,b)' = 0.$ Then $(x+c, y+d) (a,b)' = (x,y)(a,b)'.$ This means that the value of $f(ax + by)$ is determined entirely by the projection of $(x,y)'$ on the line defined by the directional vector $(a,b)'.$ Any change in a direction orthogonal to that line has no impact on the value of $f$.
We can then perceive this integral as an integral of $f$ over the line parameterized by a single variable $\lambda,$ given by a weighting function over that parameter, where the weight is given by $$p(\lambda) = \int \delta(\lambda - (a x + b y)) g(x,y) dx dy.$$
Consider what this integral is, for a given value of $\lambda.$ This is an integral of $g$ only over the line defined by $a x + b y = \lambda,$ which is perpendicular to the vector $(a,b)'.$ If $g$ is a two-dimensional probability distribution then you can think of this as marginalizing the distribution $g(x,y)$ over the direction orthogonal to $(a,b)'.$ If $g$ is a Gaussian, this marginalization is easy to do, and is done by Bishop in the passage you skipped over.
$S$ is an integral over $x$ and $y$, though you can re-imagine it as an integral over the line given by the direction $(a,b)'$, where at each incremental step, another integral is taken over the line orthogonal to $(a,b)'$ at that given point (so a rotation in the directions in which we integrate). This can be generalized in $D$ dimensions to an integral over a one-dimensional line, where at each increment a $D-1$ dimensional integral is taking place. Any function that depends only on a point's projection on the the line is a constant of the $D-1$ dimensional integral.