I know that
scaling the action with a non-zero multiplicative constant, or
adding a total divergence term to the Lagrangian density
do not change the Euler-Lagrange equations, cf. e.g. this Phys.SE post.
Apart from such trivial modifications (1&2), is Einstein-Hilbert action the unique action whose variation gives Einstein's field equations? If not, is there any other action known which differs non-trivially from Einstein-Hilbert action and whose variation gives Einstein's equations?
In some sense, yes it is, and in others it is certainly not.
The fact that in the Einstein field equations: $G_{\mu \nu} = 8 \pi T_{\mu \nu}$, you have complete 'freedom' to define $T_{\mu \nu}$ however you would like (within some kind of exceptions), means you can allow for highly non-trivial curvature dependent terms to be coupled to the stress-energy. See this really interesting paper by Capozziello et al (http://arxiv.org/abs/grqc/0703067) that shows how you can 'bunch' the additional curvature terms for an $f(R)$ gravity into the 'effective' stress-energy tensor. In this way, the uniqueness of what you may call the 'Einstein' field equations does not hold.
However, once one moves to the vacuum the answer is suddently a definitive yes. The Einstein equations are generated by having an action that contains a geometric invariant quantity. If it did not then the principles of general covariance would no longer hold since arbitrary coordinate transformations would necessarily destroy the gauge symmetries. Now, that being said, one can consider the most general possible geometric invariant Lagrangian in this sense, it would look something like.. $\mathcal{L} = f(g_{\mu \nu} R^{\mu \nu}, R^{\mu \nu} R_{\mu \nu}, R^{\alpha \beta \gamma \delta} R_{\alpha \beta \gamma \delta}, C^{\alpha \beta \gamma \delta} C_{\alpha \beta \gamma \delta}, g^{\mu \nu} R^{\alpha \beta} R_{\alpha \beta \mu \nu}...)$ and the list of possibilities goes on. Since you want to replicate the Einstein equations in themselves you can immediatly throw away all but the linear terms in the above. Any of the other variants when you perform the calculus of variations gives you higher than second order derivatives.
Now, keeping only the linear terms you see that the action will necessarily be the Einstein-Hilbert one plus (maybe) some other curvature invariants. Since these will not vanish once you perform the variation, unless they either vanish identically or vanish on the boundary as a divergence term, your field equations will necessarily still be the Einstein-Hilbert ones provided you only keep $\mathcal{L} = R^{\sigma}_{\sigma}$.
A rigorous proof of this would be interesting. I am unaware of any in the literature.
The so-called spin content of GR with respect to the restricted Lorentz group is made of weight (also called spin) = 2 tensor fields which can be taken to be the linearized gravity field of Einstein and Fierz-Pauli (this is also known as the Pauli-Fierz field, after Pauli and Fierz's articles of 1939 in Helv. Phys. Acta and PRSL). In 2000 here http://arxiv.org/abs/hep-th/0007220 it is shown (as a particular case) using BRST (antifield) perturbation theory, that the cubic vertex of the single graviton theory is unique, hence the only way gravitons can self-couple is the E-H action + the cosmological term.
Regarding point 2 of your post: boundary terms in curved space-time are a most delicate subject and a careful analysis needs to be performed. You are invited to read the relevant section of Wald's GR test.
Another useful read would be the whole book and especially the chapter 8 of http://www.amazon.com/Tensors-Differential-Variational-Principles-Mathematics/dp/0486658406/ref=sr_1_1?ie=UTF8&qid=1392672388&sr=8-1&keywords=Lovelock+Tensors
You can also add terms that allow for the nonmetricity and torsion of the connection, and then constrain these terms to zero by explicit constraints.
