Simple supersymmetric grand unified models based on the gauge group SO(10) require - in addition to gauge and matter unification - the unification of t-b-tau Yukawa couplings. Owing to sparticle contributions to fermion self-energy diagrams, the Yukawa unification however only occurs for very special values of the soft SUSY breaking parameters. We perform a search using a Markov Chain Monte Carlo (MCMC) technique to investigate model parameters and sparticle mass spectra which occur in Yukawa-unified SUSY models, where we also require the relic density of neutralino dark matter to saturate the WMAP-measured abundance. For Yukawa unified models with mu > 0, the spectrum is characterizd by three mass scales: first and second generation scalars in the multi-TeV range, third generation scalars in the TeV range, and gauginos in the similar to 100 GeV range. Most solutions give far too high a relic abundance of neutralino dark matter. The dark matter discrepancy can be rectified by i). allowing for neutralino decay to axino plus photon, ii). imposing gaugino mass non-universality or iii). imposing generational nonuniversality. In addition, the MCMC approach finds a compromise solution where scalar masses are not too heavy, and where neutralino annihilation occurs via the light Higgs h resonance. By imposing weak scale Higgs soft term boundary conditions, we are also able to generate low mu, m A solutions with neutralino annihilation via a light A resonance, though these solutions seem to be excluded by CDF/D0 measurements of the B-s -> mu(+)mu(-) branching fraction. Based on the dual requirements of Yukawa coupling unification and dark matter relic density, we predict new physics signals at the LHC from pair production of 350 - 450 GeV gluinos. The events are characterized by very high b-jet multiplicity and a dilepton mass edge around m((chi) over tilde2)(0) - m((chi) over tilde1)(0) similar to 50 - 75GeV.