Different kinds of S-matrices? It seems to me that the notion of an "S-matrix" refers to several different objects
One construction you can find in the literature is allowing the coupling constant to adiabatically approach 0 in the asymptotic past and future. This means the incoming and outgoing states belong to the free theory and the S-matrix is an operator on the free theory state space (Fock space). It is the ratio between the time evolution operators of the interacting (but adiabatically "turned off") and the free theory over an infinite time span
This construction already yields potentially several objects since there can be several free field limits. For example if the model admits something like S-duality then both the g -> 0 and the 1/g -> 0 limits are free. It is then possible to consider 2 kinds of S-matrices: in one we adiabatically set g -> 0, in the other 1/g -> 0
But it seems that there also should be something like a "fully interacting" S-matrix. This S-matrix should be an operator on the space $Fock(H_{discrete})$ where $H_{discrete}$ is the subspace of the full state space corresponding to the discrete part of mass spectrum
Is it indeed the case there are several different S-matrices? If so, why is it "carefully hidden" in the literature? Is it possible to compute the "fully interacting" S-matrix in perturbation theory using the Bethe-Salpeter equation to extract the bound state spectrum?
 A: I  don't think the S-Matrix depends on the existence of a free limit of the theory. Rather, I think it depends on cluster decomposition, which identifies the (non-perturbative) asymptotic states as approximately non-intercating multi-particle states, whose relationship to a specific Lagrangian may be complicated, e.g. they may be solitons or bound states etc. There are two dimensional models where the S-matrices can be computed exactly using algebraic methods,  both for perturbative and solitonic states, which demonstrates the distinction nicely.
On the other  hand, if you want to calculate the S-matrix in perturbation theory using the LSZ reduction, you need a concrete identification of that Fock space and those asymptotic states. In perturbation theory, for perturbative states, this is the Hilbert space of the free theory.  Note that this may no be the complete Hilbert space of the theory, and there are known examples where the perturbative S-matrix is not unitary since there is some non-zero probability to creating non-perturbative states. This is also demonstrated very concretely in two-dimensional examples.
I vaguely recall Weinberg having a nice discussion of the general definition in his QFT course, I am assuming he covers this in the first volume of his QFT series.
A: There is only one S-matrix, but it is frequently introduced in sloppy ways. 
The S-matrix is a unitary matrix between two isomorphic Fock spaces whose 1-particle sector contains precisely one particle for each bound state of the system. It can be constructed by the usual adiabatic textbook procedure if and only if there are no bound states (which is a standard requirement for ordinary perturbation theory already for nonrelativistic QM without fields). 
See a somewhat cryptic remark in Volume 1 of Weinberg's QFT book p.110 (at the anchor of the ** footnote), who is the only textbook author I know of who mentions this and hints at how to do the more general case perturbatively, though not with enough detail to be really telling. (He says a little bit more on p.461/2, but again it is quite cryptic.)
In particular, in case of QCD, only the S-matrix where the asymptotic states are hadrons, makes physical (and hence nonperturbative) sense. But this is not tractable perturbatively anyway, as it is part of the unsolved infrared problem for QCD. 
The problem is carefully hidden from textbooks because nobody really knows how to treat bound states in QFT, and talk about ignorance isn't very suitable for textbooks. 
Weinberg treats bound states in Chapter 14, but only for an electron in an external field, which begs the real question. He barely mentions the oldest (and quite unreliable) method for bound states, the Bethe-Salpeter equation, which figures on p.560, where one can find the remark ''It must be said that the theory of relativistic effects and radiative corrections in bound states is not yet in entirely satisfactory state.'' - an euphemism for the fact that it is a real mess, and nobody knows how to treat it well.
On a more positive note, people work nowadays with 

QFT in the front form, giving Schwinger-Dyson equations, 

Exciting Baryons: now and in the future
http://arxiv.org/pdf/1109.3690v2.pdf
Dyson-Schwinger equations and their application to hadronic physics
http://arxiv.org/pdf/hep-ph/9403224.pdf
or with lattice gauge theory to get bound state information.

Magnetic Properties of the Nucleon
http://pos.sissa.it/archive/conferences/139/170/Lattice 2011_170.pdf
Heavy-Baryon Spectroscopy from Lattice QCD
http://arxiv.org/pdf/1002.4710.pdf

See also https://physics.stackexchange.com/a/407913/7924
A: In the literature which is old enough, the existence of free limits in $g\to\infty$ or elsewhere – if the moduli space is more complicated – is neglected because S-duality and the multiplicity of the free limits wasn't known.
However, it's not too big an omission because all the matrices you mention differ at most by phases that depend on the external particle masses; the interacting part of the information in the S-matrix is identical. The interactions really occur in a region of the spacetime where the coupling constant has a particular finite value. The adiabatic turning-on is just a way to define the S-matrix rigorously and its details don't really matter.
So up to some transformations that are largely trivial, there is just one S-matrix in each physical system.
A: 
Is it indeed the case there are several different S-matrices?

The S-matrix is approached in a quite a few different ways.  Here are
seven (or is it six?) examples.  I personally bless only the final
three, and only the last of the blessed three is useful in practice.
A.  In reslib.com/book/Relativistic_Quantum_Mechanics#88
Bjorken and Drell follow an old approach of Feynman's ("the propagator
approach") which they acknowledge is "less compelling" than using QFT.
(Understatement!)
B. In reslib.com/book/Quantum_field_theory__Lewis_H_Ryder#175
Ryder may be echoing "the propagator approach".  Whatever it is, it
differs from the approaches below, and he attributes it to unpublished
lectures by Veltman.
C.  In reslib.com/book/Introduction_to_the_Theory_of_Quantized_Fields#207
Bogolubov and Shirkov follow an old approach (attributed to
Stueckelberg) in which the operator S is a functional of a
spacetime-cutoff function g.  (Ick.)
In sections II.3 and II.4 of
reslib.com/book/Local_quantum_physics__fields__particles__algebras#90

Haag describes (or touches on) four approaches:
D.    S = lim exp(iH₀t₂) exp(-iH(t₂ - t₁)) exp(-iH₀t₁)
Or some such.  If you split the middle exponential into two factors,
you end up with
S = Ω₊ Ω₋ 

where Ω₊ and Ω₋ are "Moller operators".  Haag mentions this only for
comparison.  It is not suitable for QFT.  But I was taught it in my
first QFT class.
D.  Haag-Ruelle collision theory.  Now we are getting somewhere.  This
is what you would (or should) call a fully-interacting S-matrix.  We
start with a Wightman field theory with a field that creates
single-particle states.  [We also assume a mass gap.  But let's put
infrared issues aside.]  By appropriately applying fields to |0> and
taking limits, we obtain "in" and "out" states |α out> and |β in>.  α
and β are descriptive of the Fock states of a free field theory.
Indeed, there are free "in" and "out" fields defined on the Hilbert
space along with the interacting field.  The S-matrix is then defined
as:
S[α,β] = <α out | β in>

Good, but computationally intractable.
E.  Araki-Haag collision theory.  I won't go into the details.  But we
can now drop the requirement that we have a field that creates
single-particle states.  This is important to algebraic QFT-ers,
because fields that create single fermions are unobservable.  Again,
computationally intractable.
F.  LSZ collision theory.
Bingo!  The starting assumption is that the interacting field weakly
approaches the "in" and "out" fields (mentioned above) in the two
asymptotic limits.  Lehmann, Symanzik, and Zimmermann then show how to
obtain
S[α,β] = <α out | β in>

from Green functions, ie, from vacuum expectation values of
time-ordered products of fields.  And Green functions can be obtained
in perturbation theory using either the Gellman-Low expansion or
Feynman path integrals.  (I vote for the latter.  In Euclidean space,
followed by Wick rotation.)
To summarize: In practice, we (should) obtain Green functions from
Feynman diagrams (by one of two methods) and then obtain S-matrix
elements from the Green functions à la LSZ.  I feel safe in saying
that this is the approach to Feynman diagrams and the S-matrix that is
found in "better" textbooks.

If so, why is it "carefully hidden" in the literature?

Physics went through a difficult phase between 1930 and 1970.  Please
read the following quote by Jost:
reslib.com/book/PCT__spin_and_statistics_and_all_that#39

The recovery may be incomplete.

Is it possible to compute the "fully interacting" S-matrix in
  perturbation theory using the Bethe-Salpeter equation to extract the
  bound state spectrum?

What you're asking for, if I understand correctly, is some way to get
the infinite sums (that are probed by the B-S equation) into the Green
functions for composite fields that create bound states.  Then you
could, ideally, compute proton-proton cross sections.
I'm not aware of such a method.  (But I probably wouldn't be aware even
if it existed.)
