Liquid-crystal phase equilibria of Lennard-Jones chain fluids and the solubility of a Lennard-Jones gas in the coexisting phases are calculated from Monte Carlo simulations. Direct phase equilibria calculations are performed using an expanded formulation of the Gibbs ensemble. Monomer densities, order parameters, and equilibrium pressures are reported for the coexisting isotropic and nematic phases of: (1) linear Lennard-Jones chains, (2) a partially-flexible Lennard-Jones chain, and (3) a binary mixture of linear Lennard-Jones chains. The effect of chain length is determined by calculating the isotropic-nematic coexistence of linear Lennard-Jones chain fluids made of 8, 10, and 12 segments (8-, 10-, 12-mer). The effect of molecular flexibility on the isotropic-nematic equilibrium is studied for a Lennard-Jones 10-mer chain fluid with one freely-jointed segment at the end of the chain. An isotropic-nematic phase split and fractionation are reported for a binary mixture of linear 7-mer and 12-mer chains. Simulation results are compared with theoretical results as obtained from a recently developed analytical equation of state based on perturbation theory. Excellent agreement between theory and simulations is observed. The solubility of a monomer Lennard-Jones gas in the coexisting isotropic and nematic phases is estimated using the Widom test-particle insertion method. A linear relationship between solubility difference and density difference at isotropic-nematic coexistence is observed. It is shown that gas solubility is independent of the nematic ordering of the fluid, at constant temperature and density conditions.