For the first question, the momentum $p = m\dot{\vec{x}}-e \vec{A}$ that is conjugate to the coordinate is also gauge-variant in exactly the way to make the canonical action gauge-invariant:
$$
\begin{aligned}
\vec{A} &\rightarrow \vec{A} + \frac{\hbar}{e}\nabla\Lambda\\
\vec{p} &\rightarrow \vec{p} - \hbar\nabla\Lambda
\end{aligned}
$$
(by the way, please take note I am using SI units. The original question does not.)
For the second question, in almost all cases, it is possible write down the transition amplitude right away. However, there are certain cases, for example the famous Duru-Kleinert path integral representation for the Coulomb problem, that required care to get.
For the third question, you can't substitute an equation of motion into the action. That would be nonsense as you have found out. (actually it looks like you already knew this, so I don't know why you asked about it anyway)
Fourth question: the reason it the answer differs in the two different ways of looking at it is because one of the ways is just wrong (see question 3).
Fifth question: If you're given the Hamiltonian, then the Lagrangian is defined/obtained by taking the Legendre transformation with respect to the momentum:
$$L = p\dot{q} - H\,,$$
and where $p$ is eliminated in favor of the conjugate velocity $\dot{q}$.
Hello handsome poster (why thank you kind stranger). The answer to your question it turns out is in Rovelli's "Quantum Gravity", at least insofar as the free scalar field is concerned. This is done in the following way. As you may recall (from Feynman and Hibbs), through various arguments about doing the path integral as a perturbation on the classical path, the path integral of a Gaussian lagrangian for a point particle
$$L = a(t) \dot x^2 + b(t) \dot x x + c(t) x^2 + d(t) \dot x + e(t) x + f(t)$$
will be
$$K(x_a, t_a; x_b, t_b) = e^{\frac{i}{\hbar} S_{cl} [x_b, x_a]} \int_0^0 \exp\{ \frac{i}{\hbar} \int_{t_a}^{t_b}[a(t) \dot y^2 + b(t) \dot y y + c(t) y^2]dt\}\mathcal Dy(t)$$
Where the second term will be some function only depending on the beginning and end time, as it does not depend on the position.
Through a rather tedious calculation, you can work out that the classical action for the harmonic oscillator is
$$S_{cl} = \frac{m\omega}{2\sin(\omega T)} ((x_a^2 + x_b^2) \cos (\omega T) - 2 x_a x_b)$$
with the remaining term of the kernel
$$\int_0^0 \exp\{ \frac{i}{\hbar} \int_{t_a}^{t_b} \frac{m}{2}[\dot y^2 - \omega y^2]dt\}\mathcal Dy(t)$$
As the function $y(t)$ goes from $0$ to $0$ in the interval $T = t_b - t_a$, it can be written as the Fourier series
$$y(t) = \sum_n a_n \sin(\frac{n\pi t}{T})$$
The action can thus be transformed into
$$\int_{t_a}^{t_b} \frac{m}{2}[\dot y^2 - \omega y^2]dt = \frac{mT}{4} \sum_n [(\frac{n\pi}{T})^2 - \omega^2] a_n^2$$
With a path integral being a simple infinite product of the individual gaussians for each $a_n$, which becomes (cf Feynman Hibbs for more details)
$$F(T) = \sqrt{\frac{m\omega}{2 \pi i \hbar \sin(\omega T)}}$$
Now the free scalar field can be decomposed as an infinite number of harmonic oscillators via some Fourier transform
\begin{eqnarray}
S &=& \int_{t_a}^{t_b} dt \int d^3x \frac{1}{2}[\partial_\mu\varphi \partial^\mu \varphi - m^2 \varphi^2]\\
&=& \int \frac{d^3k}{(2\pi)^3} \{\int_{t_a}^{t_b} dt \frac{1}{2}[\dot \varphi(k)^2 -(\vec k^2 + m^2)\vert \varphi(k) \vert^2]\}
\end{eqnarray}
(with some use of Parseval's theorem) which is just the action of the harmonic oscillator. Then, via some generous physicist magic, we can write
$$\exp[\frac{i}{\hbar}\int \frac{d^3k}{(2\pi)^3} S_{SHO}[\varphi(\vec k)]] \approx \prod_k \exp[\frac{i}{\hbar (2\pi)^3} S_{SHO}[\varphi(\vec k)]]$$
For every mode, the transition amplitude will be
$$\exp[\frac{i}{\hbar(2\pi)^3} \frac{\omega}{2\sin(\omega T)} ((\varphi_a^2(k) + \varphi_b^2(k)) \cos (\omega T) - 2 \bar \varphi_a(k) \varphi_b(k)) \sqrt{\frac{\omega}{2 \pi i \hbar \sin(\omega T)}}$$
Which will give us, once properly multiplied back,
$$\exp[\frac{i}{\hbar} \int \frac{d^3 k}{(2\pi)^3} \frac{\omega}{2\sin(\omega T)} ((\varphi_a^2(k) + \varphi_b^2(k)) \cos (\omega T) - 2 \bar \varphi_a(k) \varphi_b(k)) \prod_k [\sqrt{\frac{\omega}{2 \pi i \hbar \sin(\omega T)}}]$$
$\prod_k [\sqrt{\frac{\omega}{2 \pi i \hbar \sin(\omega T)}}]$ corresponds to a normalization factor, which Rovelli writes as being
$$\mathcal N \approx \prod_k [\sqrt{\frac{m\omega}\hbar}] \exp[-\frac{V}{2} \int \frac{d^3 k}{(2\pi)^3} \ln [\sin (\omega T)]]$$
which is apparently formally divergent (I assume due to the volume term $V$ - Rovelli does not seem to mention what it is tho), but I suppose this is similar to the Hamiltonian of the Hilbert space method also being divergent.
Best Answer
I guess you might be interested in this kind of an introduction,
http://www.math.sunysb.edu/~tony/whatsnew/column/feynman-1101/feynman1.html