Aliasing at $f_0 + kf_s$
You've made an error in this statement: &quot;I know by sampling in time domain we get shifted replicas of the original spectrum which is two delta functions at $f_0$ and $−f_0$ but how this corresponds to $f_s$ and aliasing at $f_0+kf_s$?&quot;
The two delta functions at $-f_0$ and $f_0$ are the result of taking the Fourier transform of a sinusoid $x(t)$:
$$ x(t) = \sin(2{\pi}f_0t)$$
$$\mathcal{F}(x(t)) = \frac{1}{2i} \big(\delta(f - f_0) + \delta(f + f_0) \big)$$
We haven't introduced sampling yet, this is simply the result of taking the Fourier transform. When we start sampling, then we introduce replicas of the signal's spectrum.
Let's assume we sample the $x(t)$ above at a frequency of $f_s$. We will have copies of that signal's spectrum at $f_0 + kf_s$ for all integer values of $k$. Assuming we met Nyquist, there's no problem here. The problem now comes when we have the sinusoid
$$ s(t) = \sin(2{\pi}(f_0 + f_s)t)$$
$$\mathcal{F}(s(t)) = \frac{1}{2i}\big( \delta(f - [f_0 + f_s]) + \delta(f + [f_0 + f_s]) \big)$$
When sampled, we have the spectrum copies centered at $(f_0 + f_s) + kf_s$ for all integer values of $k$. If we choose $k = -1$ then we have a copy at
$$ (f_0 + f_s) + kf_s|_{k = -1} = f_0$$
So if we had a combination of $x(t)$ and $s(t)$ and sampled it, a spectrum copy of $s(t)$ would step on $x(t)$'s copy at $f_0$. This is aliasing.