Why do I have frequency leakage in DFT after zero padding if frequency resolution is fine?
This phenomenon has nothing to do with spectral leakage. What you are observing is the effect of zero padding. Given a number of samples $N$, there is a maximum possible frequency resolution $\Delta f$ that can be achieved:

$$\Delta f=\frac{f_s}{N} $$

In your case $\Delta f$ is exactly $2\;\mathrm{Hz}$. If you zero-pad your signal, there is no extra information to retrieve - you will only decrease the frequency spacing.

In the example above, when you increase $N$ to $1000$, you get a frequency spacing of $1\;\mathrm{Hz}$.  All extra observed samples are merely an interpolation, done by the window function ($\mathrm{sinc}$ in your case). You will start observing side-lobes of the window spectrum.  Since you implicitly multiplied your signal by a rectangular window, this will result in the convolution of the spectrum of your signal (two Dirac's + DC) with the $\mathrm{sinc}$ function. 



Another way to look at it is to imagine that DFT is basically a filter bank, consisting of shifted $\mathrm{sinc}$ functions. Those are aligned in such a way, that peak of one is where zeroes of all the remaining ones are present. If you start looking in between those zeros, you will start taking those samples. Here is an example plot of such $\mathrm{sinc}$ filter bank.



Let's imagine that frequency corresponding to the blue filter is present. That will yield the amplitude in a corresponding bin. All the remaining frequencies are not present (orange and yellow), thus you multiply those $\mathrm{sinc}$'s by $0$ and get nothing in the bins. In the case of zero padding, that will be no longer the case. Samples of that blue $\mathrm{sinc}$ will fall in the intermediate bins and will be sinc-interpolated.



Here is what happens for $N=1000$ and $N=10000$:



And a zoomed part:



Things to notice:


For $N=500$, there is no leakage whatsoever. There are perfect spikes, representing each of your frequencies and DC offset.
We can also observe the FFT noise at the very bottom.
For $N=10000$ the shape of $\mathrm{sinc}$ function is clearly visible.




And obviously the code to reproduce the results:

Fs=1000; 
Ns=500;
Ns2=1000;
Ns3=10000;
t=0:1/Fs:(Ns-1)*1/Fs;
f1=10;
f2=400;
x=5+5*sin(2*pi*f1*t)+2*sin(2*pi*f2*t);

X1 = abs(fft(x))/length(x);
X2 = abs(fft(x, Ns2))/Ns;
X3 = abs(fft(x, Ns3))/Ns;

F1 = 0:Fs/Ns:Fs-Fs/Ns;
F2 = 0:Fs/Ns2:Fs-Fs/Ns2;
F3 = 0:Fs/Ns3:Fs-Fs/Ns3;

plot(F1, 20*log10(X1))
hold on
plot(F2, 20*log10(X2))
plot(F3, 20*log10(X3))
xlim([0, Fs/2])
grid on
legend({'N=500', 'N=1000', 'N=10000'})