The following simple motivation for the formula was given in my undergrad class, which I hope I'm not misremembering:
One importantTypically the way the Laplace transform arises in applications is when solving initial-value problems by transforming differential equations to polynomial equations. This technique was used extensively in the early theory of electrodynamics (as the Heaviside calculus), and it was only eventually justified using Laplace transforms.
Suppose one were looking for an integral transform $\mathcal{L}[f] = \int_a^bK(s,x)f(x)dx$ with the property that for all $s,$ $$\mathcal{L}[Df](s) = s\cdot \mathcal{L}[f](s).$$ The reason for this specific functional equation is hopefully clear: such a transform would take a differential equation to the 'obvious' associated polynomial. As to why an integral transform in the first place, one can motivate that in much the same way as the Fourier transform: we are decomposing $f$ into a sum of projections onto a one-parameter system of functions (see e.g. chapter 18 of Zorich's Mathematical Analysis for how this is done.)
A natural domain to work with for an IVP would be the positive real axis, and we then have by the above equation: $$\mathcal{L}[Df]:=\int_0^\infty K(s,x)f'(x) dx = s\cdot \int_0^\infty K(s,x)f(x) dx = s\cdot \mathcal{L}[f](s)$$ Integrating by parts and solving will yield the kernel $K(s,x) = e^{-sx}$ (and a boundary term).