Prepare some dummy data to investigate:

alist = RandomSample[Table[ToExpression["a" <> ToString[i]], i, 1, 10000]];
range = Range[10000];
Do[
 len = RandomInteger[3000, 10000];
 select = Sort[RandomSample[range, len]];
 sublists[j] = alist[[select]];
 , j, 1, 5000]

we have 5000 sublists of length between 3000 and 10000 elements "ai" in a particular order. I would like to have a function that merges all 5000 sublists such that all duplicates are discarded and any "ai" appear to the left of "aj" if they do so in any of the sublists. How can one do this in Mathematica most efficiently?

A tiny example of the above would be:

alist = a1,a3,a2,a5,a4;
sublists[1] = a1,a5;
sublists[2] = a3,a2,a5;
sublists[3] = a1,a2;

so that the function merge returns:

merge[Table[sublists[i],i,1,3]]

a1,a3,a2,a5

Note that merge was not given the actual alist.

EDIT:

Investigating the Experimental'ShortestSupersequence command, let us replace the definition of alist above by

alist = Table[a[i], i, 1, 10000];

This creates a list of strictly increasing a[i]. Generating the sublists from this as above, we get for example

seq = Fold[Experimental`ShortestSupersequence, sublists /@ Range[5000]];
seq//Length

10984

Repeating the steps seems to consistently return a seq that is longer than 10000, while

seq // DeleteDuplicates // Length

10000

To check that all sublists contain a[i] elements in strictly increasing order, we can do:

FreeQ[
 Table[
 tmp = sublists[i] /. a[x_] -> x;
 tmp = tmp[[2 ;;]] - tmp[[;; -2]];
 FreeQ[tmp/Abs[tmp], -1]
 , i, 1, 5000]
 , False]

True

The only way how some of the sublists might have some elements reversed in inconsistent order with alist is if we had ...,a[i],...,a[j],... with j<i somewhere, which is ruled out by the above test. So it seems that Experimental'ShortestSupersequence is buggy...

I think the following works but please give it a try with your data. I've tried quite a few random examples and it has worked for all of them.

First, construct a directed graph with edges defined by neighbors in the sublists:

allsublists = Array[sublists, 5000]; (* or whatever the max count is *)
G = Graph[Join @@ (BlockMap[Apply[Rule], #, 2, 1] & /@ allsublists)]

Next, look at the GraphDistanceMatrix and count how many times the symbol ∞ appears in each row:

infcount = Count[∞] /@ GraphDistanceMatrix[G]

Next, sort the graph vertices according to this infcount:

alist = SortBy[Transpose[VertexList[G], infcount], Last][[All, 1]]

All together in one function:

merge[L_] := 
 With[G = Graph[Join @@ (BlockMap[Apply[Rule], #, 2, 1] & /@ L)],
 SortBy[Transpose[VertexList[G], Count[∞] /@ GraphDistanceMatrix[G]], Last][[All, 1]]]

example

The given example is

allsublists = Array[sublists, 3]
(* a1, a5, a3, a2, a5, a1, a2 *)

Construct the directed neighbor graph:

G = Graph[Join @@ (BlockMap[Apply[Rule], #, 2, 1] & /@ allsublists),
 VertexLabels -> Automatic]

Have a look at the graph distance matrix:

VertexList[G]
(* a1, a5, a3, a2 *)
GraphDistanceMatrix[G]
(* 0, 1, ∞, 1,
 ∞, 0, ∞, ∞,
 ∞, 2, 0, 1,
 ∞, 1, ∞, 0 *)

count the infinities in each row, i.e., for each vertex we count how many other vertices are unreachable:

infcount = Count[∞] /@ GraphDistanceMatrix[G]
(* 1, 3, 1, 2 *)

sort the vertices by the number of infinities (i.e. by the number of unreachable other vertices):

alist = SortBy[Transpose[VertexList[G], infcount], Last][[All, 1]]
(* a1, a3, a2, a5 *)

Fold[Experimental`ShortestSupersequence, sublists /@ Range[3]]

a3, a1, a2, a5

Also

Fold[Experimental`ShortestSupersequence, RotateRight[sublists /@ Range[3]]]

a1, a3, a2, a5

In case this produces an output with duplicates we can use

DeleteDuplicates @ Fold[Experimental`ShortestSupersequence, sublists /@ Range[5000]]